# Three-Prey One-Predator Continuous Time Nonlinear System Model

Authors: Jingyuan Zhang; Yan Yang
Journal: Complexity
 (2020)

---

## Abstract

In this paper, we propose a new multiple-prey one-predator continuous time nonlinear system model, in which the number of teams of preys is equal to 3; namely, a continuous time three-prey one-predator model is put forward and studied. The fourth-order differential equation is established, in which the prey teams help each other. The equilibrium points and stability are analyzed. When not considering preys help each other, we study the global stability and persistence of the model without help terms. The simulation results of system solutions with help terms corresponding to locally asymptotically stable equilibrium points and without help terms corresponding to globally asymptotically stable equilibrium points are given.


---

## Body

## 1. Introduction

Nonlinear systems exist in many ways with various research directions, such as chaotic circuit [1–8], neural network [9–15], and image encryption [16–18]. Prey-predator systems also have nonlinear characteristics. In 1965, Holling proposed the functional responses of three different species to simulate predation. These functional responses describe how predators transform the captured prey into their own growth. According to the research, the prey-predator model can be divided into four categories: (i) the first category is the single-prey single-predator model; (ii) the second category is the single-prey multipredator model; (iii) the third category is the multiprey multipredator model; and (iv) the fourth category is the multiprey single-predator model. There are many research studies on the first type of the model. For example, Rosenzweig and MacArthur [19] proposed the Rosenzweig–MacArthur model (R-M model) or the predator-prey model with Holling type II functional response and gave graphical presentation and stability conditions of predator-prey interactions. Hsu and Huang [20] studied global stability of the R-M model. In [21], the conditions of internal equilibrium were given and the stability of internal equilibrium was analyzed. It also discussed certain critical situations, some of which cannot occur in the usual model. Wang and Jing [22] studied the global stability and the existence of limit cycles of a predator-prey system. The existence of at least two limit cycles was proved by qualitative analysis and Poincare–Bendixson theory. Zhang and Hou [23] investigated a nonautonomous ratio-dependent predator-prey model with exploited terms. By using the coincidence degree theory, it was proved that there are at least four positive periodic solutions. For the second category, Llibre and Xiao [24] studied the competitive exclusion of two predators and one prey. In the absence of predators, the growth rate of the prey was logistic or linear. The boundedness of the solution was proved by studying the stability of the equilibrium point of the differential equation at infinity. Through the theoretical analysis of the equations, the necessary and sufficient conditions for the existence of the exclusion principle were obtained, and the global dynamic behaviors of the three species in the first octaves were given. Wang et al. [25] studied a class of n-dimensional consumer-resource systems, in which a group of consumers compete for the same resource, and each consumer is mutually beneficial with the resource. Wang and Wu [26] considered the global dynamics of an n-dimensional Lotka–Volterra system in which (n-1) predator species compete for a single prey species. By using the theory of dynamic systems, the global dynamic behaviors of n species in the first octant were given, and sufficient and necessary conditions for the existence of the competitive exclusion principle were proved, thus extending Volterra’s principle and the work of Llibre and Xiao [24]. For the third category, Elettreby and El-Metwally [27] gave models in which two teams of predators interact with two teams of preys. The teams in each group (predators or preys) help each other. Three different versions of the multiteam predator-prey model were proposed. In addition, the equilibrium solutions, the conditions of their local asymptotic stability, and the global stability of the solution of one of the models were studied. For the fourth category model, Elettreby [28] proposed a two-prey one-predator model, in which the prey teams help each other. Its local stability was studied. In the absence of predators, there was no help between the prey teams. The global stability and persistence of the model without help were studied. However, for the continuous time multiple-prey one-predator model, so far, only the above two-prey one-predator model has been proposed. No continuous time multiple-prey one-predator model in which the number of teams of preys is greater than 2 has been proposed. But in practice, in many cases, the number of teams of preys is greater than 2. Therefore, in this paper, we propose a new multiple-prey one-predator model, in which the number of teams of preys is equal to 3; namely, a continuous time three-prey one-predator model is put forward and studied. The fourth-order differential equation is established. The equilibrium points and stability are analyzed.


