Predator and n-classes-of-prey model incorporating extended Holling type Ⅱ functional response for n different prey species

1.   Introduction

To understand and analyze the dynamic relationships between preys and predators, thousands of mathematical models have been considered. Lotka and Voltera used mathematical tools to model the prey-predator interaction, and their model is called the Lotka-Voltera model. Their model has been developed by many researchers, as they take into account many factors like time delay, harvesting factor and spreading disease among populations ^{[1,2,3,4,5,6]}. The major element in predator prey interaction is the functional response for the predator, the number of consumed prey per predator per unit time. In ^[7], Van Leeuwen, Edwin, et al. studied the generalization of functional response for a predator that switches between multiple species. Tapan Kumar in ^[8] used Holling type Ⅱ in stability analysis of a prey-predator model incorporating a prey refuge. The most useful functional response is the Holling type Ⅱ functional response, which is characterized by decelerating intake rate ^[9]. Many authors used this type of functional response for modeling the dynamics of interactions between predator and prey. Molla and Hafizul studied the dynamics of a predator prey model with Holling Type Ⅱ functional response ^[10]. Eletterby in ^[11] discussed the global stability and persistence of a two prey and one predator model. Sirisubtawee S, Khansai N, in ^[12] considered the predator prey model with generalized Holling type Ⅳ functional response. In nature, there are many predator species that consume more than one species of prey. For example, lions usually predate a number of large land-based animals, such as antelopes, buffaloes, crocodiles, giraffes, pigs, zebra, wild dogs and wildebees. There are some works modeling one predator and multiple preys ^[13,14] which do not use Holling type Ⅱ functional response in their work. In this paper, first we extend the Holling type Ⅱ functional response to n-classes (species) of prey, and then we will also consider a predator-n-prey model with the extended Holling type Ⅱ functional response. We have arranged the paper as follows: In Section 2, we derive the extended Holling type Ⅱ functional response for n species of prey. In Section 3 a mathematical model for dynamics of interaction between one predator and its n classes of prey with the extended Holling type Ⅱ for n classes of prey and a simplex invariant attractor set for the proposed model is found. In Section 4, the permanence of the proposed model is proved under certain conditions. In Section 5, an ecological threshold parameter is established and local and global stabilities of the predator free equilibrium point are studied. In Section 6 the model is also solved numerically using the Runge-Kutta method.

2.   Holling type Ⅱ functional response to n-preys

In this section, we extend the Holing type functional response for more than one species of prey. Suppose a predator consumes n different classes of prey, and for $i = \mathrm{1, 2}, \dots, n$, assume the following:

ⅰ.      ${X}_{i}$ is the number of prey in the ${i}^{th}$ class.

ⅱ.      ${N}_{i}$ is the number of prey in the ${i}^{th}$ class predated by the predator at unit time $t$.

ⅲ.      ${T}_{search}$ is the time for searching for the prey at unit time $t$.

ⅳ.      ${T}_{i}$ is the handling time for eating one prey in the ${i}^{th}$ class.

Then ${N}_{i}$ is proportional to ${X}_{i}$ and ${T}_{search} = t-\sum _{j = 1}^{n}{N}_{j}{T}_{j}$

Therefore, ${N}_{i} = {\alpha }_{i}{X}_{i}\left(t-\sum _{j = 1}^{n}{N}_{j}{T}_{j}\right)$ where ${\alpha }_{i}$ is the constant of proportionality and can be called the predation rate of the prey in the ${i}^{th}$ class.

Now,

and

Therefore,

and hence $\frac{d{N}_{i}}{dt} = \frac{{\alpha }_{i}{X}_{i}}{1+\sum _{j = 1}^{n}{\alpha }_{j}{T}_{j}{X}_{j}}$ which is the extending Holling type Ⅱ functional response for n-class of prey.

3.   Model formulation

Suppose one class of predator consumes n different classes of its prey. Then, to model the dynamics of the interactions between the predator and its prey, we assume the following:

1) Population density of the ${i}^{th}$ class of prey is denoted by ${X}_{i},$ and the predator is denoted by $Y$.

2) Each class of prey ${X}_{i}$ grows logistically with carrying capacity ${K}_{i} > 0$ and an inherent growth rate ${r}_{i} > 0$.

The predator species consumes each class of the prey according to the extended Holling type Ⅱ functional response to n preys. However, in the absence of prey populations, the predator population grows exponentially.

