# Dynamical Properties of a Delay Prey-Predator Model with Disease in the Prey Species Only

Authors: Xiangyun Shi; Xueyong Zhou; Xinyu Song
Journal: Discrete Dynamics in Nature and Society
 (2010)

---

## Abstract

A three-dimensional ecoepidemiological model with delay is considered. We first investigate the existence and stability of the equilibria. We then study the effect of the time delay on the stability of the positive equilibrium. The existence of a Hopf bifurcation at the positive equilibrium is obtained through the study of an
exponential polynomial equation with delay-dependent coefficients. Numerical simulation with a hypothetical set of data has been carried out to support the analytical findings.


---

## Body

## 1. Introduction

Recently, epidemiological models have become important tools in analyzing the spread and control of infectious diseases after the seminal model of Kermack-Mckendrick on SIRS (susceptible-infected-removed-susceptible) systems [1], and numerous mathematical models were developed to study a disease transmission, to evaluate the spread of epidemics, to understand the mechanisms of epidemics in order to prevent them or minimise the transmission of diseases and other measures. After the seminal models of Vito Volterra and Alfred James Lotka in the mid 1920s for predator-prey interactions [2, 3], predator-prey models have received a great deal of interest in different fields of science, from evolutionary and behavioural ecology to conservation and population biology and even economics [4–8].It is well known that species do not exist alone in the natural world. It is of more biological significance to study the persistence-extinction threshold of each population in systems of two or more interacting species subjected to disease. Few investigators [9–12] have paid attention to study prey-predator model with infection.In order to study the influence of disease on an environment where two or more interacting species are present. In this paper, we will put emphasis on such an ecoepidemiological system consisting of three species, namely, the sound prey (which is susceptible), the infected prey (which becomes infective by some viruses), and the predator population.We have two populations:(1)
the prey, whose total population density is denoted byN(t),(2)
the predator, whose population density is denoted byy(t).We make the following assumptions.(A1)
In the absence of infection and predation, the prey population density grows logistically with carrying capacityK(K>0) and an intrinsic birth rate constant r(r>0) [13, 14],
(1.1)dS(t)dt=rS(t)(1-S(t)K).(A2)
In the presence of disease, the total prey populationN(t) are divided into two distinct classes, namely, susceptible populations, S(t), and infected populations, I(t) [13, 14]. Therefore, at any time t, the total density of prey population is
(1.2)N(t)=S(t)+I(t).(A3)
We assume that only susceptible preyS(t) are capable of reproducing with logistic law (1.1); that is, the infected prey I(t) are removed by death (say its death rate is a positive constant μ), or by predation before having the possibility of reproducing. However, the infective population I(t) still contributes with S(t) to population growth toward the carrying capacity [13, 14].(A4)
We assume that the force of infection at timet is given by βe-μ1τS(t)I(t-τ), where β is the average number of contacts per infective per day, μ1 is disease-induced death rate and τ>0 is a fixed time during which the infectious agents develop in the vector and it is only after that time that the infected vector can infect a susceptible prey [15]. Hence, the SI model of the infected prey is
(1.3)dS(t)dt=rS(t)(1-S(t)K)-βe-μ1τS(t)I(t-τ),dI(t)dt=βe-μ1τS(t)I(t-τ)-μI(t).(A5)
Numerous field studies show that infected prey are more vulnerable to predation compared with their noninfected counterpart  [16–18]. Lafferty and Morris  [17] quantified that the predation rates on infected prey may be 31 times higher compared to that on susceptible prey. Thus, we consider the case when the predator mainly eats the infected prey. We assume that the predator eats only the infected prey with Leslie-Gower ratio-dependent schemes  [19–23]. That is to say, the predator consumes the prey according to the ratio-dependent functional response  [24–26] and the predator grows logistically with intrinsic growth rate δ and carrying capacity proportional to the prey populations size I(t).From the above assumptions we have the following model:(1.4)dS(t)dt=rS(t)(1-S(t)K)-βe-μ1τS(t)I(t-τ),dI(t)dt=βe-μ1τS(t)I(t-τ)-cy(t)I(t)my(t)+I(t)-μI(t),dy(t)dt=δy(t)(1-hy(t)I(t)).All the parameters of system (1.4) are positive. The initial conditions for system (1.4) defined in the Banach space (1.5)C+={φ∈C([-τ,0],R+3)∣φ1(θ)=S(θ),φ2(θ)=I(θ),φ3(θ)=y(θ)},
where ℝ+3={(S,I,y)∈ℝ3∣S≥0,I≥0,y≥0}. We assume that φi(0)>0(i=1,2,3) by the biological meaning.It is well known by the fundamental theorem of functional differential equations that system (1.4) has a unique solution (S(t),I(t),y(t)) satisfying initial conditions (1.5).In this paper, we will discuss the dynamical properties (e.g., the stability of equilibrium, the existence of Hopf bifurcation around the positive equilibrium) of the delay ecoepidemic model (1.4). This paper is organized as follows. Section 2 gives positivity and boundedness of solutions of system (1.4). Section 3 gives the analysis of existence and uniqueness of positive equilibrium as well as its stability, by studying the characteristic equation associated with system (1.4) which takes the form of a third-degree exponential polynomial with delay-dependent coefficients. In Section 4, some numerical simulations will be given to illustrate our theoretical results. Section 5 discusses the results.