## 2. The Model

In this section, we propose a system consisting of three teams of preys with densitiesx, y, and z, respectively, interacting with one team of predator with densities w. The assumptions of this model are as follows:(1)
In the absence of any predation, each team of preys grows logistically, that is,ax1−x, by1−y, and cz1−z.(2)
The effect of the predation is to reduce the prey growth rate by a term proportional to the prey and predator populations, that is, the −xw, −yw, and −zw terms.(3)
The teams of preys help each other against the predator, that is, axyzw term exists.(4)
In the absence of any prey for sustenance, the predator’s death rate results in inverse decay, that is, the term −dz2.(5)
The prey’s contribution to the predator growth rate is exw, fyw, and gzw, which is proportional to the available prey as well as the size of the predator population.Using the above assumptions, the following model is proposed:(1)dxdt=ax1−x−xw+xyzw,dydt=by1−y−yw+xyzw,dzdt=cz1−z−zw+xyzw,dwdt=−dw2+exw+fyw+gzw,where the coefficients a, b, c, d, e, f, and g are positive constants and the initial values x0,y0,z0,andw0 are all greater than zero. It is clear that the three teams of preys help each other, e.g., in foraging and in early warning against predation. Note that this help occurs only in the presence of predator. This is presented by the term xyzw in the prey equations.


## 3. The Analysis of the Model