Then, such interaction dynamics can be modeled mathematically as follows:

Here, the new parameters ${\theta }_{i}$ and ${\alpha }_{i}$ for $i = \mathrm{1, 2}, 3.., n$ are the conversion rate of ${X}_{i}$ and predation rate of ${X}_{i},$ respectively, where ${\theta }_{i} < {\alpha }_{i}\forall i,$ and $d$ represents the natural death rate of Y. Note that the above system contains a lot of parameters which make the analysis difficult. So, we can reduce the number of parameters by using the following dimensionless variables:

Accordingly, system (1) can be rewritten in the following non dimensional form:

where ${a}_{i} = \frac{{r}_{i}}{d}$, ${b}_{i} = {\alpha }_{i}{T}_{i}{K}_{i}$, ${p}_{i} = \frac{{\alpha }_{i}{K}_{i}}{d}$ and ${c}_{i} = \frac{{\theta }_{i}{K}_{i}}{d},$ for $i = \mathrm{1, 2}, 3, \dots, n$.

It is observed that the non-dimensional system contains $4n+1$ parameters, while the original system contains (6n+1) parameters.

In addition, in the following theorem, some properties of all solutions for the system that initiates ${R}_{+}^{n+1}$ are given.

Theorem 1. In system (2),

ⅰ. all solutions that initiate with non-negative condition are unique and non-negative;

ⅱ. the simplex set ${\mathrm{\Omega }}_{n} = \left\{\left({x}_{1}, {x}_{2}, {x}_{3}, \dots, {x}_{n}, y\right)\in {R}_{+}^{n+1};{x}_{i}\le 1andy+\sum _{j = 1}^{n}{x}_{j}\le n+\sum _{j = 1}^{n}{a}_{j}\right\}$ is positively invariant and an attractor.

Proof (ⅰ). The right hand side of each equation in system (2) is continuous and has partial derivatives on the space ${R}_{+}^{n+1}$. Therefore, they satisfy the Lipschitzian condition. Then, by uniqueness theorem, the solution is unique, and since the solution initiates at a non negative point, $\frac{d{x}_{i}}{dt} = 0$ at ${x}_{i} = 0;$ $i = \mathrm{1, 2}, \dots, n$ and $\frac{dy}{dt} = 0$ at $y = 0$. Thus, no compartment of the solution of system (2) can cross zero. Therefore, the solution is always non negative.

Proof (ⅱ). Suppose that $\left({x}_{1}\left(t\right), {x}_{2}\left(t\right), {x}_{3}\left(t\right), \dots, {x}_{n}\left(t\right), y\left(t\right)\right)$ is a solution of system (2), and initiate ${\mathrm{\Omega }}_{n}$.

For each $i = \mathrm{1, 2}, 3, \dots, n$,

Solving the above inequality, we obtain that

If ${x}_{i}\left(0\right)\in {\mathrm{\Omega }}_{n}$ then ${x}_{i}\left(t\right)\le 1$, $\forall t\ge 0$.

Let $w\left(t\right) = y\left(t\right)+\sum _{j = 1}^{n}{x}_{j}\left(t\right),$ and then

By solving the above inequality, we obtain that

If $y\left(0\right)\in {\mathrm{\Omega }}_{n}$ since ${x}_{i}\left(0\right)\in {\mathrm{\Omega }}_{n}$, $w\left(0\right)\le n+\sum _{j = 1}^{n}{a}_{j}$, and hence $w\left(t\right)\le n+\sum _{j = 1}^{n}{a}_{j},$ $\forall t\ge 0.$

Therefore, ${\mathrm{\Omega }}_{n}$ is positively invariant.

From (3) and (4), it is concluded that $\underset{t\to \infty }{\mathrm{lim}}Sup{x}_{i}\le 1$ and $\underset{t\to \infty }{\mathrm{lim}}Supw\le \sum _{j = 1}^{n}\left({a}_{j}+n\right).$

Thus ${\mathrm{\Omega }}_{n}$ is an attractor.

4.   Permanence

In ecology, determining the criteria which make the model permanent is important, because they imply that the population continues to exist. Therefore, in this section, we prove that system (2) is permanent under the conditions.

Definition. ^[2] System (2) is said to be permanent if there exist positive constants ${d}_{1}$ and ${d}_{2}$ such that

Theorem 2. For any initial values ${x}_{i}\left(0\right) < 1$ and $y\left(0\right) > 0,$ if the following conditions hold, then system (2) is permanent.

