# Analysis for a Delayed Three-Species Predator-Prey Model with Feedback Controls and Prey Diffusion

Authors: Lili Jia
Journal: Journal of Mathematics
 (2020)

---

## Abstract

In this paper, we study a class of 3-species multidelay Lotka–Volterra ratio-dependent predator-prey model with feedback controls and prey diffusion. By using the theory of delay differential inequalities and developing some new analysis methods as well as constructing a suitable Lyapunov function, some sufficient conditions are obtained to guarantee the permanence of the system and the global attractivity of the positive solution for the predator-prey system. Furthermore, the corresponding periodic system is discussed, and some conditions are established about the existence, uniqueness, and stability of the positive periodic solution for the periodic system by using the fixed-point theory and theoretical analysis. In additional, some numerical solutions of the equations describing the system are given to verify the obtained criteria. Finally, the corresponding stochastic time-delay predator-prey model with multiplicative noise sources is solved numerically, and some interesting new dynamics behavior is obtained.


---

## Body

## 1. Introduction

The dynamic relationship between predator and prey is one of the core problems in population ecology. Especially, the study of persistence and globally attractive as well as periodic of biological species is a very important research direction in population dynamics behavior research. The traditional Lotka–Volterra predator model has been studied by more and more scholars, and a lot of research results have been achieved (see [1–4] and the references cited therein). In some practical applications, such as population growth, reproductive physiology, and organic insensitivity, the growth rate of the species does not respond immediately. In other words, we have to face a time delay. The existence of delay in the system greatly affects the stability of the system. In some cases, the system is stable when the delay value is smaller. However, many scholars found that longer delays may destroy the stability of the system. Some related results are shown in [5–9] and the references cited therein. In [9], the authors considered a three-species Lotka–Volterra-type competitive-mutualism system with discrete time delays and established some conditions on the boundedness, permanence, and global attractivity for the system. Furthermore, in the real world, species can diffuse between patches. So, we should consider the effect of dispersal on the permanence and global stability of the ecosystem. Some research results can be seen, for example, in [10–13]. Song and Chen [10] studied a two-species predator-prey system with diffusion. It is proved in this paper that the system can be made persistent; furthermore, if the system is a periodic system, it can have a strictly positive periodic orbit which is globally asymptotically stable under the appropriate conditions. Moreover, predators have to search for food and have to share or compete for food, and a more suitable general Lotka–Volterra predator-prey model should be based on the so-called “ratio-dependent” theory. That is to say, the per capita predator growth rate should be a function of the ratio of prey to predator abundance (e.g., see [14–17]). In [14], Wang and Wang considered a ratio-dependent predator-prey system with two competing preys predated by one predator and showed that this system is permanent and globally asymptotically stable under some appropriate conditions. In order to study the combined effects of time delay, dispersion, and ratio-dependent on the dynamic behavior of predator-prey systems, some models have been studied by many authors [18–24]. Xu and Chaplain [18] studied a three-species predator-prey model both with time delay and ratio-dependent. They proved that the system is uniformly persistent under some appropriate conditions, and by means of constructing suitable Lyapunov functional, sufficient conditions are derived for the global asymptotic stability of the positive equilibrium of the system. In 2004, Xu et al. [19] studied a Lotka–Volterra predator-prey model with dispersion and time delays. By using Ganies and Mawhin’s continuation theorem of coincidence degree theory and constructing the appropriate Lyapunov function, a set of easily verifiable sufficient conditions are derived for the existence, uniqueness, and global stability of positive periodic solutions of the system. In [21], Sun and Yuan considered a nonautomous mixture model both with ratio-dependent and diffusion. It is shown that the system is uniformly and persistently related to the dispersion rate. Furthermore, the sufficient conditions are obtained for the global asymptotic stability of a periodic solution for the system. In [23], the authors considered the following delayed two-predator-one-prey ratio-dependent model with patch dispersion:(1)x˙1=x1ta1−a11x1t−a13x3tm13x3t+x1t−a14x4tm14x4t+x1t+D1x2t−x1t,x˙2=x2ta2−a22x2t+D2x1t−x2t,x˙3=x3t−a3+a31x1t−τ1m13x3t−τ1+x1t−τ1,x˙4=x4t−a4+a41x1t−τ2m14x4t−τ2+x1t−τ2,where xiti=1,2 represents the prey density in the i-th patch, xjtj=3,4 represents the predator density, τ1,τ2>0 are constant delays due to gestation, Dii=1,2 are positive constants that denote the dispersal rates, ai, aiji,j=1,2,3,4, and m13 and m14 are positive constants that denote interaction among species. It is shown that system (1) is uniformly persistent under some appropriate conditions, and sufficient conditions are obtained for the global stability of the positive equilibrium of the system.On the contrary, in the real world, ecosystems are continuously distributed by unpredictable forces which could cause a stable system to become unstable or cause the species to fluctuate. In order to make an unstable system or equilibrium to become stable one, we can introduce an indirect feedback control variables. The feedback control mechanisms can be implemented through harvesting or culling procedures or some biological control schemes. So, it is necessary to study models with control variables. Recently, there has been a lot of literature related to the study of Lotka–Volterra-type system with feedback controls [25–30]. In 2003, Gopalsamy and Weng [25] studied a two-species competition system with feedback controls. Some sufficient conditions are obtained for the existence of a globally attracting positive equilibrium of the system. To study the influence of the patch diffusion, ratio-dependent, and feedback control on existence of periodic solution, Xie and Weng [26] considered a predator-prey model with patch-diffusion, ratio-dependent, and feedback control. By changing the ordinary differential system into an operator equation and using Mawhin’s continuation theorem, the existence of at least one positive periodic solution for this model is proved. In order to show that whether the feedback control variables play an essential role on the persistent property and global stability of Lotka–Volterra cooperative systems or not, Liang [28] discussed a two-species system with time delays and feedback control and proved that the only positive equilibrium point of the system is globally stable. In 2016, Wang et al. [29] studied a ratio-dependent Lotka–Volterra predator-prey model with feedback control. Some sufficient conditions which guarantee the globally attractive of positive solution for the predator-prey model are obtained. In [30], Chen considered a Lotka–Volterra biological system with feedback controls and obtained some conditions for the permanence and global attractivity of this system.However, to the best of the authors’ knowledge, there are few scholars who study the three-species multidelay Lotka–Volterra ratio-dependent predator-prey model with feedback controls and prey diffusion. So, based on system (1) and motivated by the above works, in this paper, we propose and investigate the following three-species multidelay ratio-dependent predator-prey model with feedback controls and prey diffusion:(2)x˙1t=x1tr1t−a11tx1t−a13tx3tb13t+x1t−a14tx4tb14t+x1t−d1tu1t+D1tx2t−x1t,x˙2t=x2tr2t−a22tx2t+D2tx1t−x2t,x˙3t=x3t−r3t−a33tx3t−τ1+a31tx1t−τ2b13t+x1t−τ2+d2tu2t,x˙4t=x4t−r4t−a44tx4t−τ3+a41tx1t−τ4b14t+x1t−τ4+d3tu3t,u˙1t=e1t−f1tu1t+q1tx1t,u˙2t=e2t−f2tu2t−q2tx3t,u˙3t=e3t−f3tu3t−q3tx4t,where xiti=1,2 denote the prey density in the i-th patch and xjtj=3,4 represents the predator density. The predator species is confined to patch 1 while the prey species can disperse between two patches, riti=1,2 denote the intrinsic growth rate of the prey species, and rjtj=3,4 are the death rate of the predators. aiiti=1,2,3,4 denote the internal competitive coefficient of species, a13t and a14t are the ratio of prey by predator, a31t and a41t represent the nutrient absorption ratio of the predator after predation, and Dit i=1,2 are the dispersion rate of prey species. uiti=1,2,3 are the feedback control terms, eit, fit, and qiti=1,2,3 are the coefficient functions of the control equations, and τ1, τ2, τ3, and τ4 are positive constants representing delays. All of the coefficients in the model are continuous and positive-bounded functions defined on 0,+∞.Remark 1.
It is easy to see that the designed controller of this paper is similar to that in [29], and the designed controller in [29] is as follows:(3)u˙1t=e1t−f1tu1t+q1tx1t,u˙2t=e2t−f2tu2t+q2tx2t,u˙3t=e3t−f3tu3t−q3tx3t.
Namely,(4)u1t=e−∫0tf1tdt∫0te1t+q1tx1te∫0tf1tdtdt+φ10,u2t=e−∫0tf2tdt∫0te2t+q2tx2te∫0tf2tdtdt+φ20,u3t=e−∫0tf3tdt∫0te3t−q3tx3te∫0tf3tdtdt+φ30.
And the designed controller in this paper is as follows:(5)u˙1t=e1t−f1tu1t+q1tx1t,u˙2t=e2t−f2tu2t−q2tx3t,u˙3t=e3t−f3tu3t−q3tx4t.
That is,(6)u1t=e−∫0tf1tdt∫0te1t+q1tx1te∫0tf1tdtdt+φ10,u2t=e−∫0tf2tdt∫0te2t−q2tx3te∫0tf2tdtdt+φ20,u3t=e−∫0tf3tdt∫0te3t−q3tx4te∫0tf3tdtdt+φ30.
However, the state equations are different. The model in [29] contains two preys and one predator, while the model in this paper contains one prey and two predators. In particular, the new model contains prey diffusion and delays.
Due to biological interpretation of system (2), it is reasonable to consider only the positive solution of (2), in other words, to take admissible initial conditions(7)xit=ϕit≥0,t∈−τ,0,ϕi0>0,i=1,2,3,4,uj0=φj0>0,j=1,2,3,where τ=maxτ1,τ2,τ3,τ4.
Obviously, the solutions of system (2) with initial data (7) are positive for all t≥0.
This paper is organized as follows: in Section2, we give the conditions for the permanence to system (2) by means of the theory of delay differential inequalities. In Section 3, by constructing an appropriate Lyapunov function, we derive sufficient conditions for the global attractivity of positive solution for system (2). In Section 4, by using some new analysis method and the Brouwer fixed-point theory, we obtain some conditions for the existence, uniqueness, and stability of a positive periodic solution for the corresponding periodic system. Some numerical solutions of the equations describing the system are shown in Section 5 to verify the validity of the obtained criteria.