First, let us find the equilibrium point of system (1). In order to find the equilibrium points of equation (1), equation (1) can be set to zero as follows:(2)ax1−x−xw+xyzw=0,by1−y−yw+xyzw=0,cz1−z−zw+xyzw=0,−dw2+exw+fyw+gzw=0.Solving equation (2), we can obtain the equilibrium point E (x, y, z, w) as follows: E0 (0, 0, 0, 0), E1 (1, 0, 0, 0), E2 (0, 1, 0, 0), E3 (0, 0, 1, 0), E4 (0, 1, 1, 0), E5 (1, 0, 1, 0), E6 (1, 1, 0, 0), E7 (1, 1, 1, 0), E8(0,0,(dc/dc+g),gc/dc+g), E9(0,(db/db+f),0,fb/db+f), E10(da/da+e,0,0,ea/da+e), E110,bcd−cg+bg/bcd+cf+bg,bcd+cf−bf/bcd+cf+bg,bcf+bcg/bcd+cf+bg, E12(acd−cg+ag/acd+ce+ag,0,acd+ce−ae/acd+ce+ag,ace+acg/acd+ce+ag), E13(abd−bf+af/abd+be+af,abd+be−ae/abd+be+af,0,abe+abf/abd+be+af), and E141,1,1,e+f+g/d.Proposition 1.
Local stability analysis shows that system (1) has eight unstable equilibrium solutions E0, E1, E2, E3, E4, E5, E6, and E7.Proof.
The Jacobian matrix of system (1) is given by the following:(3)J=a1−2x−w1−yzxzwxyw−x1−yzyzwb1−2y−w1−xzxyw−y1−xzyzwxzwc1−2z−w1−xy−z1−xyewfwgw−2dw+ex+fy+gz.
The characteristic equation is as follows:(4)J−λE=a1−2x−w1−yz−λxzwxyw−x1−yzyzwb1−2y−w1−xz−λxyw−y1−xzyzwxzwc1−2z−w1−xy−λ−z1−xyewfwgw−2dw+ex+fy+gz−λ=0.
By substituting the pointE0 (0, 0, 0, 0) in equation (4), we get λ = 0, a, b, and c, which has three positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE1 (1, 0, 0, 0) in equation (4), we get λ = −a, b, c, and e, which has three positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE2 (0, 1, 0, 0) in equation (4), we get λ = a, −b, c, and e, which has three positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE3 (0, 0, 1, 0) in equation (4), we get λ = a, b, −c, and e, which has three positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE4 (0, 1, 1, 0) in equation (4), we get λ = a, −b, −c, and f + g, which has two positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE5 (1, 0, 1, 0) in equation (4), we get λ = −a, b, −c, and e + g, which has two positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE6 (1, 1, 0, 0) in equation (4), we get λ = −a, −b, c, and e + f, which has two positive eigenvalues. So, it is an unstable equilibrium point.
Similarly, by substituting the pointE7 (1, 1, 1, 0) in equation (4), we get λ = −a, −b, −c, and e + f + g, which has one positive eigenvalue. So, it is an unstable equilibrium point.
So, all of them are unstable equilibrium points.
By substituting the pointE8(0,0,(dc/dc+g),gc/dc+g) in equation (4), we get λ=a−gc/dc+g,b−gc/dc+g,−c,and−dgc/dc+g.
The equilibrium solutionE8(0,0,(dc/dc+g),gc/dc+g) is locally asymptotically stable under the condition: {a<gc/dc+gb<gc/dc+g.
Numerical simulations agree with these results. Leta = 0.5, b = 0.5, c = 2, d = 1, e = 1, f = 2, and g = 2. We get the stable solution (0, 0, 0.5, 1) as shown in Figure 1.
Similarly, from equilibrium pointE9(0,(db/db+f),0,fb/db+f), we can obtain the eigenvalue as follows: λ=a−fb/db+f,c−fb/db+f,−b,and−dfb/db+f.
When the following conditions are met and the real part of all eigenvalues is less than 0, it is locally asymptotically stable:{a<fb/db+fc<fb/db+f.
The above conclusion is proved by numerical simulation, witha = 0.5, b = 2, c = 0.5, d = 1, e = 1, f = 2, and g = 2. The equilibrium point (0, 0.5, 0, 1) is obtained, and the simulation results are shown in Figure 2.
Similarly, from equilibrium pointE10da/da+e,0,0,(ea/da+e), we can obtain the eigenvalue as follows: λ=−a,b−ea/da+e,c−ea/da+e,−dea/da+e.
When the following conditions are met, it is locally asymptotically stable:{b<ea/da+ec<ea/da+e.
The above conclusion is proved by numerical simulation, witha = 2, b = 0.5, c = 0.5, d = 1, e = 2, f = 2, and g = 2. The equilibrium point (0.5, 0, 0, 1) is obtained, and the simulation results are shown in Figure 3.
By substituting the pointE110,bcd−cg+bg/bcd+cf+bg,bcd+cf−bf/bcd+cf+bg,bcf+bcg/bcd+cf+bg in matrix (3), we can obtain the Jacobian matrix:(5)J=a−w1−yz000yzwb1−2y−w0−yyzw0c1−2z−w−zewfwgw−2dw+fy+gz=a−w1−yzb1−2y−w0−y0c1−2z−w−zfwgw−dw.
In addition, the Jacobian characteristic determinant is as follows:(6)J−λE=a−w1−yz−λb1−2y−w−λ0−y0c1−2z−w−λ−zfwgw−dw−λ=−A−λλ3+dw−B−Cλ2+CB+gwz+fwy−dwB−dwCλ+dwBC−gwzB+fwyC,where A=a−w1−yz,B=b1−2y−w,andC=c1−2z−w.
Using the Routh–Hurwitz condition, the following condition is necessary and sufficient that all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part:(7)A<0,dw−B−C>0,CB+gwz+fwy−dwB−dwC>0,dwBC−gwzB+fwyC>0,dw−B−CCB+gwz+fwy−dwB−dwC>dwBC−gwzB+fwyC.
The above conclusion is proved by numerical simulation, witha = 0.5, b = 2, c = 3, d = 1, e = 2, f = 2, and g = 2. The equilibrium point (0, 0.25, 0.5, 1.5) is obtained, and the simulation results are shown in Figure 4.
Similarly, by substituting the pointE12((acd−cg+ag/acd+ce+ag),0,acd+ce−ae/acd+ce+ag,ace+acg/acd+ce+ag) in matrix (3), we can obtain the Jacobian matrix:(8)J=a1−2x−wxzw0−x0b−w1−xz000xzwc1−2z−w−zewfwgw−dw=b−w1−xza1−2x−w0−x0c1−2z−w−zewgw−dw.
In addition, the Jacobian characteristic determinant is as follows:(9)J−λE=b−w1−xz−λa1−2x−w−λ0−x0c1−2z−w−λ−zewgw−dw−λ=−B1−λλ3+dw−A1−C1λ2+A1C1+zgw+ewx−dwA1−dwC1λ+dwA1C1−zgwA1+ewxC,where {A1=a1−2x−wB1=b−w1−xzC1=c1−2z−w.
Similarly, using the Routh–Hurwitz condition, the following condition is necessary and sufficient that all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part:(10)B1<0,dw−A1−C1>0,A1C1+zgw+ewx−dwA1−dwC1>0,dwA1C1−zgwA1+ewxC1>0,dw−A1−C1A1C1+zgw+ewx−dwA1−dwC1>dwA1C1−zgwA1+ewxC1.
The above conclusion is proved by numerical simulation, witha = 2, b = 0.5, c = 3, d = 1, e = 2, f = 2, and g = 2. The equilibrium point (0.25, 0, 0.5, 1.5) is obtained, and the simulation results are shown in Figure 5.
By substituting the pointE13((abd−bf+af/abd+be+af),abd+be−ae/abd+be+af,0,abe+abf/abd+be+af) in matrix (3), we can obtain the Jacobian matrix:(11)J=a1−2x−w0xyw−x0b1−2y−wxyw−y00c−w1−xy0ewfwgw−dw=c−w1−xya1−2x−w0−x0b1−2y−w−yewfw−dw.
In addition, the Jacobian characteristic determinant is as follows:(12)J−λE=c−w1−xy−λa1−2x−w−λ0−x0b1−2y−w−λ−yewfw−dw−λ=−C2−λλ3+dw−A2−B2λ2+A2B2+ewx+fwy−dwB2−dwA2λ+dwA2B2−ewxB2+fwyA2=0,where A2=a1−2x−w,B2=b1−2y−w,andC2=c−w1−xy
Using the Routh–Hurwitz condition, the following condition is necessary and sufficient that all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part:(13)C2<0,dw−A2−B2>0,A2B2+ewx+fwy−dwB2−dwA2>0,dwA2B2−ewxB2+fwyA2>0,dw−B2−A2A2B2+ewx+fwy−dwB2−dwA2>dwA2B2−ewxB2+fwyA2.
The above conclusion was proved by numerical simulation, witha = 2, b = 3, c = 0.5, d = 1, e = 2, f = 2, and g = 2. The equilibrium point (0.25, 0.5, 0, 1.5) is obtained, and the simulation results are shown in Figure 6.
By substituting the pointE141,1,1,e+f+g/d in matrix (3), we can obtain the Jacobian matrix and the Jacobian characteristic determinant as follows:(14)J−λE=−a−λww0w−b−λw0ww−c−λ0ewfwgw−e−f−g−λ=−e−f−g−λ−a−λwww−b−λwww−c−λ=−−e−f−g−λλ3+a+b+cλ2+ab+ac+bc−3w2λ+abc−a+b+cw2−2w3=0.
It is known from above equation thatλ1 = -e-f-g < 0. Using the Routh–Hurwitz condition, the following condition is necessary and sufficient that all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part:(15)ab+ac+bc−3w2>0,abc−a+b+cw2−2w3>0,abc−a+b+cw2−2w3<a+b+cab+ac+bc−3w2.
When above conditions are met,E141,1,1,e+f+g/d is locally asymptotically stable.
The above conclusion is proved by numerical simulation, witha = 4, b = 4, c = 4, d = 3, e = 2, f = 2, and g = 2. The equilibrium point (1, 1, 1, 2) is obtained, and the simulation results are shown in Figure 7.
Now, we consider the above system but without team interaction help terms:(16)dxdt=ax1−x−xw,dydt=by1−y−yw,dzdt=cz1−z−zw,dwdt=−dw2+exw+fyw+gzw.
It is easy to get its equilibrium point as follows:x=bag−bgc−cfb+afc+dbca/dbca+bag+bec+afc, y=dbca+bag+bec−agc−aec/dbca+bag+bec+afc, z=dbca−afb−bea+afc+bec/dbca+bag+bec+afc, and w=cabf+e+g/dbca+bag+bec+afc.
To prove the global stability of the above equilibrium point (the internal solution) of system (16), we use the following propositions.Figure 1
The valuesa = 0.5, b = 0.5, c = 2, d = 1, e = 1, f = 2, and g = 2 are used. So, we get the stable solution (x = 0, y = 0, z = 0.5, and w = 1).Figure 2
The valuesa = 0.5, b = 2, c = 0.5, d = 1, e = 1, f = 2, and g = 2 are used. So, we get the stable solution (x = 0, y = 0.5, z = 0, and w = 1).Figure 3
The valuesa = 2, b = 0.5, c = 0.5, d = 1, e = 2, f = 2, and g = 2 are used. So, we get the stable solution (x = 0.5, y = 0, z = 0, and w = 1).Figure 4
The valuesa = 0.5, b = 2, c = 3; d = 1, e = 2, f = 2, and g = 2 are used. So, we get the stable solution (x = 0, y = 0.25, z = 0.5, and w = 1.5).Figure 5
The valuesa = 2, b = 0.5, c = 3, d = 1, e = 2, f = 2, and g = 2 are used. So, we get the stable solution (x = 0.25, y = 0, z = 0.5, and w = 1.5).Figure 6
The valuesa = 2, b = 3, c = 0.5, d = 1, e = 2, f = 2, and g = 2 are used. So, we get the stable solution (x = 0.25, y = 0.5, z = 0, and w = 1.5).Figure 7
The valuesa = 4, b = 4, c = 4, d = 3, e = 2, f = 2, and g = 2 are used. So, we get the stable solution (x = 1, y = 1, z = 1, and w = 2).Proposition 2.
(see [28]). The system du/dt=ru1−u has two equilibrium solutions, u = 0 and u = 1, where the constant r > 0. The nonzero equilibrium solution u = 1 is globally stable for all solutions u with the initial condition u0>0.Proposition 3.
Whena, b, c, d, e, f, and g are all positive, the internal solutions of (16) are globally asymptotically stable for all solutions with the initial conditions x0>0, y0>0, z0>0, and w0>0.Proof.
From the first equation in system (16), it is known that dx/dt≤ax1−x. According to Proposition 2, we can obtain that x≤U1x,whereU1x=1.
From the second equation in system (16), it is known that dy/dt≤by1−y. According to Proposition 2, we can obtain that y≤U1y,whereU1y=1.
From the third equation in system (16), it is known that dz/dt≤cz1−z. According to Proposition 2, we can obtain that z≤U1z,whereU1z=1. Therefore, from the fourth equation in system (16), it is known that dw/dt≤−dw2+eU1x+fU1y+gU1zw anddw/dt≤dweU1x+fU1y+gU1z/d−w.
According toProposition 2, we can obtain that w≤U1w,whereU1w=eU1x+fU1y+gU1z/d.
On the other hand,(17)dxdt≥ax1−x−xU1w,dxdt≥ax1−U1wa−x.
UsingProposition 2, we can obtain that x≥L1x,whereL1x=1−U1w/a.
By the same way,y≥L1y,whereL1y=1−U1w/b, and z≥L1z,whereL1z=1−U1w/c.
Similarly,dw/dt≥−dw2+eL1x+fL1y+gL1zw and dw/dt≥dweL1x+fL1y+gL1z/d−w.
Using Proposition2, we can obtain that w≥L1w,whereL1w=eL1x+fL1y+gL1z/d.
Continuing the above process, one gets a sequence of upper limitsUix,Uiy,Uiz,andUiw and lower limits Lix,Liy,Liz,andLiw, where i = 2, 3, related by the relations:(18)Uix=1−Liwa,Uiy=1−Liwb,Uiz=1−Liwc,Uiw=eUix+fUiy+gUizd,Lix=1−Uiwa,Liy=1−Uiwb,Liz=1−Uiwc,Liw=eLix+fLiy+gLizd,i=2,3,...,n,whose solutions are Uix=Lix=bag−bgc−cfb+afc+dbca/dbca+bag+bec+afcUiy=Liy=dbca+bag+bec−agc−aec/dbca+bag+bec+afcUiz=Liy=dbca−afb−bea+afc+bec/dbca+bag+bec+afcUiw=Liw=cabf+e+g/dbca+bag+bec+afc.
This completes the proof.
The conclusion is proved by numerical simulation. When arbitrarily taking positive number asa = 4, b = 6, c = 8, d = 2, e = 3, f = 2, and g = 2 and x0=0.2, y0=0.1, z0=0.3, and w0=1.4, we can obtain internal solution in equation (16) as (0.475, 0.65, 0.7375, 2.1), and the simulation results are shown in Figure 8, which shows that the simulation result agrees with theory result of internal solution (equilibrium point). When arbitrarily taking positive number as a = 2, b = 4, c = 5, d = 5, e = 2, f = 2, and g = 2 and x0=1.5, y0=0.1, z0=1.3, and w0=1.4, we can obtain internal solution in equation (16) as (0.5652, 0.7826, 0.8261, 0.8696), and the simulation results are shown in Figure 9, which shows that the simulation result agrees with theory result of internal solution (equilibrium point).Figure 8
Positive parametersa = 4, b = 6, c = 8, d = 2, e = 3, and f = 2 corresponding to internal solution in equation (16) as (0.475, 0.65, 0.7375, 2.1).Figure 9
Positive parametersa = 2, b = 4, c = 5, d = 5, e = 2, f = 2, and g = 2 corresponding to internal solution in equation (16) as (0.5652, 0.7826, 0.8261, 0.8696).Proposition 4.
System (16) is persistent.Proof .
According to Proposition3, we have proved that for all t > 0, all solutions with positive initial conditions satisfy {0<Lix≤x≤UixLiy≤y≤UiyLiy≤z≤UizLiw≤w≤Uiw.
Thus, system (16) is persistent.


## 4. Comparison with Reference [28]

Reference [28] studied the two-prey one-predator model, while this paper studies the three-prey single-predator model. Compared with reference [28], our model in this paper has the following characteristics: (i) the number of differential equations for the system increases by one; (ii) the number of equilibrium points increases by half, from 8 in reference [28] to 15 in this paper; and (iii) the equilibrium point converges faster, and the convergence time of numerical simulation is about 10 s in reference [28] and 2 s in this paper.


## 5. Conclusions

In this paper, a three-prey single-predator continuous time nonlinear system is studied and a model of fourth-order nonlinear differential equation is proposed. The equilibrium point, the characteristic root, and the stability of the system solution of the differential equation are analyzed. When considering preys help each other, there are 15 equilibria in the system in which 8 equilibria are unstable and other 7 equilibria are locally asymptotically stable. When not considering preys help each other, all equilibria are globally asymptotically stable.



---

*Source: 8869989-2020-10-10.xml*