# Integrability and Multiple Limit Cycles in a Predator-Prey System with Fear Effect

Authors: Xinli Li; Ming Zhang
Journal: Journal of Function Spaces
 (2019)

---

## Abstract

In this paper, a class of predator-prey systems with fear effect is investigated, the integrability conditions of the origin and the positive equilibrium are obtained, and the fact that three limit cycles can be bifurcated from the positive equilibrium is proved, so bistable phenomenon can occur for this system.


---

## Body

## 1. Introduction

Predator-prey systems which were considered to be among the most important models in ecological systems have been investigated intensively and extensively since they were proposed; there have been many good results for these systems. The first model for predator-prey problem was Lotka-Volterra model; this system was improved gradually year by year. Many new models were constructed with the development of biological system.Functional response which was introduced by Holling has a great effect on predator-prey system. In 1965, the famous Holling type II functional response of predators [1] was proposed by Holling. As application, the population dynamics of predator-prey systems with the Holling type II functional response have been studied by many authors; for such a model, the result about the existence of a unique stable limit cycle has also been proved. Considering the reality, many more complicated functional responses were suggested one by one, and predator-prey systems have also been modelled gradually. For example, some monotone response functions dependent on prey have been investigated by May [2], Seo and DeAngelis [3], and Huang et al. [4]. Furthermore, Zhu et al. [5], Ruan and Xiao [6], Freedman and Wolkowicz [7], and Wolkowicz [8] studied some nonmonotone response functions dependent on prey. Some interesting bifurcation phenomenon was found for those systems. Functional responses dependent on prey only have been considered intensively; some functional responses dependent on both prey and predators were also considered in recent years, and the famous functional responses dependent on both prey and predators are the Beddington-DeAngelis functional responses, see [9–13]), and ratio dependent functional response [14, 15].Although the functional responses become more and more complicated for predator-prey systems, they can only reflect the direct killing. There is profound lack of the functional responses and many realities which cannot be concluded in those systems. For example, in 2016, Zou and Wang considered a predator-prey model incorporating the cost of fear to explore the impact that fear can have on population dynamics in predator-prey systems, dynamical behaviors of this system was discussed, and a limit cycle which can be bifurcated from the positive equilibrium by Hopf bifurcation was proved [16]. But the number of limit cycles is not considered in their paper. In fact there can be many limit cycles in predator-prey system which can lead to multiple stable phenomenons. For example, the system of Leslie-type predator-prey schemes with a nonmonotonic functional response and Allee effect allows the existence of three limit cycles under certain conditions over the parameters; by Hopf bifurcation [17], the first two cycles are infinitesimal ones generated, and then the third one can be bifurcated from a homoclinic loop. The singular perturbation theory which was developed by Dumortier and Roussarie can be used to study the canard phenomenon. The authors studied the canard phenomenon for predator-prey systems with response functions of Holling types [18] by using this method. For a Gause type predator-prey system with Holling type III functional response and Allee effect on prey, multiple limit cycles were considered in [19] by multiple Hopf bifurcations. In [20], a criterion to calculate the multiplicity of a multiple focus was developed by the authors for general predator-prey systems. Recently, the number of limit cycles in a Leslie-Gower-type predator-prey system with weak Allee effect on prey was considered in [21].The number of limit cycles is closely related not only to the Hilbert 16th problem but also to bistable phenomenon or multiple stable phenomenon. For instance, in [21], the existence conditions for which three limit cycles appear and surround different equilibrium points were proved. Above all, the three limit cycles have different stability which means the phenomenon of multistability. For finitely smooth planar autonomous differential systems, in [22], the author considered the bifurcation theory. Theory of rotated equations is discussed and applied to a population model in [23]. For a family of real planar polynomial ODE systems depending on parameters, the problem of how to determine the systems in the family which become time-reversible after some affine transformation was considered in [24]. Recently, an improvement on the number of limit cycles bifurcating from a nondegenerate center of homogeneous polynomial systems was obtained in [25] by normal form method. Some new perturbation methods were given in [26–28]. The integrability of equivalent systems was also considered in [29, 30]. Therefore, in this paper, we reinvestigated the predator-prey system with fear effect which is written as follows:(1)dudt=r0u1+kv-du-au2-puv1+qu,dvdt=puv1+qu+mv.By transformation(2)u¯=cp-mqmu,v¯=kv, and time scale(3)dt=1+qu1+kvmdt¯, system (1) can be reduced to(4)dudt=ua1+a2u-a3v-a4uv-a5u2-a6v2-a5u2v,dvdt=vu-11+v.And system (4) has three singular points(5)0,0,a2+a22+4a1a52a5,0,1,v¯1 where(6)v¯1=-a3+a4+a5-a3+a4+a52+4a1+a2-a5a62a6.The classification of the three singular points and stability have been discussed in [16], but the integrability and center problem remained. The rest of the paper is organized as follows. In the next section, we prove that the origin of system (4) is integral. In Section 3, the first four Lyapunov constants will be computed; bifurcation of limit cycles and center conditions of (4) are investigated. Section 4 is devoted to discussing the kind of the boundary equilibrium. At last, a conclusion is drawn to illustrate that bistable phenomenon can occur in this system.