## 2. Positivity and Boundedness of Solutions

It is important to show positivity and boundedness for the system (1.4) because the variables S, I and y represent populations. Positivity implies that the populations survive and boundedness may be interpreted as a natural restriction to growth as a consequence of limited resources.Let us put system (1.4) in a vector form by setting X=col(S,I,y)∈ℝ3, (2.1)G(X)=(G1(X)G2(X)G3(X))=(rS(1-SK)-βe-μ1τSI(t-τ)βe-μ1τSI(t-τ)-cyImy+I-μIδy(1-hyI)),
where G:𝒞+→ℝ3 and G∈𝒞∞(ℝ3). Then system (1.4) becomes (2.2)Ẋ=G(X),
with X(θ)=(ϕ1(θ),ϕ2(θ),ϕ3(θ))∈𝒞+ and ϕi(0)>0(i=1,2,3). It is easy to check in (2.1) that whenever choosing X(θ)∈𝒞+ such that xi=0, then Gi(X)|xi(t)=0,X∈ℝ+3≥0, (i=1,2,3), where x1=S(t), x2=I(t) and x3=y(t).Due to lemma in [27] any solutions of (2.2) with X(θ)∈𝒞+, that is, X(t)=X(t,X(θ)), is such that X(θ)∈ℝ3 for all t>0.Because the analysis of the boundedness of solutions of system (1.4) is similar to the one in Section 2 in  [28], we omit it. Let M≥0, and we can get that all solutions of system (1.4) with initial condition enter the region B={(S(t),I(t),y(t))∣0≤S(t),I(t),y(t)≤M}.


## 3. Stability Analysis and Hopf Bifurcation

