# Stability and Hopf Bifurcation in a Computer Virus Model with Multistate Antivirus

Authors: Tao Dong; Xiaofeng Liao; Huaqing Li
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

By considering that people may immunize their computers with countermeasures in susceptible state, exposed state and using anti-virus software may take a period of time, a computer virus model with time delay based on an SEIR model is proposed. We regard time delay as bifurcating parameter to study the dynamical behaviors which include local asymptotical stability and local Hopf bifurcation. By analyzing the associated characteristic equation, Hopf bifurcation occurs when time delay passes through a sequence of critical value. The linerized model and stability of the bifurcating periodic solutions are also derived by applying the normal form theory and the center manifold theorem. Finally, an illustrative example is also given to support the theoretical results.


---

## Body

## 1. Introduction

As globalization and development of communication networks have made computers more and more present in our daily life, the threat of computer viruses also becomes an increasingly important issue of concern. In 2003, a virus, called worm king, rapidly spread and attacked the global world, which results the network of the internet to be seriously congested and server to be paralyzed [1]. In 2010, the report of pestilence about computer virus in China revealed that more than 90% computers in China are infected computer virus.Computer viruses are small programs developed to damage the computer systems erasing data, stealing information. Their action throughout a network can be studied by using classical epidemiological models for disease propagation [2–6]. In [7–9], based on SIR classical epidemic model, Mark had proposed the dynamical models for the computer virus propagation, which provided estimations for temporal evolutions of infected nodes depending on network parameters [10–12]. In [13], Richard and Mark propose a modified propagation model named SEIR (susceptible-exposed-infected-recover) model to simulate virus propagation. In [14], on this basis of the SIR model, Yao et al. proposed a SIDQV model with time delay which add a quarantine state to clean the virus. However, both above models assume the viruses are cleaned in the infective state. In fact, in addition to clean viruses in state I, people may immunize their computers with countermeasures in state S and state E in the real world. Moreover there may be a time lag when the node uses antivirus software to clean the virus.In this paper, in order to overcome the above-mentioned limitation, we present a new computer virus model with time delay which is depending on the SEIR model [15];  time delay can be considered the period of the node uses antivirus software to clean the virus. This model provides an opportunity for us to study the behaviors of virus propagation in the presence of antivirus countermeasures, which are very important and desirable for understanding of the virus spread patterns, as well as for management and control of the spread. The remainder of this paper is organized as follows. In Section 2, the stability of trivial solutions and the existence of Hopf bifurcation are discussed. In Section 3, a formula for determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions will be given by using the normal form and center manifold theorem introduced by Hassard et al.  in [16]. In Section 4, numerical simulations aimed at justifying the theoretical analysis will be reported.


## 2. Mathematical Model Formulation

Our model is based on the traditional SEIR model [7–9, 15, 17]. The SEIR model has four states: susceptible, exposed (infected but not yet infectious), infectious, and recovered. Our assumptions on the dynamical model are as follows.(1)
In the real world, in addition cleaning viruses in state I, people may immunize their computers with countermeasures in state S and state E (after virus being cleaned), which may result in new state transition paths in comparison with SIR model:
S-R: using countermeasure of real-time immunization,E-R: using real-time immunization after virus codes cleaning.(2)
In state S, when people install the antivirus software on their computer, we assume that their computer can be immunized at a unit time.(3)
In state E, since the computer is infected by the virus, the antivirus software may use a period to search the document and clean the viruses.(4)
Denote the period of time of killing viruses when users find that their computers are infected by viruses.(5)
While the computer is installed the antivirus software, it will not be quarantine or replacement. On the basis of the above hypotheses (1)–(5), the dynamical model can be formulated by the following equations:(2.1)dS(t)dt=uN-βI(t)S(t)-(ρSR+μ)S(t),dE(t)dt=βI(t)S(t)-(α+μ)E(t)-ρERE(t-τ),dI(t)dt=αE(t)-(γ+μ)I(t),dR(t)dt=ρSRS(t)+ρERE(t-τ)+γI(t-τ)-μR(t),
where ρSR describes the impact of implementing real-time immunization, ρER describes the impact of cleaning the virus and immunizing the nodes, and μ describes the impact of quarantine or replacement. α is the transition rate from E to I, and γ is the recovery rate from I to R. τ is the time delay that the node usees antivirus software to clean the virus. β is the transition rate from S to E.