## 2. Integrability Conditions of(0,0)

For the integrability of origin of system (4), the authors have proved that it is a saddle. In this section, we will compute the saddle values; furthermore, the integrability conditions can be obtained at the same time. We only consider a special case; namely, let a1=1; we study the integrability conditions theoretically. The system (4) becomes(7)dudt=u1+a2u-a3v-a4uv-a5u2-a6v2-a5u2v,dvdt=vu-11+v.Theorem 1.
The first five saddle values at the origin of system (7) are(8)μ1=a2a3-a4,μ2=-a5-3a3a5-2a32a5+a2a6-a5a6, If 1+3a3+2a32+a6≠0(9)μ3=-a2a61+3a3+2a32+a6-2-5a3+2a2a3-3a32+5a2a32+3a2a33+8a6+a2a6+9a3a6+a2a3a6, (1) (1+a3)(2a3+3a32+a6)≠0(10)μ4=-a622+5a3+3a32-8a6-9a3a631+a321+3a3+2a32+a62a3+3a32+a63h1,μ5=a622+5a3+3a32-8a6-9a3a6241+a331+3a3+2a32+a622a3+3a32+a64h2,μ6=-a622+5a3+3a32-8a6-9a3a61201+a341+3a3+2a32+a622a3+3a32+a65h3; (2) a3=-1(11)μ3=a2a62; (3) a6=-2a3-3a32(12)μ3=a2a32+3a322-1-a3+a321+9a3+9a32,μ4=-a2a32+3a3243-1-a3+a3272+74a2+a22-144a3+36a2a3,μ5=a2a32+3a348-1-a3+a322163092+310204a2+142507a22+2052a23+259668a3+444930a2a3+179790a22a3+2752a23a3. If a6=-(1+3a3+2a32), (13)μ2=-a21+a31+2a3, (4) a2=0(14)μ3=1+a325+9a3a5,μ4=-131+a32274+717a3+407a32a5; (5) a3=-1/2(15)μ3=182+a21+a3a5. Here h1,h2,h3 are given in the appendix.Thus, we can get the following integrability conditions of system (7) easily.Theorem 2.
The origin of system (7) is integrability if and only if one of the following conditions holds: (16)Ia2=a4=a5=0;IIa4=a2a3,a5=a6=0;IIIa4=-a2,a3=-1,a6=0;IVa2=-2,a3=-12,a4=1,a6=0.Proof.
When condition (I) holds, system (7) becomes(17)dudt=u1-a3v-a6v2,dvdt=vu-11+v, which has an integrating factor(18)μ1u,v=1uv1+v.
When condition (II) holds, system (7) can be rewritten as(19)dudt=u1+a2u-a3v-a2a3uv,dvdt=vu-11+v,which admits an integrating factor(20)1-1+uuv1+v,a2=-1,uv-1-2a2/1+a21+v-1+a2a3/1+a2,a2≠-1.
When condition (III) holds, system (7) can be changed into(21)dudt=u1+a2u+v+a2uv-a5u2-a5u2v,dvdt=vu-11+v,which admits an integrating factor(22)μ2u,v=1u-1-a2u+a5u2v1+v.
When condition (IV) holds, system (7) can be transformed into(23)dudt=u1-2u+12v-uv-a5u2-a5u2v,dvdt=vu-11+v,which has an integrating factor(24)μ3u,v=1u3v3.


## 3. Hopf Bifurcation at(1,v¯1)