In this section, we focus on investigating the stability of the equilibria and Hopf bifurcation around the positive equilibrium of the system (1.4). System (1.4) has the boundary equilibrium E1(μ/βe-μ1τ,(r/βe-μ1τ)(1-μ/βe-μ1τK),0)≜(S1,I1,y1) and the positive equilibrium E2(S2,I2,y2), where S2=(μh+c+μm)/βe-μ1τ(m+h), I2=r(βe-μ1τmK+βe-μ1τhK-μh-c-μm)/(βe-μ1τ)2K(m+h), y2=r(βe-μ1τmK+βe-μ1τhK-μh-c-μm)/(βe-μ1τ)2hK(m+h). Clearly, if 1-μ/βe-μ1τK>0, then E1 exists and remains positive, and E2 exists and remains positive if β>(1/Ke-μ1τ)(μ+c/(m+h))≜β0.LetE*(S*,I*,y*) be any arbitrary equilibrium. Then the characteristic equation about E* is given by (3.1)|r-2rS*K-βe-μ1τI*-λ-βe-μ1τS*e-λτ0βe-μ1τI*βe-μ1τS*e-λτ-cmy*2(my*+I*)2-μ-λ-cI*2(my*+I*)20δhy*2I*2δ-2δhy*I*-λ|=0.For equilibriumE1, (3.1) reduces to (3.2)|-rS1K-βe-μ1τI1-λ-βe-μ1τS1e-λτ0βe-μ1τI1βe-μ1τS1e-λτ-μ-λ-c00δ-λ|=0.It is easy to see that the equilibriumE1 is unstable.For equilibriumE2, (3.1) reduces to (3.3)P(λ,τ)+Q(λ,τ)e-λτ=0
with (3.4)P(λ,τ)=λ3+A1(τ)λ2+A2(τ)λ+A3(τ),Q(λ,τ)=B1(τ)λ2+B2(τ)λ+B3(τ),
where (3.5)A1(τ)=δ+μ+cmy22(my2+I2)2+rS2K,A2(τ)=δ(μ+cmy22(my2+I2)2)+rS2K(δ+μ+cmy22(my2+I2)2)+δhcI22(my2+I2)2,A3(τ)=rS2K[δ(μ+cmy22(my2+I2)2)+δhcI22(my2+I2)2],B1(τ)=-βe-μ1τS2,B2(τ)=-δβe-μ1τS2-rβe-μ1τS22K+(βe-μ1τ)2S2I2,B3(τ)=-rδβe-μ1τS22K+δ(βe-μ1τ)2S2I2.
For τ=0, the transcendental equation (3.3) reduces to (3.4) (3.6)λ3+(A1(0)+B1(0))λ2+(A2(0)+B2(0))λ+A3(0)+B3(0)=0.
We can easily get (3.7)A1(0)+B1(0)=rS2K+δ-δ*=rβ̅0β+δ-δ*>0,A2(0)+B2(0)=rS2K(δ-δ*)+β2S2I2=r(μ+cm+h)(1+δ-δ*Kβ-β̅0β),A3(0)+B3(0)=δβ2S2I2=rδ(μ+cm+h)(1-β̅0β)>0
if δ>δ* and β>β0 (obviously, β0>β̅0), where δ*=ch/(m+h)2 and β̅0=(1/K)(μ+c/(m+h)).By Routh-Hurwitz criterion, we can get the following theorem.Theorem 3.1.
Whenτ=0, the positive equilibrium E2 is locally asymptotically stable provided that the conditions(H1)
rβ0>δ*β,(H2)
δ>δ*,(H3)
β>β0.In the following, we investigate the existence of purely imaginary rootsλ=iω(ω>0) to (3.3). Equation (3.3) takes the form of a third-degree exponential polynomial in λ, which all the coefficients of P and Q depending on τ. Beretta and Kuang [29] established a geometrical criterion which gives the existence of purely imaginary of a characteristic equation with delay-dependent coefficients.In order to apply the criterion due to Beretta and Kuang [29], we need to verify the following properties for all τ∈[0,τmax), where τmax is the maximum value which E2 exists. We can easily obtain that τmax=(1/μ1)ln(Kβ(m+h)/μ(m+h)+c).(a)
P(0,τ)+Q(0,τ)≠0;(b)
P(iω,τ)+Q(iω,τ)≠0;(c)
limsup{|P(λ,τ)/Q(λ,τ)|:|λ|→∞,Reλ≥0}<1;(d)
F(ω,τ)=|P(iω,τ)|2-|Q(iω,τ)|2 has a finite number of zeros;(e)
Each positive rootω(τ) of F(ω,τ)=0 is continuous and differentiable in τ whenever it exists.Here,P(λ,τ) and Q(λ,τ) are defined as in (3.4).Letτ∈[0,τmax). It is easy to see that
(3.8)P(0,τ)+Q(0,τ)=A3(τ)+B3(τ)>0.
This implies that (a) is satisfied, and (b) is obviously true because (3.9)P(iω,τ)+Q(iω,τ)=-iω3-A1(τ)ω2+iA2(τ)ω+A3(τ)-B1(τ)ω2+iB2(τ)ω+B3(τ)=[-(A1(τ)+B1(τ))ω2+A3(τ)+B3(τ)]+i[-ω3+(A2(τ)+B2(τ))ω]≠0.From (3.4) we know that (3.10)lim⁡|λ|→+∞|Q(λ,τ)P(λ,τ)|=0.
Therefore (c) follows.LetF be defined as in (d). From (3.11)|P(iω,τ)|2=(-ω3+A2(τ)ω)2+(-A1(τ)ω2+A3(τ))2=ω6+(A12(τ)-2A2(τ))ω4+(A22(τ)-2A1(τ)A3(τ))ω2+A32(τ),|Q(iω,τ)|2=B12(τ)ω4+[b22(τ)-2B1(τ)B3(τ)]ω2+B32(τ),
we have (3.12)F(ω,τ)=ω6+a1(τ)ω4+a2(τ)ω2+a3(τ),
where (3.13)a1(τ)=A12(τ)-2A2(τ)-B12(τ)=δ2+μ2+c2m2(m+h)2+(rS2K)2+2μcm(m+h)2-2δch(m+h)2-β2(e-μ1τ)2S22,a2(τ)=A22(τ)+2B1(τ)B3(τ)-2A1(τ)A2(τ)-B22(τ)=[δ(μ+cmy22(my2+I2)2)+rS2K(δ+μ+cmy22(my2+I2)2)+δhcI22(my2+I2)2]×[δ(μ+cmy22(my2+I2)2)+rS2K(δ+μ+cmy22(my2+I2)2)+δhcI22(my2+I2)2-2(δ+μ+cmy22(my2+I2)2+rS2K)]-[δβe-μ1τS2]2-[rβe-μ1τS22K-(βe-μ1τ)2S2I2]2,a3(τ)=A32(τ)-B32(τ)=(rS2K)2[δ(μ+cmy22(my2+I2)2)+δhcI22(my2+I2)2]2-[δ(βe-μ1τ)2S2I2-rδβe-μ1τS22K]2.It is obvious that property (d) is satisfied, and by Implicit Function Theorem, (e) is also satisfied.Now letλ=iω(ω>0) be a root of (3.3). Substituting it into (3.3) and separating the real and imaginary parts yields (3.14)-A3(τ)+A1(τ)ω2=B2(τ)ωsinωτ+(B3(τ)-B1(τ)ω2)ωcos⁡ωτ,ω3-A2(τ)ω=B2(τ)ωcos⁡ωτ+(B3(τ)-B1(τ)ω2)ωsinωτ.From (3.14) it follows that(3.15a)sinωτ=B1(τ)ω5+[A1(τ)B2(τ)-A2(τ)B1(τ)-B3(τ)]ω3B12(τ)ω4(B22(τ)-2B1(τ)B3(τ))ω2+B32(τ)+(A2(τ)B3(τ)-A3(τ)B3(τ))ωB12(τ)ω4(B22(τ)-2B1(τ)B3(τ))ω2+B32(τ),(3.15b)cos⁡ωτ=[B2(τ)-A1(τ)B1(τ)]ω4+[A1(τ)B3(τ)+A3(τ)B1(τ)-A2(τ)B2(τ)]ω2B12(τ)ω4(B22(τ)-2B1(τ)B3(τ))ω2+B32(τ)-A3(τ)B3(τ)B12(τ)ω4(B22(τ)-2B1(τ)B3(τ))ω2+B32(τ).By the definitions ofP(λ,τ), Q(λ,τ) as in (3.3), and applying the property (a), (3.7) can be written as (3.16)sinωτ=Im⁡P(iω,τ)Q(iω,τ),cos⁡ωτ=-ReP(iω,τ)Q(iω,τ),
which yields (3.17)|P(iω,τ)|2=|Q(iω,τ)|2.Assume that𝒥⊆ℝ+0={ζ∣ζ>0} is the set where ω(τ) is a positive root of (3.18)F(ω,τ)=|P(iω,τ)|2-|Q(iω,τ)|2
and for τ∉𝒥, ω(τ) is not definite. Then for all τ in 𝒥, ω(τ) satisfied (3.19)F(ω,τ)=0.
The polynomial function F can be written as F(ω,τ)=f(ω2,τ), where f is a third-degree polynomial, defined by (3.20)f(z,τ)=z3+a1z2+a2z+a3.Depending on the determinant of equation(3.21)f(z,τ)=z3+a1z2+a2z+a3=0,M=(q/2)2+(p/3)3, where p=a2-(1/3)a12, q=(2/27)a13-(1/3)a1a2+a3, there are three cases for the solutions of (3.21).(i)
IfM>0, (3.21) has a real root and a pair of conjugate complex roots. The real root is positive and is given by
(3.22)ν1=-q2+M3+-q2-M3-13a1.(ii)
IfM=0, (3.21) has three real roots, of which two are equal. In particular, if a1>0, there exists only one positive root, ν1=2-q/23-a1/3; If a1<0, there exists a positive root ν1=2-q/23-a1/3 for -q/23>-a1/3, and there exist three positive roots for a1/6<-q/23<-a1/3, ν1=2-q/23-a1/3, ν2=μ3=--q/23-a1/3.(iii)
IfM<0, there are three distinct real roots, ν1=2(|p|/3)cos(φ/3)-a1/3, ν2=2(|p|/3)cos(φ/3+2π/3)-a1/3, ν3=2(|p|/3)cos(φ/3+4π/3)-a1/3, where cosφ=-q/2(|p|/3)3. Furthermore, if a1>0, there exists only one positive root. Otherwise, if a1<0, there may exist either one or three positive real roots. If there is only one positive real root, it is equal to max(ν1,ν2,ν3).Clearly, the number of positive real roots of (3.21) depends on the sign of a1. When a1≥0, (3.21) has only one positive real root. Otherwise, there may exist three positive real roots.It is easy to know that(3.23)a1=A12(τ)-2A2(τ)-B12(τ)=δ2+μ2+c2m2(m+h)2+(rS2K)2+2μcm(m+h)2-2δch(m+h)2-β2(e-μ1τ)2S22>0if condition(H4):
(3.24)δ2+μ2+c2m2(m+h)2+(rS2K)2+2μcm(m+h)2>2δch(m+h)2+β2(e-μ1τ)2S22holds.Hence, (3.21) has only one positive real root if condition (H4) holds. We denote by z+ this positive real root. Thus, (3.21) has only one positive real root ω=z+, and the critical values of τ and ω(τ) are impossible to solve explicitly, so we will use the procedure described in Beretta and Kuang [29]. According to this procedure, we define θ(τ)∈[0,2π) such that sinθ(τ) and cosθ(τ) are given by the right-hand sides of (3.15a) and (3.15b), respectively, with θ(τ) given by (3.21). This define θ(τ) in a form suitable for numerical evaluation using standard software.And the relation between the argumentθ(τ) and ω(τ)τ in (3.19) for τ>0 must be (3.25)ω(τ)τ=θ(τ)+2nπ,n=0,1,2,….
Hence we can define the maps: τn:𝒥→ℝ+0 given by (3.26)τn(τ):=θ(τ)+2nπω(τ),τn>0,n=0,1,2,…,
where a positive root ω(τ) of (3.19) exists in 𝒥.Let us introduce the functionsMn(τ):𝒥→ℝ, (3.27)Mn(τ)=τ-θ(τ)+2nπω(τ),n=0,1,2,….
that are continuous and differentiable in τ. Thus, we give the following theorem which is due to Beretta and Kuang [29].Theorem 3.2.
Supposeβ>β0 and condition (H4) hold. Assume that ω(τ) is a positive root of (3.4) defined for τ∈𝒥, 𝒥⊆ℝ+0, and at some τ*∈𝒥, Mn(τ*)=0 for some n∈ℕ0={1,2,3,…}. Then a pair of simple conjugate pure imaginary roots λ=±iω exists at τ=τ* which crosses the imaginary axis from left to right if δ(τ*)>0 and crosses the imaginary axis from right to left if δ(τ*)<0, where
(3.28)δ(τ*)=sign⁡{∂F∂ω(ωτ*,τ*)}sign⁡{dMn(τ)dτ|τ=τ*}.
Since(3.29)∂F∂ω(ω,τ)=2ω∂f∂z(ω2,τ),
condition (3.28) is equivalent to
(3.30)sign⁡{dReλdτ|λ=iω(τ*)}=sign⁡{∂f∂z(ω2(τ*),τ*)}sign⁡{dMn(τ)dτ|τ=τ*}.
We can easily observe thatMn(0)<0. Moreover, for all τ∈𝒥, Mn(τ)>Mn+1(τ), when n∈ℕ. Therefore, if M0 has no root in 𝒥, then the Mn functions have no root in 𝒥 and, if the function Mn(τ) has positive roots τ∈𝒥 for some n∈ℕ, there exists at least one root satisfying(3.31)dMndτ(τ)>0.Using Theorem3.1, we can conclude the existence of a Hopf bifurcation as stated in the following theorem.Theorem 3.3.
Assume the following.(1)
If the functionM0(τ) has no positive root in 𝒥, then the positive equilibrium E2 is locally asymptotically stable for all τ≥0.(2)
If the functionM0(τ) has at least one positive root in 𝒥, then there exists τ*∈𝒥, such that the positive equilibrium E2 is locally asymptotically stable for 0≤τ<τ* and becomes unstable for τ≥τ*, with a Hopf bifurcation occurring when τ=τ*, if and only if
(3.32)sign⁡{∂f∂z(ω2,τ)}>0.