## 3. Local Stability of the Equilibrium and Existence of Hopf Bifurcation

We may see that the first three equations in (2.1) are independent of the fourth equation, and therefore, the fourth equation can be omitted without loss of generality. Hence, model (2.1) can be rewritten as(3.1)dS(t)dt=uN-βI(t)S(t)-(ρSR+μ)S(t),dE(t)dt=βI(t)S(t)-(α+μ)E(t)-ρERE(t-τ),dI(t)dt=αE(t)-(γ+μ)I(t).
For the convenience of description, we define the basic reproduction number of the infection as(3.2)R0=μNβα(ρSR+μ)(α+ρER+μ)(γ+μ).
Clearly, we have the following results with respect to the stable state of system (3.1). Here, the proof is omitted (see [17] for the details).Theorem 3.1.
IfR0<1, system (3.1) has only the disease-free equilibrium E0=(μN/(ρSR+μ),0,0) and is globally asymptotically stable. If R0>1, E0 becomes unstable and there exists a unique positive equilibrium Eve, where Eve=(μN/(ρSR+μ)R0,μN(R0-1)/R0(α+μ+ρER),αE*/(γ+μ)). Furthermore, for any τ>0, E0 is asymptotically stable if R0<1 and unstable if R0>1.To investigate the qualitative properties of the positive equilibriumE* with τ>0, it is necessary to make the following assumption:(H1)R0>1.Under hypothesis (H1), the Jacobian matrix of the system (3.1) about Eve is given by(3.3)J(Eve)=[-a10-a2a3-a4-a7e-λτa20a5-a6],
where a1=βI+ρSR+μ,a2=βS,a3=βI,a4=α+μ,a5=α,a6=γ+μ,a7=ρER.We can obtain the following characteristic equation:(3.4)λ3+b1λ2+b2λ+b3+e-λτ(b4λ2+b5λ+b6)=0,
where(3.5)b1=a1+a4+a6,b2=a1a6+a4(a1+a6)-a2a5,b3=a1a4a6-a1a2a5+a2a3a5,b4=a7,b5=a7(a1+a6),b6=a1a6a7.
If iω(ω>0) is a root of (3.4), then(3.6)-iω3-b1ω2+b2iω+b3+e-iωτ(-ω2b4+b5iω+b6)=0.
Separating the real and imaginary parts of (3.6), we have(3.7)b5ωsin⁡ωτ+(b6-b4ω2)cos⁡ωτ=b1ω2-b3,b5ωcos⁡ωτ-(b6-b4ω2)sin⁡ωτ=ω3-b2ω.
Adding up the squares of (3.7) yields (3.8)ω6+(b12-2b2-b42)ω4+(b22-2b1b3+2b4b6-b52)ω2+(b32-b62)=0.
Letting z=ω2,c1=b12-2b2-b42,c2=b22-2b1b3+2b4b6-b52,c3=b32-b62, then (3.8) becomes(3.9)z3+c1z2+c2z+c3=0.
Letting z*=(1/3)(-c1+c12-3c2),h(z*)=(z*)3+c1(z*)2+c2z*+c3, then we have the following results (see [18–22] for details) about the distributions of the positive roots of (3.9).Lemma 3.2 (see [18–22]).
(i) Ifc3<0, then (3.9) has at least one positive root.
(ii) Ifc3≥0 and c12-3c2≤0, then (3.9) has no positive root.
(iii) Ifc3≥0 and c12-3c2>0, then (3.9) has positive roots if and only if z*>0 and h(z*)≤0.Suppose (3.9) has positive roots; without loss of generality, we assume that it has three positive roots defined by ωk=zk, k=1,2,3. By (3.7), we have(3.10)cos⁡(ωkτ)=(b1ωk2-b3)(b6-b4ωk2)+b5ωk2(ωk2-b2)b52ωk2+(b6-b4ωk2)2.
Thus, denoting(3.11)τkj=1ωkarccos⁡(b1ωk2-b3)(b6-b4ωk2)+b5ωk2(ωk2-b2)b52ωk2+(b6-b4ωk2)2+2jπωk,
where k=1,2,3; j=0,1,…, then ±iω is a pair of purely imaginary roots of (3.4) with τkj. Define
(3.12)τ0=τk00=min⁡k=1,2,3{τk00},ω0=ωk0.
Note that when τ=0, (3.4) becomes(3.13)λ3+(b1+b4)λ2+(b2+b5)λ+(b3+b6)=0.
In addition, Routh-Hurwitz criterion [13] implies that, if the following condition holds, then all roots of (3.13) have negative real parts.(H2)(b1+b4)>0,(b1+b4)(b2+b5)-(b3+b6)>0.Till now, we can employ a result from Ruan and Wei [23] to analyze (3.4), which is, for the convenience of the reader, stated as follows.Lemma 3.3 (see [23]).
Consider the exponential polynomial(3.14)P(λ,e-λτ,…,e-λτm)=λn+p1(0)λn-1+⋯+pn-1(0)λ+pn(0)+[p1(1)λn-1+⋯+pn-1(1)λ+pn(1)]e-λτ+⋯+[p1(m)λn-1+⋯+pn-1(m)λ+pn(m)]e-λτm,
where τi≥0 (i=1,2,…,m) and pj(i)(j=1,2,…,m) are constants. As(τ1,τ2,…,τm) vary, the sum of the order of the zeros of P(λ,e-λτ,…,e-λτm) on the open right half plane can change only if a zero appears on or crosses the imaginary axis.Using Lemmas3.2 and 3.3 we can easily obtain the following results on the distribution of roots of the transcendental (3.4).Lemma 3.4.
(2.1) If c3>0 and c12-3c2≤0, then all roots with positive real parts of (3.4) have the same sum as those of the polynomial (3.13) for all τ≥0.
(3.1) If either c3<0 or c3≥0 and c12-3c2>0, z*>0,h(z*)≤0, then all roots with positive real parts of (3.4) have the same sum as those of the polynomial (3.13) for τ∈[0,τ0).Lemma 3.5.
If3wk4+c1wk1+c2≠0, then the following transversality condition holds:
(3.15)sgn⁡{Re{(dλdτ)-1}}≠0whenτ=τ0.Proof.
Differentiating (3.4) with respect to τ yields
(3.16)[3λ2+2b1λ+b2+(2b4λ+b5-τ(b4λ2+b5λ+b6)e-λτ)]dλdτ=λ(b4λ2+b5λ+b6)e-λτ.
For the sake of simplicity, denoting ω0 and τ0 by ω,τ respectively, then
(3.17)(dλdτ)-1=3λ2+2b1λ+b2λ(b4λ2+b5λ+b6)e-λτ+2b4λ+b5λ(b4λ2+b5λ+b6)-τλ=-2λ3+b1λ2-b3λ2(λ3+b1λ2+b2λ+b3)+b4λ2-b6λ(b4λ2+b5λ+b6)-τλ=-2iω3-b1ω2-b3ω2(b3-b1ω2-i(ω3-b2ω))+b4ω2+b6ω2(b6-b3ω2+b4iω)-τiω.
Then we get
(3.18)Re{(dλdτ)-1}=-1ω2[b33-2ω6-(b12-2b2)ω4(b3-b1ω2)2+(ω3-b2ω)2+b62-b42ω4(b6-b4ω2)2+b52ω2]=2ω6+c1ω4-c3ω2((b6-b4ω2)2+b52ω2)=3ω4+c1ω2+c2((b6-b4ω2)2+b52ω2).
Then, if 3ω4+c1ω2+c2≠0, we have sgn{Re{(dλ/dτ)-1}}≠0, we complete proof.Thus from Lemmas3.2, 3.3, 3.4, and 3.5, and we have the following.Theorem 3.6.
Suppose that (H1) and (H2) hold, then the following results hold.(1)
The positive equilibrium of (3.1) is asymptotically stable, if c3>0 and c12-3c2≤0;(2)
if eitherc3<0 or c3≥0 and c12-3c2>0, z*>0,h(z*)≤0, system (3.1) is asymptotically stable for τ∈[0,τ0) and system (3.1) undergoes a Hopf bifurcation at the origin when τ=τ0.