where ${R}_{i}$ are given in the proof.

Proof. From Theorem 1 we have that

Thus, $\mathrm{max}\{\underset{t\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}Sup{x}_{i}\left(t\right), \underset{t\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}Supy\left(t\right)\}\le n+\sum _{j = 1}^{n}{a}_{j}$.

If $t\to \infty$, then

Let us define ${M}_{i} = 1-\frac{{p}_{i}}{{a}_{i}}\left(n+\sum _{j = 1}^{n}{a}_{j}\right)$. We let inequality (5) is satisfied, ${M}_{i} > 0,$ and ${\in }_{i} > 0$.

By solving differential inequality (7), we obtain that

Therefore, there is ${N}_{i} > 0$ such that

Let ${H}_{i} = min\left\{{x}_{i}\left(t\right); t\in [0, {N}_{i}]\right\}$, and then ${H}_{i} = 0$ if and only if there exist $s\in [0, {N}_{i}$ ] such that ${x}_{i}\left(s\right) = 0.$

However, at ${x}_{i} = 0$, $\frac{d{x}_{i}}{dt} = 0$, and ${x}_{i}\left(t\right) = 0$, $\forall t > {N}_{i}.$ This is contrary to (8), and hence ${H}_{i}$ $\ne 0$.

Let ${R}_{i} = min\left\{\frac{{M}_{i}}{2}, {H}_{i}\right\}$, and then

In Theorem 1 we proved that if ${x}_{i}\left(0\right) < 0$, then

Hence, $\frac{dy}{dt}\ge \left(\frac{\sum _{j = 1}^{n}{c}_{j}{R}_{j}}{1+\sum _{j = 1}^{n}{b}_{j}}-1\right)y$ $\forall t\ge 0$, and $\frac{\sum _{j = 1}^{n}{c}_{j}{R}_{j}}{1+\sum _{j = 1}^{n}{b}_{j}}-1 > 0$ due to condition (6).

Therefore, if $\left(0\right) > 0$, then $y\left(t\right) > y\left(0\right),$ $\forall t > 0$.

Hence $\underset{t\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{i}\mathrm{n}\mathrm{f}y\left(t\right) > y\left(0\right)$. The proof is completed.

5.   Predator free equilibrium points and Ecological threshold parameter

In this section we discuss a threshold parameter ${R}_{0}$ which is known as the ecological basic reproduction number, and it may be thought about as the number of predators one predator gives rise to during its life ^[10]. Note that the system (2) has ${2}^{n}$ equilibrium points which do not contain a predator compartment. The points are ${E}_{f}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)$ or ${E}_{p}^{n} = \left({x}_{p1}, {x}_{p2}, {x}_{p3}, .., {x}_{pn}, 0\right),$ where ${x}_{pi} = 1$ or ${x}_{pi} = 0$. At least one of ${x}_{fi}$ is zero where each of the above equilibrium points is locally asymptotically stable.

Theorem 3. The equilibrium point ${E}_{p}^{n} = \left({x}_{p1}, {x}_{p2}, {x}_{p3}, .., {x}_{pn}, 0\right)$ is a saddle point, and the predator free equilibrium point ${E}_{f}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)$ is locally asymptotically stable if ${R}_{0} = \frac{\sum _{j = 1}^{n}{c}_{j}}{1+\sum _{j = 1}^{n}{b}_{j}} < 1$.

Proof. For any point in ${R}_{+}^{n+1}$ the variational matrix of the system (2) is

where

Eigenvalues of $V\left({E}_{p}^{n}\right)$ for $p = \mathrm{2, 3}, 4, \dots, {2}^{n}$ in the direction of compartment whose value in ${E}_{p}^{n}$ is positive and the eigenvalues in the directions of compartment $f$ the system whose value is in ${E}_{p}^{n}$, is negative. Therefore ${E}_{p}^{n} = \left({x}_{f1}, {x}_{f2}, {x}_{f3}, .., {x}_{fn}, 0\right)$ is a saddle point.

However, the eigenvalues of $J\left({E}_{pf}^{n}\right)$ are ${\lambda }_{i} =$ $-{a}_{i}$ for $i = \mathrm{1, 2}, 3, \dots, n$ and ${\lambda }_{n+1} = \frac{\sum _{j = 1}^{n}{c}_{j}}{1+\sum _{j = 1}^{n}{b}_{j}}-1$.