According to results in [16], there exists a positive singular point (1,v¯1). Let x=u-1,y=v-v¯1; system (4) can be transformed into(25)dxdt=J11x+J12y+f1x,y,dydt=J21x+J22y+f2x,y. Here (26)J11=a2-2a5-a4v1-2a5v1;J12=-a3-a4-a5-2a6v1;J21=v11+v1;J22=0. Suppose there are two characteristic values λ1=i,λ2=-i; then(27)λ1+λ2=J11=0,λ1λ2=J12J21=1, namely,(28)a2-2a5-a4v1-2a5v1=0, and we have(29)v¯1=a2-2a5a4+2a5. Without loss of generality, let(30)a4=-2a5+β,a2=2a5+αβ, and the singular point (1,v¯1) becomes (1,α). From the first equation J11=0, we get(31)a6=a1+a5-a3α+a5αα2.Meanwhile, the second equation J12J21=1 means that(32)-a3-a4-a5-2a6v1v11+v1=1,which yields(33)a3=2a1+2a5+2a1α+3a5α+a5α2+αβ+α2β-1α1+α.Accordingly, system (7) can be changed into(34)dxdt=-1+uα21+α-a5u2α21+α2-vα1+uα1+αa5u+β+v2-1+a11+α+1+αa5+αβ,dydt=uv+α1+v+α.The Jacobin matrix at the(1,α) of system (34) can be written as follows:(35)J0,01=0-1α1+αα1+α0, and by linear transformation(36)x=v,y=-α1+αu system can be shifted to(37)dudt=-vu+α1+u+αα1+α,dvdt=--v+α+α2α21+α-u2+a1u2+a5u2-a5v2-a5uv2-uα-u2α+2a1u2α+2a5u2α-a5v2α-uα2+a1u2α2+a5u2α2+u2αβ+uvαβ+2u2α2β+uvα2β+u2α3β.The Jacobin matrix at the origin of system (37) can be written as follows:(38)J0,02=0-110,and, with the help of computer algebraic system Maple, we can compute the first five Lyapunov constants for system (37) by using formal series method which was given in the following.Theorem 3.
For system (37), the first five Lyapunov constants at the origin can be given as follows:(39)μ1=i4α21+αa52+α21+αβ+α1+αβ-1+1+αa1+αβ,μ2=ia5-2+αβ+α2β24α61+α3β2320a52-6a5αβ-4a52αβ-10a5α2β-4a52α2β+α3β2+a5α3β2+α4β2+a5α4β2-α3β3-α4β3,μ3=-ia5-2+αβ+α2β1536α71+α5β-5+α1+αβ3a5-24012+5α5+α-40α1+α-18+α79+57α+30α2β+4α21+α2480+α1126+α426+205αβ2-2α31+α396+α216+α54+73αβ3+2α41+α4-12+α-20+α3+10αβ4-α61+α5-4+α2β5+α21+α-1+ββ-480+α1+α200+β720+α-420+α-360+1+αβ204+154α-21+α6+α2+7αβ+α31+α2β2,μ4=ia5-1+β-2+αβ+α2β46080α61+α5-5+αβ+α2β2h6h4,μ5=ia5-1+β-2+αβ+α2β13271040α81+α7-5+αβ+α2β4h6h5, where h4,h5,h6 are given in the appendix.Theorem3 yields the following.Proposition 4.
For system (37), the first five Lyapunov constants at the origin are zero if and only if one of the following conditions holds:(40)a1=--1+αβ+α2β1+α,a5=0;(41)a1=-1+a5+a5α1+α,β=2α1+α.Furthermore, we have the following.Theorem 5.
The origin of system (37) is a center if and only if one of the conditions in Proposition 4 holds.Proof.
When condition (40) holds, system (37) becomes(42)dudt=-vu+α1+u+αα1+α,dvdt=-u-v+α+α2-1+vβα1+α,which has four invariant algebraic curves (43)L1=u+α,L2=u+1+α,L3=-v+α+α2,L4=-1+vβ, and admits an integrating factor(44)μ4x,y=L1-1+αβ-1+α1+αβL2-1-1+αβ-1+α1+αβL3-1-α1+αβ.
When condition (41) holds, system is simplified into(45)dudt=-vu+α1+u+αα1+α,dvdt=-v+α+α2-2uv+a5v2+a5uv2+αu+a5αv2+α2uα21+α2,which has three invariant algebraic curves L1,L2,L3 and admits an integrating factor(46)μ5x,y=L11+2a5+2/1+αL2-α+2/αL3-3.Two examples of phase plane of system (42) when condition (40) or (41) is satisfied are given in Figures 1 and 2.Figure 1
The phase plane of system (42) when α=β=1,a1=-1/2,a5=0.Figure 2
The phase plane of system (45) when α=β=1,a1=-3/2,a5=1.Based on Theorem3, we can conclude the following.Theorem 6.
The origin of system (37) is fifth order weak foci if and only if one of the following conditions holds: (47)H1:a1=-0.2282,a5=0.337573,α=-9.39076168286412,β=0.02411227315630107;H2:a1=-1.57238,a5=0.688478,α=4.0444054384945956,β=0.28989931115360074;H3:a1=2.5,a5=0,α=-0.5000000000000074,β=0.9999999999999943.Proof.
Letμ1=μ2=μ3=0, and we obtain(48)a1=2a5-αβ-α2β+a5α2β+a5α3β+α2β2+2α3β2+α4β2α1+α2β;a52=-6a5αβ-10a5α2β+α3β2+a5α3β2+α4β2+a5α4β2-α3β3-α4β34-5+αβ+α2β;a5=α21+α-1+ββh6h4; submitting a5 into a52, we have(49)h7=-28802+β+α-80102+α74+5α5+α+80-249+α-412+α-216-66α+α3β+41+α-1260+α150+α-2188+α1576+α1295+α350+69αβ2-4α1+α2-1620+α-2730+α-3805+2α-526+α-79+α57+20αβ3+2α21+α3684+α522+α244+α-468+α-342+α24+19αβ4-2α31+α4-36+α48+α110+α-74+α-58+α5+3αβ5-α51+α512+14α+α3β6+α71+α6β. Solve(50)h4=0,h7=0 and h5≠0, and we can obtain the results.Considering the reality, all parameters are positive in a predator-prey system, so we cannot get five limit cycles for this system.In fact,h5=0 means that(51)α=-110742+151145β106746, and, submitting it into a1,a5, we get Figure 3 which shows that a1,a5 cannot be larger than zero simultaneously.Figure 3
a 1 , a 5 with β.We can adjust the parameters to satisfyμ1=μ2=0,μ3≠0, and all parameters are positive. For example,(52)a1=0.380667,a2=0.268514,a3=0.735621,a4=0.481486,a5=0.009,a6=0.10697. Then three limit cycles can be bifurcated from the origin; namely, bistable phenomenon can occur for this system; this is not discussed in previous discussion.