## 4. Direction of the Hopf Bifurcation

In this section, we derive explicit formulae for computing the direction of the Hopf bifurcation and the stability of bifurcation periodic solution at critical valuesτ0 by using the normal form theory and center manifold reduction.Lettingx1=S-S*,x2=E-E*,x3=I-I*,x̅i(t)=xi(τt),τ=τ0+μ, and dropping the bars for simplification of notation, system (3.1) is transformed into an FDE as(4.1)ẋ(t)=Lμ(xt)+f(μ,xt),
with(4.2)Lμφ=(τ0+μ)[B1φ(0)+B2φ(-1)],
where(4.3)B1=[-a10-a2a3-a4a20a5-a6],B2=[0000-a7e-λτ0000],f(μ,φ)=(τ0+μ)(-βφ1(0)φ2(0)βφ1(0)φ2(0)0).
Using the Riesz representation theorem, there exists a function η(θ,μ) of bounded variation for θ∈[-1,0], such that(4.4)Lμφ=∫-10dη(θ,μ)φ(θ)φ∈C.
In fact, we can choose(4.5)η(θ,μ)=(τ0+μ)[B1δ(θ)+B2δ(θ+1)],
where δ(θ) is Dirac delta function.In the next, forφ∈[-1,0], we define(4.6)A(μ)φ={dφdθ,θ∈[-1,0),∫-10dη(θ,μ)φ(θ),θ=0,(4.7)R(μ)φ={0,θ∈[-1,0],f(μ,φ),θ=0.
Then system (4.2) can be rewritten as(4.8)ẋ(t)=A(μ)xt+R(μ)xt,
where xt(θ)=x(t+θ).The adjoint operatorA* of A is defined by(4.9)A*(μ)ψ={-dψ(s)dθ,s∈(0,1],∫-10dηT(t,0)ψ(-t),s=0,
where ηT is the transpose of the matrix η.Forφ∈C1[-1,0] and ψ∈C1[0,1], we define(4.10)〈ψ,φ〉=ψ̅(0)⋅φ(0)-∫θ=-10∫ξ=0θψ̅(ξ-θ)dη(θ)φ(ξ)dξ,
where η(θ)=η(θ,0). We know that ±iτ0ω0 is an eigenvalue of A(0), so ±iτ0ω0 is also an eigenvalue of A*(0). We can get(4.11)q(θ)=(1q1q2)eiτ0ω0θ,-1<θ≤0.
From the above discussion, it is easy to know that(4.12)Aq(0)=iτ0ω0q(0).
Hence we obtain(4.13)q1=iω0q2a5,q2=-iω0+a1a2.
Suppose that the eigenvector q* of A* is(4.14)q*(s)=(1q1*q2*)eiτ0ω0s,
Then the following relationship is obtained:(4.15)A*q(0)=-iτ0ω0q*(0).
Hence we obtain(4.16)q1*=a1-iω0a3,q2*=a4+a7eiω0τ0a5q1*.
Let(4.17)〈q*,q〉=1.
One can obtain(4.18)〈q*,q〉=q̅*(0)⋅q(0)-∫θ=-10∫ξ=0θq̅*T(ξ-θ)dη(θ)φ(ξ)dξ=1ρ(1+q1q̅1*+q2q̅2*)-∫θ=-10∫ξ=0θτ01ρ(1q̅1*q̅2*)[(-a10-a2a3-a4a20a5-a6)δ(θ)+(0000-a70000)δ(θ+1)](1q1q2)eiτ0ω0θdξdθ=1ρ(1+q1q̅1*+q2q̅2*)-1ρτ0e-iω0τ0a7q1q̅1*.Hence we obtain(4.19)ρ=(1+q1q̅1*+q2q̅2*)-τ0e-iω0τ0a7q1q̅1*.In the remainder of this section, by using the same notations as in Hassard et al. [16], we first compute the coordinates for describing the center manifold C0 at μ=0. Leting  xt be the solution of (4.1) with μ=0, we define(4.20)z(t)=〈q*,xt〉,W(t,θ)=xt-2Re{z(t)q(θ)}.
On the center manifold C0 we have(4.21)W(t,θ)=W(z,z̅,t),
where(4.22)W(z,z̅,t)=W20(θ)z22+W11(θ)zz̅2+W02(θ)z̅22+⋯.In fact,z and z̅ are local coordinate for C0 in the direction of q and q*. Note that, if xt is, we will deal with real solutions only. Since μ=0(4.23)ż(t)=〈q*,ẋt〉=〈q*,A(μ)xt+R(μ)xt〉=〈q*,Axt〉+〈q*,Rxt〉=iτ0w0z+q̅*(0)⋅f(0,W(t,0)+2Re[z(t)q(0)]).
Rewrite (4.23) as(4.24)ż(t)=iτ0ω0z+g(z,z̅),
where(4.25)g(z,z̅)=g20z22+g11zz̅+g02z̅22+g21z2z̅2+⋯.
From (4.1) and (4.24), we have(4.26)Ẇ=ẋt-żq-z̅̇q̅̇={AW-2Re[q̅*(0)⋅f(z,z̅)q(θ)],θ∈[-2τ,0),AW-2Re[q̅*(0)⋅f(z,z̅)q(θ)]+f0(z,z̅),θ=0.
Let(4.27)Ẇ=AW+H(z,z̅,θ),
where(4.28)H(z,z̅,θ)=H20(θ)z22+H11(θ)zz̅2+H02(θ)z̅22+⋯.
Expanding the above series and comparing the corresponding coefficients, we obtain(4.29)(A-2iw0)W20(θ)=-H20(θ),AW11(θ)=-H11(θ),(A+2iw0)W02(θ)=-H02(θ).
Since xt=x(t+θ)=W(z,z̅,θ)+zq+z̅·q̅, we have(4.30)xt=(W(1)(z,z̅,θ)W(2)(z,z̅,θ)W(3)(z,z̅,θ))+z(1q1q2)eiω0θ+z̅(1q̅1q̅2)e-iω0θ.
Thus, we can obtain(4.31)φ1(0)=z+z̅+W20(1)(0)z22+W11(1)(0)zz̅+W02(1)(0)z̅22,φ2(0)=zq1+z̅q̅1*+W20(1)(0)z22+W11(1)(0)zz̅+W02(1)(0)z̅22.
So(4.32)φ1(0)φ2(0)=q1z2+q̅1z̅2+(q1+q2)zz̅+(W11(2)+12W20(2)+W11(1)q1+12W20(1)q̅)z2z̅.
It follows from (4.24) and (4.25) that(4.33)f(φ,μ)=(K11z2+K12zz̅+K13z̅2+K14z2z̅K21z2+K22zz̅+K23z̅2+K24z2z̅0),
where(4.34)K11=-βq1,K12=-βq̅1,K13=-β(q1+q̅1),K14=-β(W11(2)+12W20(2)+W11(1)q1+12W20(1)q̅),K21=βq1,K22=βq̅1,K23=β(q1+q̅1),K24=β(W11(2)+12W20(2)+W11(1)q1+12W20(1)q̅).
Since q*(0)=(1/ρ)(1,q1*¯,q2*¯)T, we have(4.35)g(z,z̅)=1ρ(1,q1*¯,q2*¯)(K11z2+K12zz̅+K13z̅2+K14z2z̅K21z2+K22zz̅+K23z̅2+K24z2z̅0).
Comparing the coefficients of the above equation with those in (4.27), we have(4.36)g20=1ρ(K11+K21q1*¯),g11=1ρ(K12+K22q1*¯),g02=1ρ(K13+K23q1*¯),g21=1ρ(K14+K24q1*¯).
In what follows, we focus on the computation of W20(θ) and W11(θ). For the expression of g21, we have(4.37)H(z,z̅,θ)=-2Re[q̅*(0)⋅f(z,z̅)q(θ)]=-(g20z22+g11zz̅+g02z̅22+⋯)q(θ)-(g̅20z̅22+g̅11zz̅+g̅02z22+⋯)q(θ).
Comparing the coefficients of the above equation, we can obtain that(4.38)H20(θ)=-g20q(θ)-g̅02q̅(θ),θ∈[-1,0),(4.39)H11(θ)=-g11q(θ)-g̅11q̅(θ),θ∈[-1,0).
Substituting (4.39) into (4.27) and (4.38) into (4.27), respectively, we get(4.40)Ẇ20(θ)=2iτ0ω0W20(θ)+g20q(θ)+g̅20q̅(θ),Ẇ11(θ)=+g11q(θ)+g̅11q̅(θ).
So(4.41)W20(θ)=ig20τ0ω0q(0)eiτ0ω0θ-g̅023iτ0ω0q̅(0)e-iτ0ω0θ+E1e2iτ0ω0θ,W11(θ)=g11iτ0ω0q(0)eiτ0ω0θ-g̅11iτ0ω0q̅(0)e-iτ0ω0θ+E2.
In the sequel, we will determine E1 and E2. Form the definition of A in (4.8), we have(4.42)∫-10dη(θ)W20(θ)=2iτ0ω0W20(0)-H20(0),(4.43)∫-10dη(θ)W11(θ)=-H11(0).
From (4.6) and (4.38)-(4.39), we have(4.44)H20(θ)=-g20q(θ)-g̅02q̅(θ)+(K11,K21,0)T,(4.45)H11(θ)=-g11q(θ)-g̅11q̅(θ)+(K12,K22,0)T.
Substituting (4.41) and (4.44) into (4.42) and noticing that(4.46)(iω0I-∫-10eiω0θdη(θ))q(0)=0,(-iω0I-∫-10e-iω0θdη(θ))q̅(0)=0,
we can obtain(4.47)(2iω0I-∫-10e2iτ0ω0θdη(θ))E1=(K11K210)T,
which leads to(4.48)(2iω0+a10a2-a32iω0+a4+a7e-2iω0τ0-a20-a52iω0+a6)(E1(1)E1(2)E1(3))=(K11K210),(a10a2-a3a4+a7-a20-a5a6)(E2(1)E2(2)E2(3))=(K12K220).
It follows that(4.49)E1(1)=K11-a2E1(3)2iω0+a1,E1(2)=2iω0+a6a5E1(3),E1(3)=K22-K11/(2iω0+a1)a2a3/(2iω0+a1)+(2iω0+a4+a7e-2iω0τ0)((2iω0+a6)/a5)-2iω0-a6,E2(1)=K12-a2E2(3)a1,E2(2)=a6a5E2(3),E2(3)=a1a5K22-a5K12a2a3a5+(a4+a7)a1a6-a1a5a6.
Based on the above analysis, we can see each gij in (4.37) is determined by parameters and delays in (3.1). Thus, we can compute the following quantities:(4.50)μ2=-ReC1(0)Reλ′(τ0),T2=-Im⁡C1(0)+μ2Im⁡λ′(0)ω0,β2=-2ReC1(0).Theorem 4.1.
In (4.50), the following results hold.(1)
The sign ofμ2 determines the directions of the Hopf bifurcation: if μ2>0(μ2<0) then the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for τ>τ0(τ<τ0).(2)
The sign ofβ2 determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if β2<0(β2>0).(3)
The sign ofT2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2>0(T2<0).


## 5. Numerical Examples

In this section, some numerical results of system (3.1) are presented to justify the Previous theorem above. As an example, considering the following parameters: μ=0.01,N=10000,γ=0.08,α=0.1,β=0.01,ρSR=0.2,ρER=0.2, then R0=1.706,c3=-3.6284e-5, and E*=(279,133.6,148.4). According to the Lemma 3.2, (3.9) has one positive real root ω=0.1194. Correspondingly, by (3.13), we obtain τ0=14.05. First, we choose τ=13<τ0, the corresponding wave form and phase plots are shown in Figure 1; it is easy to see from Figure 1 that system (3.1) is asymptotically stable. Finally, we choose τ=14.15>τ0 the corresponding wave form and phase plots are shown in Figure 2; it is easy to see that Figure 2 undergoes a Hopf bifurcation.τ=13<τ0. The positive equilibrium E* of system (3.1) is asymptotically stable.
(a)(b)(c)(d)τ=14.15>τ0. The bifurcation periodic solution for system (3.1) is stable.
(a)(b)(c)(d)


## 6. Conclusions

In this paper, considering that in addition to cleaning viruses in state I, people may immunize their computers with countermeasures in stateS and state E, and since using antivirus software will take a period of time, we have constructed a computer virus model with time delay depending on the SEIR model. The theoretical analyses for the computer virus models are given. Furthermore, we have proved that when time cross through the critical value, the system exist a Hopf bifurcation. Finally, simulation clarifies our results.



---

*Source: 841987-2012-04-17.xml*