## 4. Numerical Simulations

To determine whether it is possible that there are parameter values such that a stability switch is possible, we choose a set of parameters as follows:r=0.1, K=500, β=0.001, μ=0.03, c=8, δ=0.2, m=150, h=0.5, μ1=0.04.Using Maple 10, we compute the coefficientsa1(τ), a2(τ), a3(τ) for the above values. We can get a1(τ)=0.2829795477+0.004706297728e0.04τ>0. Hence the polynomial function f, defined in (3.20), has only one positive real root.The functionM0 is drawn for τ∈𝒥 in Figure 1. One can see that there are two critical values of the delay τ for which stability switches occur, that is, τ*=4.69 and τ**=13.78. In particulary from Theorem 3.3, a Hopf bifurcation occurs when τ is approximately equal to 4.69. Thus, periodic solutions appear.Figure 1
The graph of the functionM0(τ). Two critical values of τ (i.e., τ* and τ**), for which stability switches can occur, appear.Using DDE23, a Matlab Solver for delay differential equation, we can compute the solutions of (1.4) for the above-mentioned values of the parameters, see Figures 2–4. Before the Hopf bifurcation occurs (0≤τ<τ*), solutions are asymptotically stable and converge to the positive equilibrium. When τ*≤τ≤τ**, periodic solutions appear (See Figure 3). When τ increases, longer-periods oscillations exist. For τ=τ**(=13.78), a stability switch occurs, that is, the positive equilibrium becomes asymptotically stable again and solutions converge to the positive equilibrium (Figure 4). In Figure 5, using the time delay τ as the bifurcation parameter, we can obtain the bifurcation diagram by Matlab package DDE-BIFTOOL.Time evolution of all the population for the model (1.4) with τ=3.3 and initial value (50,80,160). (a), (b) and (c) are S, I and y versus t, respectively, and (d) is the projected S-I-y phase plane.
(a)(b)(c)(d)Time evolution of all the population for the model (1.4) with τ=5.6 and initial value (50,80,160). (a), (b) and (c) are S, I and y versus t, respectively, and (d) is the projected S-I-y phase plane.
(a)(b)(c)(d)Time evolution of all the population for the model (1.4) with τ=15.2 and initial value (50,80,160). (a), (b) and (c) are S, I and y versus t, respectively, and (d) is the projected S-I-y phase plane.
(a)(b)(c)(d)Figure 5
Bifurcation diagram forτ.