## 4. Boundary Singular Point

In this section, we consider the boundary singular point of system (4). The boundary singular point is(53)u1,0=a2+a22+4a1a52a5,0, and the Jacobin matrix at (u1,0) of system (4) is (54)J11=-2a1-a2u1,J12=-u1a3+u1a4+a5u1,J21=0,J22=u1-1. It is easy to testify that when u1>1, system (4) has two eigenvalues with different signs which means (u1,0) is a saddle. When u1<1, system (4) has two eigenvalues with negative signs which means (u1,0) is a stable node. When u1=1, namely, a5=a1+a2(u1,0), system (4) becomes (1,0) and one eigenvalue is zero. System (4) becomes(55)dudt=ua1+a2u-a1+a2u2-a3v-a4uv-a1+a2u2v-a6v2,dvdt=vu-11+v.Shifting the (1,0) of system (4) to origin, system (4) becomes(56)dudt=1+u-va3+a4+a4u+a6v-a21+uu+v+uv-a1v+u2+u1+v,dvdt=vu1+v.the Jacobin matrix at origin system (56) is (57)J11=-2a1-a2,J12=-a1-a2-a3-a4,J21=0,J22=0. By transformation(58)x=v,y=u+a1+a2+a3+a42a1+a2v, and time scale(59)τ=-2a1+a2t, system (56) becomes(60)dxdσ=x1+xa1x+a2x+a3x+a4x-2a1y-a2y2a1+a22=ϕx,y,dydσ=y+α2x2+α3x3+⋯=ψx,y.Based on the theorem in [31], we have the following.Theorem 7.
The origin of system (56) is a saddle-node singular point.Proof.
Fromψ(x,y)=0, we get a solution y=φ(x) by implicit function theorem; submitting it into ϕ(x,y), we have(61)ϕx,y=a1+a2+a3+a42a1+a22x2+h.o.t. Based on the theorem in [31], the origin of system (56) is a saddle-node singular point.An numerical simulation of phase portrait of system (56) when a1=a2=a3=a4=1 is given in Figure 4 which shows that the origin of system (56) is a saddle-node singular point.Figure 4
Whena1=a2=a3=a4=1, the (1,0) is a saddle node.


## 5. Conclusion

In this paper, a class of predator-prey systems with fear effect is investigated; the integrability conditions of(0,0) and positive equilibrium are obtained and proved by finding their integrating factors. Furthermore, we compute the first five focal values at the positive equilibrium; considering the reality, when we choose all parameters to be positive, the fact that at least three limit cycles can be bifurcated from the positive equilibrium is proved. Consequently, bistable phenomenon can occur for this system.



---

*Source: 3948621-2019-06-27.xml*