## 2. Permanence

Firstly, we introduce the following notations and definitions. For a functiongt defined on t0,+∞, we set(8)gm=supgtt0<t<∞,gl=infgtt0<t<∞.Definition 1.
System (2) is called permanent if there exist positive constants Mi, Nj, mi, nj (i=1,2,3,4, j=1,2,3), and T such that mi≤xit≤Mi and nj≤ujt≤Nj for any positive solution Zt=x1t, x2t,x3t,x4t,u1t,u2t,u3t of system (29) as t>T.
As a direct corollary of Lemma 2.1 of Chen [31], we have the following.Lemma 1.
Ifa>0,b>0, and x˙≥b−ax, when t≥0 and x0>0, we have liminft⟶+∞xt≥b/a. If a>0, b>0, and x˙≤b−ax, when t≥0 and x0>0, we have limsupt⟶+∞xt≤b/a.As a direct corollary of Lemma 2.2 of Chen [31], we have the following.Lemma 2.
Ifa>0,b>0, and x˙≥b−ax, when t≥0 and x0>0, we have liminft⟶+∞xt≥b/a. If a>0, b>0, and x˙≤xb−ax, when t≥0 and x0>0, we have limsupt⟶+∞xt≤b/a.Lemma 3 (see [32]).
Assume that, foryt>0, it holds that y˙t≤ytλ−∑k=0mμkyt−kτ+D, with initial conditions yt=ϕt≥0 for t∈−mτ,0 and ϕ0>0, where λ>0, μk≥0k=0,1,…,m, μ=∑k=0mμk>0, and D≥0 are constants. Then, there exists a positive constant My<+∞ such that(9)limsupt⟶+∞yt≤My=−Dλ+Dλ+y∗expλmτ<+∞,where y=y∗ is the unique solution of equation yλ−μy+D=0.Lemma 4 (see [32]).
Assume that, foryt>0, it holds that y˙t≥ytλ−∑k=0mμkyt−kτ. If (9) holds, then, there exists a positive constant my>0 such that, for μ=∑k=0mμk>0, we have(10)liminft⟶+∞yt≥my=λμexpλ−μMymτ>0.
For system (2), we let(11)M1=M2=maxr1ma11l,r2ma22l,N1=e1m+q1mM1f1l,N2=e2mf2l,N3=e3mf3l,M3=a31mM1/b13l+M1+d2mN2−r3la33l×expa31mM1b13l+M1+d2mN2−r3lτ1,M4=a41mM1/b14l+M1+d3mN3−r4la44l×expa41mM1b14l+M1+d3mN3−r4lτ3,m1=r1l−a13mM3/b13l−a14mM4/b14l−d1mN1−D1ma11m,m2=r2l−D2ma22m,n1=e1l+q1lm1f1m,n2=e2l−q2mM3f2m,n3=e3l−q3mM4f3m,m3=a31lm1/b13m+m1+d2ln2−r3ma33m×expa31lm1b13m+m1+d2ln2−r3m−a33mM3τ1,m4=a41lm1/b14m+m1+d3ln3−r4ma44m×expa41lm1b14m+m1+d3ln3−r4m−a44mM4τ3.Next, we state and prove our main result with respect to the permanent for system (2).Theorem 1.
Suppose that system (2) satisfies the initial conditions (7) and following assumptions:(H1): r3l<a31mM1/b13l+M1+d2mN2(H2): r4l<a41mM1/b14l+M1+d3mN3(H3): r1l>a13mM3/b13l+a14mM4/b14l+d1mN1+D1m(H4): r2l>D2m(H5): e2l>q2mM3(H6): e3l>q3mM4(H7): r3m<a31lm1/b13m+m1+d2ln2(H8): r4m<a41lm1/b14m+m1+d3ln3
Then, system (2) is permanent.Proof.
According to the first and the second equations of system (2), we define W1t=maxx1t,x2t and calculate the upper right derivative of W1t along the positive solution of system (2), and we have that(P1): ifx1t≥x2t, then(12)D+W1t=x˙1t=x1tr1t−a11tx1t−a13tx3tb13t+x1t−a14tx4tb14t+x1t−d1tu1t+D1tx2t−x1t≤x1tr1t−a11tx1t≤x1tr1m−a11lx1t=W1tr1m−a11lW1t.(P2): ifx1t≤x2t, then(13)D+W1t=x˙2t=x2tr2t−a22tx2t+D2tx1t−x2t≤x2tr2t−a22tx2t≤x2tr2m−a22lx2t=W1tr2m−a22lW1t.
It follows from (P1) and (P2) that(14)D+W1t≤W1trim−aiilW1t,i=1,2.
By (14), we can also deduce the following:(A)
IfW10=maxx10,x20≤M1, then we have maxx1t,x2t≤M1,t≥0(B)
IfW10=maxx10,x20>M1, take appropriate α>0, we have the following three possibilities:(a)
W10=x10>M1,x10>x20(b)
W10=x20>M1,x10<x20(c)
W10=x10=x20>M1
If (a) holds, then there existsε>0,t∈0,ε such that W1t=x1t>M1. Then, we get(15)D+W1t=x˙1t≤a11lW1tr1ma11l−W1t≤−α<0.
Similarly, if (b) holds, there existsε>0,t∈0,ε, and W1t=x2t>M1 to derive(16)D+W1t=x˙2t≤a22lW1tr2ma22l−W1t≤−α<0.
Moreover, in case (c), we also have(17)D+W1t=x˙it≤aiilW1trimaiil−W1t≤−α<0,i=1or2.
From what we have discussed above, it yields that ifW10>M1, then W1t is strictly monotone decreasing with speed at least α, so there exists T1>0 such that, for any time t≥T1, W1t=maxx1t,x2t≤M1, which admits(18)limsupt⟶+∞x1t≤M1=M2=maxr1ma11l,r2ma22l,(19)limsupt⟶+∞x2t≤M1=M2=maxr1ma11l,r2ma22l.
According to the fifth equation of system (2), we have(20)u˙1t=e1t−f1tu1t+q1tx1t≤e1m−f1lu1t+q1mM1.
Using Lemma1, we can deduce(21)limsupt⟶+∞u1t≤e1m+q1mM1f1l=N1.
Similarly, in light of the last two equations of system (26), we then have(22)u˙2t=e2t−f2tu2t−q2tx3t≤e2t−f2tu2t≤e2m−f2lu2t,u˙3t=e3t−f3tu3t−q3tx4t≤e3t−f3tu3t≤e3m−f3lu3t.
By utilizing Lemma1 again, we can get(23)limsupt⟶+∞u2t≤e2mf2l=N2,(24)limsupt⟶+∞u3t≤e3mf3l=N3.
From the third and fourth equations of system (2), we have(25)x˙3t=x3t−r3t−a33tx3t−τ1+a31tx1t−τ2b13t+x1t−τ2+d2tu2t≤x3t−r3l−a33lx3t−τ1+a31mM1b13l+M1+d2mN2,x˙4t=x4t−r4t−a44tx4t−τ3+a41tx1t−τ4b14t+x1t−τ4+d3tu3t≤x4t−r4l−a44lx4t−τ3+a41mM1b14l+M1+d3mN3.
By Lemma3, combining with H1andH2, we have(26)limsupt⟶+∞x3t≤a31mM1/b13l+M1+d2mN2−r3la33l×expa31mM1b13l+M1d2mN2−r3lτ1=M3,(27)limsupt⟶+∞x4t≤a41mM1/b14l+M1+d3mN3−r4la44l×expa41mM1b14l+M1+d3mN3−r4lτ3=M4.
On the contrary, from the first equation of system (2), it follows that(28)x˙1t=x1tr1t−a11tx1t−a13tx3tb13t+x1t−a14tx4tb14t+x1t−d1tu1t+D1tx2t−x1t≥x1tr1t−a11tx1t−a13tx3tb13t−a14tx4tb14t−d1tu1t−D1t≥x1tr1l−a11mx1t−a13mM3b13l−a14mM4b14l−d1mN1−D1m.
By Lemma2 and assumption H3, we have(29)liminft⟶+∞x1t≥r1la11m−a13mM3b13la11m−a14mM4b14la11m−d1mN1a11m−D1ma11m=m1.
For the second equation of system (2), it holds that(30)x˙2t=x2tr2t−a22tx2t+D2tx1t−x2t≥x2tr2t−a22tx2t−D2t≥x2tr2l−a22mx2t−D2m.
From assumptionH4 and by Lemma 2, we can get(31)liminft⟶+∞x2t≥r2l−D2ma22m=m2.
For the fifth equation of system (2), we have(32)u˙1t=e1t−f1tu1t+q1tx1t≥e1l−f1mu1t+q1lm1.
By Lemma1, it holds that(33)liminft⟶+∞u1t≥e1l+q1lm1f1m=n1.
For the sixth and the seventh equations of system (2), we can obtain(34)u˙2t=e2t−f2tu2t−q2tx3t≥e2l−f2mu2t−q2mM3,u˙3t=e3t−f3tu3t−q3tx4t≥e3l−f3mu3t−q3mM4.
By Lemma1, H5, and H6, we have(35)liminft⟶+∞u2t≥e2l−q2mM3f2m=n2,(36)liminft⟶+∞u3t≥e3l−q3mM4f3m=n3.
According to the third and the fourth equations of system (2), we have(37)x˙3t=x3t−r3t−a33tx3t−τ1+a31tx1t−τ2b13t+x1t−τ2+d2tu2t≥x3t−r3m−a33mx3t−τ1+a31lm1b13m+m1+d2ln2,x˙4t=x4t−r4t−a44tx4t−τ3+a41tx1t−τ4b14t+x1t−τ4+d3tu3t≥x4t−r4m−a44mx4t−τ3+a41lm1b14m+m1+d3ln3.
Using Lemma4, H7, andH8, we also derive that(38)liminft⟶+∞x3t≥a31lm1/b13m+m1+d2ln2−r3ma33m×expa31lm1b13m+m1+d2ln2−r3m−a33mM3τ1=m3,(39)liminft⟶+∞x4t≥a41lm1/b14m+m1+d3ln3−r4ma44m×expa41lm1b14m+m1+d3ln3−r4m−a44mM4τ3=m4.
Summarizing (18)–(39), we complete the induction of Theorem 1.Remark 2.
It is obvious that the proof method for Lemma 2.1 in [29] is the differential inequality theory, while that of this article is the delay differential inequality theory. From the result of Theorem 1, the feedback control items and delays have influence on the permanence of system (2). In other words, the densities of the prey and predator populations are uniformly bounded when system (2) satisfies initial conditions (7) and assumptions H1−H8.