## 5. Discussion

In this paper, we consider the ecoepidemiology model with time delay and Leslie-Gower schemes. The existence of a Hopf bifurcation at a positive steady-state is obtained through the study of an exponential polynomial characteristic equation with delay-dependent coefficients. That is to say, the positive steady-state of this model depends on the maturation periodτ through a factor e-μ1τ, where μ1 is the rate of disease-induced death. We can find stability switch occurs from stability to instability. Our computer simulation can easily confirm this. So the time delay plays a stabilizing role for our model.In addition, ifμ1=0, the steady state does not depend on time delay τ. From [28], we can see that the stability switch occurs toward instability if it is stable at τ=0, and further increase of time delay does not cause the occurrence of stability switch. This implies that both time delay and disease-induced death rate are necessary for the asymptotical stability of system (1.4).Finally, although we give the reason why we just think of such situation by citing that the rate of the predation on infected prey may be 31 times higher compared to that on susceptible prey [17], it is not usual in real world. In fact, if we deal with a huge population of the prey, then the susceptible population may be too large to ignore their role as prey of the predator. Thus, we can present the more realistic model as follows:
(5.1)dS(t)dt=rS(t)(1-S(t)K)-βe-μ1τS(t)I(t-τ)-μ1S(t),dI(t)dt=βe-μ1τS(t)I(t-τ)-cy(t)I(t)my(t)+I(t)-μ2I(t),dy(t)dt=δy(t)(1-hy(t)I(t)+σS(t)),
where μ1<μ2 and 0≤σ<1. Study on the dynamic behavior of system (5.1) will be obtained in the future.



---

*Source: 196204-2010-12-30.xml*