Consequently, the predator free equilibrium point ${E}_{f}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)$ is locally asymptotically stable if and only if ${R}_{n} < 1.$

Note that the corresponding equilibrium point of ${E}_{1}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)$ in the original system (1) is ${E}_{k}^{n} = \left({K}_{1}, {K}_{2}, {K}_{3}, \dots, {K}_{n}, 0\right),$ and ${R}_{n} = \frac{\sum _{j = 1}^{n}{\theta }_{j}{K}_{j}}{d\left(1+\sum _{j = 1}^{n}{\alpha }_{j}{T}_{j}{K}_{j}\right)}$, where $\frac{\sum _{j = 1}^{n}{\theta }_{j}{K}_{j}}{\left(1+\sum _{j = 1}^{n}{\alpha }_{j}{T}_{j}{K}_{j}\right)}$ is the growth rate of the new predator at ${E}_{k}^{n} = \left({K}_{1}, {K}_{2}, {K}_{3}, \dots, {K}_{n}, 0\right),$ and $\frac{1}{d}$ is the life span of the predator. Subsequently, their product gives the mean number of new born of the predator. Therefore, ${R}_{n}$ can be interpreted as the ecological basic reproduction number ^[14].

In the following theorem, global stability analysis of the predator prey free equilibrium point is given.

Theorem 4. ${E}_{f}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)$ is globally asymptotically stable if the following conditions hold:

Proof. Consider the function

It is easy to show that ${L}_{n}\left({x}_{1}, {x}_{2}, {x}_{3}, \dots {x}_{n}, y\right)\in {C}^{1}\left({R}_{+}^{n+1}, R\right)$ with ${L}_{n}\left({E}_{f}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)\right) = 0$ and ${L}_{n}\left({x}_{1}, {x}_{2}, {x}_{3}, \dots {x}_{n}, y\right) > 0,$ $\forall \left({x}_{1}, {x}_{2}, {x}_{3}, \dots {x}_{n}, y\right)\in {R}_{+}^{n+1}$.

Furthermore,

Conditions (11) and (12) guarantee that

Thus, using condition (12), it is obtained that $\frac{d{L}_{n}}{dt}$ is negative. Hence, ${L}_{n}$ is a Lyapunov function with respect to ${E}_{f}^{n} = \left(\mathrm{1, 1}, 1, \dots, \mathrm{1, 0}\right)$, so ${E}_{f}^{n}$ is globally asymptotically stable.

Note that the condition (11) guarantees that ${R}_{n} < 1$.

6.   Numerical simulation

In this section, we apply model (2) on the interactions between three prey species and one predator with the parameter values being

The model is also solved numerically by using the Runge-Kutta method.

Observe that free equilibrium point ${E}_{f}^{3} = \left(\mathrm{1, 1}, \mathrm{1, 0}\right)$ of the predator is globally asymptotically stable, as shown in Figure 1.

It is observed that ${R}_{0} = 0.45 < 1$ for the parameters given in $†$, which confirms the analytical result in Theorem 3. If we exchange the values of ${c}_{i}$ and ${b}_{i}$ and solve the model numerically, we see that ${E}_{f}^{3} = \left(\mathrm{1, 1}, \mathrm{1, 0}\right)$ is unstable, as shown in Figure 2.

After this exchanging the ecological basic reproduction number increased to ${R}_{0} =$1.0714$> 1$. Additionally, this numerical solution confirms the analytical results in Theorem 3.

7.   Conclusions and discussion

In this paper, interaction dynamics of a predator and all its prey species are modeled with extended Holling type Ⅱ functional response for more than one prey species. We also have shown that the proposed model is uniformly bounded for any parameter values and the permanence condition is provided. We have studied the stability for positive equilibrium points numerically, whereas it is difficult to find enough conditions for the analytical part for the positive equilibrium points. The ecological basic reproduction of the model is established by using the criteria which make the predator free equilibrium point to locally asymptotically stable. Local and global asymptotic stability of the predator free equilibrium points is studied. We also conclude that if the ecological basic reproduction number ${R}_{0} < 1,$ then the predator free equilibrium point is asymptotically stable, which leads to extinction of the predators. On the other hand, if ${R}_{0} > 1,$ then the predator free equilibrium point is unstable, which implies the persistence of the predators. We solved system (2) in the case $n = 3$ numerically by using the Runge-Kutta method. The analytical results are confirmed by the numerical simulations.

Conflict of interest

The authors declare that there are no competing interests.