## 3. Global Attractivity

In this section, the global attractivity of system (2) will be studied. To get the sufficient conditions for global attractivity of system (2), we firstly give the following definition and Lemma.Definition 2.
System (2) is said to be globally attractive if there exists a positive solution(40)Xt=x1t,x2t,x3t,x4t,u1t,u2t,u3t,of system (2) such that(41)limt⟶+∞xit−yit=0,limt⟶+∞uit−vit=0,i=1,2,3,4,j=1,2,3,for any other positive solution Yt=y1t,y2t,y3t,y4t,v1t,v2t,v3t of system (2).Lemma 5 (see [33]).
If the functionft:R+⟶R is uniformly continuous and the limit limt⟶+∞∫0tfsds exists and is finite, then limt⟶+∞ft=0.
Next, we state and prove our main result with respect to the global attractivity for system (2).Theorem 2.
Assume that system (2) satisfies H1−H8 and the following conditions:
(H9) Ai>0,Bj>0i=1,2,3,4,j=1,2,3, where(42)B1=f1l−d1m,B2=f2l−d2m−M3d2ma33mτ1,B3=f3l−d3m−M4d3ma44mτ3,A1=a11l−a13mM3b13l+m12−a14mM4b14l+m12−D2mm2−a31mb13mb13l+m12−M3a31ma33mτ1b13l+m1−a41mb14mb14l+m12−M4a41ma44mτ3b14l+m1−q1m,A2=a22l−D1mm1,A3=a33l−r3m+a33mM3+a31mM1b13l+M1+d2mN2a33mτ1−M3a33m2τ1−a13mb13m+a13mM1b13l+m12−q2m,A4=a44l−r4m+a44mM4+a41mM1b14l+M1+d3mN3a44mτ3−M4a44m2τ3−a14mb14m+a14mM1b14l+m12−q3m.
Then, system (2) is globally attractive.Proof.
Suppose thatx1t,x2t,x3t,x4t,u1t,u2t,u3t and y1t,y2t,y3t,y4t,v1t,v2t,v3t are any two different positive solutions of system (2). Then, from Theorem 1, there exist positive constants Mi, mi, Nj, nji=1,2,3,4,j=1,2,3, and T such that(43)mi≤xit,yit≤Mi,i=1,2,3,4,nj≤ujt,vjt≤Nj,j=1,2,3.
We define(44)V1t=lnx1t−lny1t,V2t=lnx2t−lny2t.
Calculating the upper-right derivative ofV1t and V2t along the solution of system (2), we have(45)∑i=12D+ViT=∑i=12D+lnxit−lnyit≤sgnx1t−y1t−a11tx1t−y1t−a13tb13t+a13ty1tx3t−y3tb13t+x1tb13t+y1t−a13ty3tx1t−y1tb13t+x1tb13t+y1t−a13ty3tx1t−y1tb13t+x1tb13t+y1t−a14tb14t+a14ty1tx4t−y4tb14t+x1tb14t+y1t−a14ty4tx1t−y1tb14t+x1tb14t+y1t−d1tu1t−v1t+D1tx2tx1t−y2ty1t+sgnx2t−y2t−a22tx2t−y2t+D2tx1tx2t−y1ty2t≤−a11t−a13ty3tb13t+x1tb13t+y1t−a14ty4tb14t+x1tb14t+y1t×x1t−y1t−a22tx2t−y2t+a13tb13t+a13ty1tb13t+x1tb13t+y1tx3t−y3t+a14tb14t+a14ty1tb14t+x1tb14t+y1tx4t−y4t+d1tu1t−v1t+D˜1t+D˜2t,where(46)D˜1t=D1tx2tx1t−y2ty1t,x1t>y1t,D1ty2ty1t−x2tx1t,x1t<y1t,D˜2t=D2tx1tx2t−y1ty2t,x2t>y2t,D2ty1ty2t−x1tx2t,x2t<y2t.
Next, we discussD˜1t under the following three cases:(a)
Ifx1t>y1t, then(47)D˜1t≤D1tx1tx2t−y2t≤D1mm1x2t−y2t.(b)
Ifx1t<y1t, then(48)D˜1t≤D1ty1ty2t−x2t≤D1mm1x2t−y2t.(c)
Ifx1t=y1t, the same conclusion as (a) and (b) holds.
Combining the conclusions in (a)–(c), we have(49)D˜1t≤D1mm1x2t−y2t.
A similar argument shows that(50)D˜2t≤D2mtm`2x1t−y1t.
By means of (45)–(50), we have(51)∑i=12D+Vit≤−a11l−a13mM3b13l+m12−a14mM4b14l+m12−D2mm2×x1t−y1t−a22l−D1mm1x2t−y2t+a13mb13m+a13mM1b13l+m12x3t−y3t+a14mb14m+a14mM1b14l+m12x4t−y4t+d1mu1t−v1t.
Define(52)V31t=lnx3t−lny3t.
Calculating the upper right derivative ofV31t along the solution of system (2), it follows that(53)D+V31t=sgnx3t−y3t−a33tx3t−τ1−y3t−τ1+a31tb13tx1t−τ2−y1t−τ2b13t+x1t−τ2b13t+y1t−τ2+d2tu2t−v2t=sgnx3t−y3t−a33tx3t−y3t+d2tu2t−v2t+a31tb13tx1t−τ2−y1t−τ2b13t+x1t−τ2b13t+y1t−τ2+a33t∫t−τ1tx˙3θ−y˙3θdθ=sgnx3t−y3t−a33tx3t−y3t+d2tu2t−v2t+a31tb13tx1t−τ2−y1t−τ2b13t+x1t−τ2b13t+y1t−τ2+a33t∫t−τ1tx3θ−r3θ−a33θx3θ−τ1+a31θx1θ−τ2b13θ+x1θ−τ2+d2θu2θ−y3θ−r3θ−a33θy3θ−τ1+a31θy1θ−τ2b13θ+y1θ−τ2+d2θv2θdθ=sgnx3t−y3t−a33tx3t−y3t+d2tu2t−y2t+a31tb13tx1t−τ2−y1t−τ2b13t+x1t−τ2b13t+y1t−τ2+a33t∫t−τ1tx3θ−y3θ×−r3θ−a33θy3θ−τ1+a31θy1θ−τ2b13θ+y1θ−τ2+d2θv2θ+x3θ−a33θx3θ−τ1−y3θ−τ1+a31θb13θ+x1θ−τ2x1θ−τ2−y1θ−τ2+d2θu2θ−v2θdθ≤−a33tx3t−y3t+d2tu2t−v2t+a31tb13tb13t+x1t−τ2b13t+y1t−τ2×x1t−τ2−y1t−τ2+a33t∫t−τ1tr3θ+a33θy3θ−τ1+a31θy1θ−τ2b13θ+y1θ−τ2+d2θv2θx3θ−y3θ+x3θa33θx3θ−τ1−y3θ−τ1+a31θb13θ+x1θ−τ2x1θ−τ2−y1θ−τ2+d2θu2θ−v2θdθ.
Define(54)V32t=∫t−τ1t∫sta33s+τ1r3θ+a33θy3θ−τ1+a31θy1θ−τ2b13θ+y1θ−τ2+d2θv2θ×x3θ−y3θ+x3θa33θx3θ−τ1−y3θ−τ1+a31θb13θ+x1θ−τ2x1θ−τ2−y1θ−τ2+d2θu2θ−v2θdθds.
From (53) and (54), we have(55)∑i=12D+V3it≤−a33tx3t−y3t+d2tu2t−v2t+a31tb13tb13t+x1t−τ2b13t+y1t−τ2×x1t−τ2−y1t−τ2+∫t−τ1ta33s+τ1dsr3t+a33ty3t−τ1+a31ty1t−τ2b13t+y1t−τ2+d2tv2t×x3t−y3t+M3∫t−τ1ta33s+τ1ds×a33tx3t−τ1−y3t−τ1+M3∫t−τ1ta33s+τ1ds×a31tb13t+x1t−τ2x1t−τ2−y1t−τ2+M3∫t−τ1ta33s+τ1ds×d2tu2t−v2t.
We also define(56)V33t=∫t−τ2ta31w+τ2b13w+τ2b13w+τ2+x1wb13w+τ2+y1w×x1w−y1wdw+M3∫t−τ1t∫ww+τ1a33s+τ1a33w+τ1×x3w−y3wdsdw+M3∫t−τ2t∫w+τ2−τ1w+τ2a33s+τ1a31w+τ2b13w+τ2+x1w×x1w−y1wdsdw,and let(57)V3t=V31t+V32t+V33t.
According to (55) and (56), calculating the upper right derivative of V3t, we have(58)D+V3t≤−a33l−r3m+a33mM3+a31mM1b13l+M1+d2mN2∫t−τ1ta33ms+τ1ds−M3a33m∫tt+τ1a33ms+τ1dsx3t−y3t+a31mb13mb13l+m12+M3a31mb13l+m1×∫t+τ2−τ1t+τ2a33ms+τ1dsx1t−y1t+d2m+M3d2m∫t−τ1ta33ms+τ1ds×u2t−v2t≤−a33l−r3m+a33mM3+a31mM1b13l+M1+d2mN2a33mτ1−M3a33m2τ1x3t−y3t+a31mb13mb13l+m12+M3a31ma33mτ1b13l+m1x1t−y1t+d2m+M3d2ma33mτ1u2t−v2t.
Similarly, we define(59)V41t=lnx4t−lny4t.
Calculating the upper right derivative ofV41t along the solution of system (2), we get(60)D+V41t=sgnx4t−y4t−a44tx4t−τ3−y4t−τ3+a41tb14tx1t−τ4−y1t−τ4b14t+x1t−τ4b14t+y1t−τ4+d3tu3t−v3t=sgnx4t−y4t−a44tx4t−y4t+d3tu3t−v3t+a41tb14tx1t−τ4−y1t−τ4b14t+x1t−τ4b14t+y1t−τ4+a44t∫t−τ3tx˙4θ−y˙4θdθ=sgnx4t−y4t−a44tx4t−y4t+d3tu3t−v3t+a41tb14tx1t−τ4−y1t−τ4b14t+x1t−τ4b14t+y1t−τ4+a44t∫t−τ3tx4θ−r4θ−a44θx4θ−τ3+a41θx1θ−τ4b14θ+x1θ−τ4+d3θu3θ−y4θ−r4θ−a44θy4θ−τ3+a41θy1θ−τ4b14θ+y1θ−τ4+d3θv3θdθ=sgnx4t−y4t−a44tx4t−y4t+d3tu3t−v3t+a41tb14tx1t−τ4−y1t−τ4b14t+x1t−τ4b14t+y1t−τ4+a44t∫t−τ3tx4θ−y4θ×−r4θ−a44θy4θ−τ3+a41θy1θ−τ4b14θ+y1θ−τ4+d3θv3θ+x4θ−a44θx4θ−τ3−y4θ−τ3+a41θb14θ+x1θ−τ4x1θ−τ4−y1θ−τ4+d3θu3θ−v3θdθ≤−a44tx4t−y4t+d3tu3t−v3t+a41tb14tb14t+x1t−τ4b14t+y1t−τ4×x1t−τ4−y1t−τ4+a44t∫t−τ3tr4θ+a44θy4θ−τ3+a41θy1θ−τ4b14θ+y1θ−τ4+d3θv3θx4θ−y4θ+x4θa44θx4θ−τ3−y4θ−τ3+a41θb14θ+x1θ−τ4x1θ−τ4−y1θ−τ4+d3θu3θ−v3θdθ.
We define(61)V42t=∫t−τ3t∫sta44s+τ3r4θ+a44θy4θ−τ3+a41θy1θ−τ4b14θ+y1θ−τ4+d3θv3θx4θ−y4θ+x4θa44θx4θ−τ3−y4θ−τ3+a41θb14θ+x1θ−τ4x1θ−τ4−y1θ−τ4+d3θu3θ−v3θdθds.
By (60) and (61), we have(62)∑i=12D+V4it≤−a44tx4t−y4t+d3tu3t−v3t+a41tb14tb14t+x1t−τ4b14t+y1t−τ4×x1t−τ4−y1t−τ4+∫t−τ3ta44s+τ3ds×r4t+a44ty4t−τ3+a41ty1t−τ4b14t+y1t−τ4+d3tv3t×x4t−y4t+M4∫t−τ3ta44s+τ3ds×a44tx4t−τ3−y4t−τ3+M4∫t−τ3ta44s+τ3ds×a41tb14t+x1t−τ4x1t−τ4−y1t−τ4+M4∫t−τ3ta44s+τ3ds×d3tu3t−v3t.
Let(63)V4t=V41t+V42t+V43t,where(64)V43t=∫t−τ4ta41w+τ4b14w+τ4b14w+τ4+x1wb14w+τ4+y1w×x1w−y1wdw+M4∫t−τ3t∫ww+τ3a44s+τ3a44w+τ3×x4w−y4wdsdw+M4∫t−τ4t∫w+τ4−τ3w+τ4a44s+τ3a41w+τ4b14w+τ4+x1w×x1w−y1wdsdw.
By (62) and (64), we have(65)D+V4t≤−a44l−r4m+a44mM4+a41mM1b14l+M1+d3mN3×∫t−τ3ta44ms+τ3ds−M4a44m∫tt+τ3a44ms+τ3dsx4t−y4t+a41mb14mb14l+m12+M4a41mb14l+m1×∫t+τ4−τ3t+τ4a44ms+τ3dsx1t−y1t+d3m+M4d3m∫t−τ3ta44ms+τ3dsu3t−v3t.
Moreover, defineV5t,V6t,andV7t by(66)V5t=lnu1t−lnv1t,V6t=lnu2t−lnv2t,V7t=lnu3t−lnv3t.
Calculating the upper right derivative ofV5t, V6t, and V7t along the solutions of system (2), it holds that(67)D+V5t≤sgnu1t−v1t−f1tu1t−v1t+q1tx1t−y1t≤−f1ltu1t−v1t+q1mtx1t−y1t,(68)D+V6t≤sgnu2t−v2t−f2tu2t−v2t−q2tx3t−y3t≤−f2ltu2t−v2t+q2mtx3t−y3t,(69)D+V7t≤sgnu3t−v3t−f3tu3t−v3t−q3tx4t−y4t≤−f3ltu3t−v3t+q3mtx4t−y4t.
Finally, we define a Lyapunov function as follows:(70)Vt=∑i=17Vit.
From (51), (58), (65), and (67)–(69), we can get that, for all t≥T+τ,(71)D+Vt≤−∑i=14Aixit−yit−∑j=13Bjujt−vjt.
In view of conditionsH9 of Theorem 2, there exists a constant α>0 and T∗>T+τ such that, for all t≥T∗, it holds that(72)Ai≥α>0,Bj≥α>0,i=1,2,3,4,j=1,2,3.
Integrating fromT∗ to t on both sides of (71) and by (72), we have(73)Vt+α∫T∗t∑i=14xis−yis+∑j=13ujs−vjsds≤VT∗<+∞.
Therefore,Vt is bounded on T∗,+∞, and we have(74)∫T∗∞∑i=14xit−yit+∑j=13ujt−vjtds≤VTα<+∞.
By (74), we also have(75)∑i=14xit−yit+∑j=13ujt−vjt∈L1T,+∞.
Applying Theorem1, we can obtain that xit−yit, ujt−vjt, i=1,2,3,4, j=1,2,3, and their derivatives remain bounded on T∗,+∞, and xit−yit and ujt−vjt, i=1,2,3,4, j=1,2,3, are uniformly continuous on T∗,+∞. By Lemma 5, we can conclude that(76)limt⟶+∞xit−yit=0,i=1,2,3,4,limt⟶+∞ujt−vjt=0,j=1,2,3.
This completes the proof, and the solution of system (2) is globally attractive.Remark 3.
From the result of Theorem2, the feedback control items and delays have influence on the globally attractive of system (2), namely, the population densities for the system will approach to a positive equilibrium when system (2) satisfies initial conditions (7) and assumptions H1−H9.Remark 4.
It should not be difficult to see that the method and skills of constructing a suitable Lyapunov function in the paper are very different with those of [29], and the method of constructing Lyapunov function in [29] cannot be used in this manuscript.


## 4. Periodic Solution

Assuming that coefficients of system (2) are positive continuous and ω−periodic functions, then system (2) is changed to the ω−periodic system. In this section, we shall obtain conditions for the existence, uniqueness, and stability of a positive periodic for system (2) by using the fixed-point theory and some new analysis method. For convenience, we give firstly the following lemma.Lemma 6 (see [34]).
LetS⊂Rn be convex and compact. If mapping T:S⟶S is continuous, then there exists a fixed point, i.e., there exists x∗∈S such that Tx∗=x∗.
Next, we state and prove our main result with respect to the periodic solution for system (2).Theorem 3.
Assume that system (2) is a ω-periodic system and satisfies conditions H1−H9; then, system (2) has a positive unique ω-periodic solution, which is globally asymptotically stable.Proof.
According to the existence and uniqueness theorem of solutions of functional differential equations [35–38], we can define a Poincare mapping T:R+7⟶R+7 as follows:(77)TX0=Xt,ω,X0,where Xt,ω,X0=x1t,x2t,x3t,x4t,u1t,u2t,u3t be a positive solution of system (2) with initial conditions (7). And define(78)S=x1,x2,x3,x4,u1,u2,u3∈R+7mi≤xi≤Mi,i=1,2,3,4nj≤uj≤Nj,j=1,2,3.Then, it is obvious that S⊂R+7 is a convex and compact set. By Theorem 1 and the continuity of solution of system (2) with respect to initial conditions (7), the mapping T:S⟶S is continuous. Furthermore, it is not difficult to show that system (2) has a positive unique ω−periodic solution, which is globally asymptotic stability by using Lemma 6 and Theorem 2.Remark 5.
From the result of Theorem3, the feedback control items and delays have influence on the periodic of ω−periodic system (2). That is, the population densities for the periodic system will tend to a positive ω−periodic equilibrium when system (2) satisfies ω−periodic initial conditions (7) and assumptions H1−H9.


## 5. Numerical Simulation

In this section, we give some numerical simulations to illustrate the feasibility of our main results. In order to facilitate the verification of Theorem3 where the system is a periodic system, we select the periodic function as the coefficient of system (2) and consider the following multidelay periodic predator-prey model with feedback control and prey diffusion:(79)x˙1t=x1t5+cosπt−6.25+1.25sinπtx1t−0.75+0.25sinπtx3t1.8+0.2sinπt+x1t−1+0.5sinπtx4t1.5+0.5sinπt+x1t−0.015+0.005sinπtu1t+0.002+0.001sinπtx2t−x1t,x˙2t=x2t2.5+0.5sinπt−2.25+0.25cosπtx2+0.0015+0.0005sinπtx1t−x2t,x˙3t=x3t−0.15+0.01cosπt−5.5+0.5sinπtx3t−τ1+5.5+0.5sinπtx1t−τ21.8+0.2sinπt+x1t−τ2+0.5+0.1sinπtu2t,x˙4t=x4t−0.025+0.005cosπt−7.25+0.05sinπtx4t−τ3+7.5+0.5sinπtx1t−τ41.5+0.5sinπt+x1t−τ4+0.45+0.05sinπtu3t,u˙1t=1.5+0.5cosπt−1.6+0.4sinπtu1t+0.15+0.05sinπtx1t,u˙2t=2+0.5cosπt−1.05+0.05sinπtu2t−1.5+0.5sinπtx3t,u˙3t=3.55+0.05cosπt−3.45+0.45sinπtu3t−2.25+0.25sinπtx4t,where time delay τ1=0.01, τ2=0.02, τ3=0.03, and τ4=0.04, and the initial conditions are as follows:(80)x1t=sint+1,t∈−0.04,0,x20=1,x3t=sint+1,t∈−0.01,0,x4t=sint+1,t∈−0.03,0,ui0=1,i=1,2,3.From the conditions of Theorem3, we select some specific parameter values shown in models (79) and (80) according to the calculation. Of course, the choice of parameters in the model is not unique. It is not difficult to verify that system (79) satisfies the conditions of Theorems 1–3. It holds from Theorems 1 and 2 that system (79) is permanent and globally attractive. By means of Theorem 3, system (79) has a positive unique periodic solution, which is globally asymptotic stability. By using the software package MATLAB 7.1, we can solve the numerical solutions of systems (79) and (80), as shown in Figures 1–3. Figure 1 shows that the permanence of system (79) with initial conditions (80) and time delay τ1=0.01,τ2=0.02,τ3=0.03,andτ4=0.04. It is not difficult to see from Figure 2 that system (79) is globally attractive. Figure 3 shows the dynamic behavior of system (79).Figure 1
The numerical solution of system (79) with initial condition (60) (80).
(a)
(b)Figure 2
Numerical solution of system (79) with different initial conditions.Figure 3
Dynamic behavior of system (79).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)In general, ecosystems represented by deterministic Lotka–Volterra models are not suitable for describing the real behavior of population dynamics. The feedback control variable we described above is strictly related to the environmental noise effect. It is well known that the environmental noise is a very important factor in modeling natural phenomena in ecological systems. At present, it is impossible to model population dynamics without considering the interaction between the ecological system and environmental noise, always present in nature. In fact, the environmental noise acts constructively, producing counterintuitive effects and interesting noise-induced phenomena, such as stochastic resonance, noise-enhanced stability, and resonant activation (cf. [39–43]). In addition, the study of nonlinear dynamical systems with external noise has found some counterintuitive phenomena that play a constructive role in noise and are highly fundamental and practical interests in many scientific fields (cf. [44–47]). The presence of noise sources can alter the stability of the ecosystem (cf. [48–51]). In addition, this noise in general breaks the symmetry of the ecosystem, see, for example, [52], and moreover, it gives rise to extinction phenomena because the multiplicative noise has the effect of an absorbing barrier, see, for example, [53].As a response to this, we consider the following corresponding stochastic predator-prey model with the multiplicative noise sources:(81)x˙1t=x1t5+cosπt−6.25+1.25sinπtx1t−0.75+0.25sinπtx3t1.8+0.2sinπt+x1t−1+0.5sinπtx4t1.5+0.5sinπt+x1t+Dξ1t+0.002+0.001sinπtx2t−x1t,x˙2t=x2t2.5+0.5sinπt−2.25+0.25cosπtx2t+0.0015+0.0005sinπtx1t−x2t,x˙3t=x3t−0.15+0.01cosπt−5.5+0.5sinπtx3t−τ1+5.5+0.5sinπtx1t−τ21.8+0.2sinπt+x1t−τ2+Dξ2t,x˙4t=x4t−0.025+0.005cosπt−7.25+0.05sinπtx4t−τ3+7.5+0.5sinπtx1t−τ41.5+0.5sinπt+x1t−τ4+Dξ3t,where ξit are δ-correlated Gaussian white noise sources with zero mean, that is, δit=0 and δitδit′=Diδt−t′δij, with Di, the noise intensities (here, we choose all equal noise intensities, that is, Di=D) and time delays τ1=0.01, τ2=0.02, τ3=0.03, and τ4=0.04, and the initial conditions are as follows:(82)x1t=sint+0.5,t∈−0.04,0,x20=0.5,x3t=sint+0.5,t∈−0.01,0,x4t=sint+0.5,t∈−0.03,0.By using MATLAB7.1 software, we can obtain some numerical solutions of systems (81) and (82) with different noise intensities, as shown in Figures 4–6. From Figure 4, it is easy to find that the new stochastic equations (81) and (82) have an almost deterministic periodic solution when the noise intensity is very low. Figure 5 shows that the periodicity of the prey density in the 1-st patch x1t and the predator densities xit,i=3,4 are destroyed when the noise intensities increase to 0.1, and the smoothness of the solution becomes worse as the noise density increases. Moreover, as can be seen from Figure 6, when the noise intensities increase to 2.5, all predator species tend to become extinct. In order to further discover the influence of noise on the property of the solution of the system, we continue to solve numerically equations (81) and (82) by changing the parameters and initial conditions of the system, as shown in Figures 7–11. From Figures 7–9, it is easy to find that the initial value has little effect on the properties of the solution of equations (81) and (82). Figure 10 shows that the noise intensity plays an important role in the periodicity and persistence of the solution. Moreover, as can be seen from Figure 11, time delays play an important role in the stability and periodicity of solutions. By comparing the results shown in Figures 4–11 with those shown in Figures 1–3, it can be seen that stochastic environmental noise has obvious influence on the persistence, periodicity, and stability of the predator-prey system.Figure 4
The numerical solutions of systems (81) and (82) with very low noise intensities.
(a)
(b)
(c)
(d)Figure 5
The numerical solutions of systems (81) and (82) with low noise intensities.
(a)
(b)
(c)
(d)Figure 6
The numerical solutions of systems (81) and (82) with high noise intensities.
(a)
(b)
(c)
(d)Figure 7
The numerical solutions of systems (81) and (82) with very low noise intensities and different initial values.
(a)
(b)
(c)
(d)Figure 8
The numerical solutions of systems (81) and (82) with low noise intensities and different initial values.
(a)
(b)
(c)
(d)Figure 9
The numerical solutions of systems (81) and (82) with high noise intensities and different initial values.
(a)
(b)
(c)
(d)Figure 10
The numerical solutions of systems (81) and (82) with different noise intensities.
(a)
(b)
(c)
(d)Figure 11
The numerical solutions of systems (81) and (82) with different time delays, where τ1=0.01,0.02,0.03,0.04, τ2=0.3,0.03,0.4,0.05, τ3=0.6,0.04,0.7,0.06, τ4=0.8,0.05,0.9,0.07, τ5=1,0.06,1.1,0.07, τ6=1.1,0.08,1.2,0.08, and τ7=1.2,0.09,1.3,0.09.
(a)
(b)
(c)
(d)


## 6. Conclusion

In this paper, we propose and investigate a class of 3-species multidelay ratio-dependent predator-prey model with feedback controls and prey diffusion. By using the theory of delay differential inequalities and developing some new analysis methods as well as constructing a suitable Lyapunov function, some sufficient conditions are obtained to ensure the permanence and global stability of the system. Furthermore, by defining a Poincare mapping and using the Brouwer fixed-point theorem, some conditions for the existence, uniqueness, and stability of positive periodic solutions for the corresponding periodic systems are obtained. In addition, some numerical solutions of the equations describing the system are given to verify the obtained criteria. Finally, the corresponding stochastic predator-prey model with multiplicative noise sources is solved numerically, and some interesting new dynamics are obtained.Time delay and diffusion are very common phenomena in the ecosystem, but some scholars have not studied these problems deeply. The time delay and diffusion terms added in this paper are relatively simple and cannot reflect more general ecosystems, but they are the focus of our future research, including the extension of multiple time delays to infinite time delays. By analyzing the system after adding the feedback control, it can be found that the feedback control item has influence on the persistence and global stability of the original system. Therefore, in some ecosystems, some species can be controlled to maintain the balance and sustainable development of the ecosystem, which is also the practical significance of this project. In addition, by comparing the numerical solutions of equations (79) and (80) with those of equations (81) and (82), it is found that the stochastic environmental noise will affect the persistence, periodicity, and stability of the periodic predator-prey system.Remark 6.
The main contribution and innovation of this paper are as follows: (1) we introduce the control variables, ratio-dependent functions, prey diffusion term, and delays to the known model, and the new model can better describe the interactions among multispecies. (2) To study the new model, we obtain some new methods and skills, such as the structure of the Lyapunov function and the applications of delay differential inequalities, which can also be used to research other related models. (3) In this paper, the research contents are more rich than the related references. We study not only the permanence and global attractivity for the new model but also study the periodic solution. In addition, the corresponding stochastic predator-prey model with multiplicative noise sources is solved numerically, and some interesting new dynamics are obtained. (4) The sufficient conditions obtained herein are new, general, and easily verifiable, which provide flexibility for the application and analysis of the 3-species multidelay Lotka–Volterra ratio-dependent predator-prey model.



---

*Source: 5703859-2020-04-11.xml*