# Dynamics of SCIR Modeling for COVID-19 with Immigration

Authors: Haiyin Li; Yan Wu
Journal: Complexity
 (2022)

---

## Abstract

In this study, for COVID-19, we divide people into four categories: susceptibleSt, closely contacted Ct, infective It, and removed Rt according to the current epidemic situation and then investigate two models: the SCIR models with immigration (Model (2) and without immigration (Model 1). For the former, Model 1, we obtain the condition for global stability of its disease-free equilibrium. For the latter, Model 2, we establish the local asymptotic stability of its endemic equilibrium by constructing Lyapunov function. Afterwards, by the bifurcation theory, we qualitatively analyze the properties of its Hopf bifurcations of the latter. Finally, numerical simulations are given to illustrate the obtained results of two models. The results imply the importance of finding closely contacted and overseas imports on epidemic control. It indicates that not only the incubation delay τ is crucial for the containment of the COVID-19 but also the scientific and rigorous containment measures are the key factors of the success of the containment.


---

## Body

## 1. Introduction

The 2019 novel coronavirus disease (COVID-19) is running rampantly in the world, which has caused a great concern on the global public health. Research on infectious diseases has never stopped [1–10]. In [11], the authors investigated global stability for the SEIR model in epidemiology. The dynamical behavior of an epidemic model with nonlinear incidence was studied in Ruan and Wang [12]. Peng [13] considered asymptotic curve of the positive steady state of the SIS epidemic reaction diffusion system. From 2020 to 2021, to defeat the epidemic, scientists in different fields cooperated to investigate COVID-19 from different points of view. Sex ratio, media reports, super-spreaders, cytokine storm, mask wearing, healthcare, vaccination, age structure, impact on energy demand and consumption, etc., are covered [14–26].In order to prevent the spread of the epidemic, the suspected cases and occasional local cases should be checked, including nasopharyngeal and oropharyngeal swabs, CT scans, serum testing, bronchoalveolar lavage, and other different methods of inspection [27]. In [28], their model is based on the usual SIR model, dividing the total proportion of infected individuals into two parts: patients who have not yet been detected and patients who have been detected by tests. In addition, because a large number of people travel from one country to another, the disease is mainly spread between different countries through air travel. For countries with better epidemic control, imported cases have become the main source of cases. Therefore, it is important to consider foreign inputs.Based on the global pandemic COVID-19, the epidemic lasts for such a long time; personnel exchanges between countries are inevitable. In this study, we discuss the impact of immigration between countries on the prevention and control of epidemic. The total population of each group is divided into four compartments, that is, susceptible, closely contacted, infective, and removed compartments. We will investigate two models: the SCIR models with immigration (Model (2) and without immigration (Model 1), respectively.The arrangement of this study is as follows. In Section2, we establish the global stability of disease-free equilibrium for Model 1. In Section 3, we obtain the local asymptotic stability of endemic equilibrium and analyze the Hopf bifurcations for Model 2. In Section 4, this study is concluded.


## 2. Model 1

### 2.1. Construction of Model 1

Take one country as a whole, not taking into account the number of imported people, the birth rate, or the normal death rate in the countries.The model (Model 1) reads(1)S′t=aCt−bStItN,C′t=bStItN−aCt−βCt−τ,I′t=βCt−τ−rIt,R′t=rIt,where St means the number of susceptible individuals, Ct means the closely contacted individuals, that is the exposed, It means the infected at time t, Rt is removed from the infected system including recovered and dead, and τ indicates the incubation delay, i.e., the duration between the time when the people who is exposed to the infected cases, but not yet infectious, and later the time when they become infectious. And a means transition rate of closely contacted individuals to susceptible by testing and quarantining, b means contact rate of transmission per contact from infected class, β means transition rate of closely contacted individuals to the infected case, and r means removal rate including mortality of infectious disease and recovery rate of infected individuals. Furthermore, N is nearly a constant, the total population of a country or region.SinceRt does not intervene in the dynamics of St,Ct, and It, the following system can be separated from system (1):(2)S′t=aCt−bStItN,C′t=bStItN−aCt−βCt−τ,I′t=βCt−τ−rIt.The initial conditions of the delayed system (2) are(3)S0θ=φ1θ,C0θ=φ2θ,I0θ=φ3θ−τ≤θ≤0;Sθ>0,Cθ≥0,Iθ≥0,where φ=φ1,φ2,φ3T∈C−τ,0,R+3 is the initial function, C−τ,0,R+3 is the Banach space of continuous functions mapping the interval −τ,0 into R+3, R+3=S,C,I:S≥0,C≥0,I≥0 and φ=maxφθ:θ∈−τ,0, and φθ is any norm in R+3. Usually, we use regular symbols Ctθ=Ct+θ for θ∈−τ,0.


### 2.2. Global Asymptotical Stability for Model 1

System (2) does not have the endemic equilibrium and has only one equilibrium, i.e., the disease-free equilibrium E0S0,0,0, where S0 is the number of the initial susceptible individuals.Theorem 1.
Ifb<r/2, then E0S0,0,0 of system (2) has global asymptotical stability.Proof .
AssumeVt=St−S0+Ct+It. If St>S0, Vt=St−S0+Ct+It, we have(4)V′t=S′t+C′t+I′t=−rI<0.
IfSt<S0,Vt=S0−St+Ct+It, we have(5)V′t=−2aCt+2bStN−rIt.
Because ofSt/N≤1, if 2b<r, we can get V′t≤−2aCt<0.
IfSt=S0,Vt=Ct+It, we have(6)V′t=−aCt+bStN−rIt.
Because ofSt/N≤1, if b<r, we can get V′t≤−aCt<0.
Therefore, whenb<r/2, then the unique equilibrium, i.e., the disease-free equilibrium.
E0 has global asymptotical stability.Example 1.
Leta=0.1,b=0.01,β=0.058,r=0.025,andN=1.0×107; then, system (2) reads(7)S′t=0.1Ct−0.01StIt1.0×107,C′t=0.01StIt1.0×107−0.1Ct−0.058Ct−τ,I′t=0.058Ct−τ−0.025It.For system (7), apparently b−r/2<0, by Theorem 1, the equilibrium E0S0,0,0 has global asymptotical stability. We have Figure 1 by simulations. Note that, in Figure 1, the time histories of St,Ct,andIt in system (7), for τ=4, with initial functions S−τ,C−τ,I−τ=107,30,5 and S0,C0,I0=107,2000,100 and for t:0−4000a,c,e and t:7000−8000b,d,f are presented, respectively, and all converge to 107,0,0 as t tends to +∞.Figure 1
The time histories ofSt,Ct,andIt in system (7), for τ=4, with initial functions S−τ,C−τ,I−τ=107,30,5 and S0,C0,I0=107,2000,100, and t:0−4000a,c,e and t:7000−8000b,d,f, respectively.
(a)(b)(c)(d)(e)(f)


## 2.1. Construction of Model 1

Take one country as a whole, not taking into account the number of imported people, the birth rate, or the normal death rate in the countries.The model (Model 1) reads(1)S′t=aCt−bStItN,C′t=bStItN−aCt−βCt−τ,I′t=βCt−τ−rIt,R′t=rIt,where St means the number of susceptible individuals, Ct means the closely contacted individuals, that is the exposed, It means the infected at time t, Rt is removed from the infected system including recovered and dead, and τ indicates the incubation delay, i.e., the duration between the time when the people who is exposed to the infected cases, but not yet infectious, and later the time when they become infectious. And a means transition rate of closely contacted individuals to susceptible by testing and quarantining, b means contact rate of transmission per contact from infected class, β means transition rate of closely contacted individuals to the infected case, and r means removal rate including mortality of infectious disease and recovery rate of infected individuals. Furthermore, N is nearly a constant, the total population of a country or region.SinceRt does not intervene in the dynamics of St,Ct, and It, the following system can be separated from system (1):(2)S′t=aCt−bStItN,C′t=bStItN−aCt−βCt−τ,I′t=βCt−τ−rIt.The initial conditions of the delayed system (2) are(3)S0θ=φ1θ,C0θ=φ2θ,I0θ=φ3θ−τ≤θ≤0;Sθ>0,Cθ≥0,Iθ≥0,where φ=φ1,φ2,φ3T∈C−τ,0,R+3 is the initial function, C−τ,0,R+3 is the Banach space of continuous functions mapping the interval −τ,0 into R+3, R+3=S,C,I:S≥0,C≥0,I≥0 and φ=maxφθ:θ∈−τ,0, and φθ is any norm in R+3. Usually, we use regular symbols Ctθ=Ct+θ for θ∈−τ,0.


## 2.2. Global Asymptotical Stability for Model 1

System (2) does not have the endemic equilibrium and has only one equilibrium, i.e., the disease-free equilibrium E0S0,0,0, where S0 is the number of the initial susceptible individuals.Theorem 1.
Ifb<r/2, then E0S0,0,0 of system (2) has global asymptotical stability.Proof .
AssumeVt=St−S0+Ct+It. If St>S0, Vt=St−S0+Ct+It, we have(4)V′t=S′t+C′t+I′t=−rI<0.
IfSt<S0,Vt=S0−St+Ct+It, we have(5)V′t=−2aCt+2bStN−rIt.
Because ofSt/N≤1, if 2b<r, we can get V′t≤−2aCt<0.
IfSt=S0,Vt=Ct+It, we have(6)V′t=−aCt+bStN−rIt.
Because ofSt/N≤1, if b<r, we can get V′t≤−aCt<0.
Therefore, whenb<r/2, then the unique equilibrium, i.e., the disease-free equilibrium.
E0 has global asymptotical stability.Example 1.
Leta=0.1,b=0.01,β=0.058,r=0.025,andN=1.0×107; then, system (2) reads(7)S′t=0.1Ct−0.01StIt1.0×107,C′t=0.01StIt1.0×107−0.1Ct−0.058Ct−τ,I′t=0.058Ct−τ−0.025It.For system (7), apparently b−r/2<0, by Theorem 1, the equilibrium E0S0,0,0 has global asymptotical stability. We have Figure 1 by simulations. Note that, in Figure 1, the time histories of St,Ct,andIt in system (7), for τ=4, with initial functions S−τ,C−τ,I−τ=107,30,5 and S0,C0,I0=107,2000,100 and for t:0−4000a,c,e and t:7000−8000b,d,f are presented, respectively, and all converge to 107,0,0 as t tends to +∞.Figure 1
The time histories ofSt,Ct,andIt in system (7), for τ=4, with initial functions S−τ,C−τ,I−τ=107,30,5 and S0,C0,I0=107,2000,100, and t:0−4000a,c,e and t:7000−8000b,d,f, respectively.
(a)(b)(c)(d)(e)(f)


## 3. Model 2

### 3.1. Construction of Model 2

With the joint efforts of the whole country, the epidemic in some countries has been brought under control. However, due to the immigration of imported cases and the virus on imported cold-chain food, the epidemic has occurred many times. For the convenience of discussion, we will collectively refer to imported personnel and imported cold-chain food as imported from abroad.The corresponding model (Model 2) is(8)S′t=Λ1+aCt−bStItN−dSSt,C′t=Λ2+bStItN−aCt−βCt−τ−dCCt,I′t=Λ3+βCt−τ−rIt−dIIt,R′t=rIt−dRRt,where Λ1,Λ2,andΛ3 present the immigration rates of susceptible class, closely contacted class, and infected class from the foreign countries, respectively. dS,dC,dI,anddR present the death rates of susceptible individuals, closely contacted individuals, infected, and removed, respectively. Moreover, N=S+C+I+R is a variable, which keeps increasing with the input of Λ1,Λ2, and Λ3 and holds decreasing with output of the population which we do not consider. Therefore, combining the two cases into consideration, we can assume that the population is nearly a constant, the total population of a country or region.System (8) has no disease-free equilibrium 0,0,0,0 or other boundary equilibrium; there exists only the endemic equilibrium E∗S∗,C∗,I∗,R∗. From (8), it follows that(9)C∗=Λ1+Λ2−dSS∗β+dC,I∗=Λ3r+dI+βΛ1+Λ2−dSS∗r+dIβ+dC,R∗=rI∗dR,where S∗ satisfies(10)m1S∗2−m2S∗+m3=0,m1=bβdS>0,m2=bβΛ1+Λ2+bΛ3β+dC+NdSr+dIβ+dC+a>0,m3=Nr+dI.Λ1β+dC+aΛ1+Λ2>0. We analyze the root of (10), by calculation:(11)Δ=m22−4m1m3=aNdSr+dI−bβΛ1+Λ22+NdSr+dIβ+dC2+bdCΛ32+2b2βdCΛ3Λ1+Λ2+Λ3+2abβNdSΛ3r+dI+2abNdSdCΛ3r+dI+2NdSr+dIβ+dCa+bβΛ2+Λ3+bdCΛ3−bβΛ1.WhenΔ>0, equation (10) has two positive roots S1 and S2. Furthermore, if Λ1+Λ2−dSSi>0,i=1,2, then system (8) has two positive equilibria E1∗S1∗,C1∗,I1∗,R1∗ and E2∗S2∗,C2∗,I2∗,R2∗.WhenΔ>0, equation (10) has two positive roots S1 and S2. However, if Λ1+Λ2−dSS1<0, or Λ1+Λ2−dSS2<0, then system (8) has the unique equilibrium, i.e., the positive equilibrium E2∗S2∗,C2∗,I2∗,R2∗ or E1∗S1∗,C1∗,I1∗,R1∗.WhenΔ=0, equation (10) has one positive root S∗. Furthermore, if Λ1+Λ2−dSS∗>0, then system (8) has the unique equilibrium, i.e., the positive equilibrium E∗S∗,C∗,I∗,R∗.


### 3.2. Locally Asymptotically Stability for Model 2

The Jacobian matrix [29] on E∗ of sysem (8) is(12)JE∗=−bI∗N−dSa−bS∗N0bI∗N−a−dC−βe−λτbS∗N00βe−λτ−r−dI000r−dR.Hence, we can get the characteristic equation as(13)λ+dRλ+bI∗N+dSλ+a+dC+βe−λτλ+r+dR+βe−λτbI∗NbS∗N−βe−λτbS∗Nλ+bI∗N+dS−abI∗Nλ+r+dR=0.So, we have the characteristic rootλ=−dR, and the other characteristic roots satisfying(14)λ3+a1λ2e−λτ+a2λ2+a3λe−λτ+a4λ+a5e−λτ+a6=0,with(15)a1=β,a2=a+dC+r+dI+bI∗N+dS,a3=βr+dI+bI∗N+dS−bS∗N,a4=bI∗NdC+r+dI+dSa+dC+r+dI+r+dIa+dC,a5=βr+dIbI∗N+dS−dSbS∗N,a6=dSr+dIa+dC+dCbI∗Nr+dI.By (14), when τ=0, we have(16)λ3+a1+a2λ2+a3+a4λ+a5+a6=0,with(17)a1+a2=β+a+dC+r+dI+bI∗N+dS>0,a3+a4=βr+dI+bI∗N+dS−bS∗N+bI∗NdC+r+dI+dSa+dC+r+dI+r+dIa+dC,a5+a6=βr+dIbI∗N+dS−dSbS∗N+dSr+dIa+dC+dCbI∗Nr+dI.Assume(18)Δ1=a1+a2,Δ2=a1+a21a5+a6a3+a4=a1+a2a3+a4−a5+a6,Δ3=a1+a210a5+a6a3+a4a1+a200a5+a6=a5+a6Δ2.Using Routh–Hurwotz criterion, when(19)Δ1>0,Δ2>0,Δ3>0.Then, all roots of (16) have negative real parts. So, we can obtain the following result.Theorem 2.
If condition (19) holds, all roots of (16) have negative real parts. The endemic equilibrium E∗ of system (8) is asymptotically stable for τ=0 if condition (19) holds.Next, we suppose condition (19) is established. In order to determine if the real parts of some roots of (9) reaches to zero and ultimately become positive as τ varies continuously, let λ=iωω≠0 be the root; we have(20)a5−a1ω2cosωτ+a3ωsinωτ=a2ω2−a6,a3ωcosωτ+a1ω2−a5sinωτ=ω3−a4ω.Thus,a5−a1ω22+a3ω2=a2ω2−a62+ω2ω2−a42, in other words,(21)ω6+a22−a12−2a4ω4+a42−a32+2a1a5−2a2a6ω2+a62−a52=0.Assumez=ω2, equation (21) can be converted to(22)z3+pz2+qz+g=0,with p=a22−a12−2a4,q=a42−a32+2a1a5−2a2a6,andg=a62−a52.Assumehz=z3+pz2+qz+gandΓ=p2−3q. We have the following result.Lemma 1 (see [30]).
(i) Wheng<0, equation (22) has at least one positive root.(i)
Wheng≥0 and Γ=p2−3q<0, equation (22) has no positive root.\(ii)
Wheng≥0 and Γ=p2−3q≥0, equation (22) has positive root if and only if z1=−p+Γ/3>0 and hz1≤0Therefore, the following theorem is obtained by Lemma1.Theorem 3.
When(23)g<0.Equation (22) has at least one positive root.Moreover, we suppose condition (23) is established. Equation (14) has at least one positive root, i.e., equation (22) has at least one positive root. And equation (22) is a cubic, so it may have one, two, or three positive roots. Without loss of generality, we suppose that it has three positive roots, expressed as z1, z2, and z3 Then, equation (21) has three positive roots, that is, ω1=z1,ω2=z2, and ω3=z3.Assume(24)τkj=1ωkθk+2jπ,k=1,2,3,j=0,1,2,….Withθk∈0,2π, meeting(25)sinθk=ωka3a2ωk2−a6+ωk2−a4a1ωk2−a5a3ωk2+a1ωk2−a52,cosθk=a2ωk2−a6a5−a1ωk2+a3ωk2ωk2−a4a3ωk2+a1ωk2−a52,so ±iωk are a pair of purely imaginary roots of equation (14) with τ=τkj,k=1,2,3,j=0,1,2,…; obviously,(26)limj⟶∞τkj=∞,k=1,2,3.So, we can assume(27)τ0=τk00=min1≤k≤3τk0,ω0=ωk0.When condition (23) is established, owing to Theorem 2, equation (14) always has positive root and has roots having zero real part for τ=τkj, and τ0 is the minimum of τ. Therefore, for all τ∈0,τ0, the roots of equation (14) have negative real parts. We can have the following theorem.Theorem 4.
Whenτ∈0,τ0, all roots of equation (14) have negative real parts. That is, when τ∈0,τ0, the endemic equilibrium E∗of system (8) is locally asymptotically stable.Assumeλτ=ατ+iωτ is the root of equation (14) satisfying ατ0=0andωτ0=ω0. We suppose h′z0≠0 in order to ensure that ±iω0 are simple pure imaginary roots of equation (14) for τ=τ0 and λτ satisfies the transversality condition.The two sides of Equation (14) are differentiated with respect toτ; we can obtain(28)dλτdτ=a1λ2+a3λ+a5λe−λτ3λ2+2a2λ+a4−e−λτa1λ2τ+a3λτ−2a1λ+a5τ−a3=Aτ+BτiDτ+Eτi,with(29)Aτ=e−ατa1α2−ω2+a3α+a5αcosωτ+ωsinωτ−e−ατ2a1αω+a3ωωcosωτ−αsinωτ,Bτ=e−ατa1α2−ω2+a3α+a5ωcosωτ−αsinωτ+e−ατ2a1αω+a3ωαcosωτ+ωsinωτ,Dτ=3α2−ω2+2a2α−a4−e−ατcosωτa1τα2−ω2+a3ατ−2a1α+a5τ−a3+e−ατsinωτ2a1ταω+a3τω−2a1ω,Eτ=6αω+2a2ω−e−ατcosωτ2a1ταω+a3τω−2a1ω−e−ατsinωτa1τα2−ω2+a3ατ−2a1α+a5τ−a3.Puttingτ=τ0 into the above equation, separating the real and imaginary parts, and by Equation (20) and the definition of hz, we obtain(30)dReλτdτ|τ=τ0=Aτ0Dτ0+Bτ0Eτ0D2τ0+E2τ0=ω023ω04+2a22−a12−2a4ω02+a42−a32+2a1a5D2τ0+E2τ0=ω02h′ω02D2τ0+E2τ0.WhendReλτ/dτ<0 for τ<τ0 and close to τ0, we obtain equation (14) has a root λτ=αλ+iωλ with αλ. This contradicts with Theorem 4, so we have the following conclusion.Theorem 5.
Whenτ=τ0, a pair of simple purely imaginary roots of equation (14) are ±iω0 and all other roots have negative real parts. In addition, dReλτ0/dτ>0.Hence,E∗ of system (8) is asymptotically stable (unstable) when τ<τ0τ>τ0, respectively, and τ0 satisfying Equation (27).


### 3.3. Local Asymptotic Stability of Model 2

AssumeSt=S∗+S˜Ct=C∗+C˜It=I∗+I˜Rt=R∗+R˜; then, system (8) is linearized as(31)S′t=−bI∗N+dSS˜t+aC˜t−bS∗NI˜t,C˜′t=bI∗NS˜t−a+dCC˜t−βC˜t−τ+bS∗NI˜t,I˜′t=βC˜t−τ−r+dII˜t,R˜′t=rI˜t−dRR˜t.For the simplicity of the symbol, we omit the tildes; system (31) can be rewritten as(32)S′t=−bI∗N+dSSt+aCt−bS∗NIt,C′t=bI∗NSt−a+dCCt−βCt−τ+bS∗NIt,I′t=βCt−τ−r+dIIt,R′t=rIt−dRRt.We rewrite the second equation of system (32) as(33)C′t=bI∗NSt−a+dC+βCt+bS∗NIt+β∫t−τtC′xdx=bI∗NSt−a+dC+βCt+bS∗NIt+β∫t−τtbI∗NSx−a+dCCx−βCx−τ+bS∗NIxdx.By the same method, from the third equation of the linearized system (32), we have(34)C′t=βCt−r+dIIt−β∫t−τtbI∗NSx−a+dCCx−βCx−τ+bS∗NIxdx.In order to prove the local asymptotic stability of the equilibrium0,0,0,0 of system (32), it is sufficient to consider the existence of a strict Lyapunov functional.LetV1t=S2t; the full-time derivative of V1t is(35)V1′t=−2bI∗N+dSS2+2aSC−2bS∗NSI.ConsideringV21t=C2t, from (33) and the inequality ab≤a2+b2/2, we obtain(36)V21′t=2bI∗NStCt−2a+dC+βC2t+2bS∗NCtIt+2βCt∫t−τtbI∗NSx−a+dCCx−βCx−τ+bS∗NIxdx≤2bI∗NStCt−2a+dC+βC2t+2bS∗NCtIt+βbI∗N+a+dC+β+bS∗NτC2t+β∫t−τtbI∗NS2x+a+dCC2x+βC2x−τ+bS∗NI2xdx.ConsideringV22t=β∫t−τt∫ztbI∗/NS2x+a+dCC2x+βC2x−τ+bS∗/NI2xdxdz, we have(37)V22′t=βτbI∗NS2t+a+dCC2t+βC2t−τ+bS∗NI2t−β∫t−τtbI∗NS2x+a+dCC2x+βC2x−τ+bS∗NI2xdx.In addition, consideringV23t=β2τ∫t−τtC2xdx and V2t=V21t+V22t+V23t, we obtain(38)V2′t≤2bI∗NStCt−2a+dC+βC2t+2bS∗NCtIt+βbI∗N+a+dC+β+bS∗NτC2t+βτbI∗NS2t+a+dCC2t+βC2t+bS∗NI2t.Similarly, from (34) and ab≤a2+b2/2, let V31t=I2t, V32t=V22t, V33t=V23t, and V3t=V31t+V32t+V33t; we obtain(39)V3′t≤2βCtIt−2r+dII2t+βbI∗N+a+dC+β+bS∗NτI2t+βτbI∗NS2t+a+dC+βC2t+bS∗NI2t.AssumeV4t=R2t and Vt=V1t+V2t+V3t+V4t. Obviously, Vt≥0 and Vt=0 if and only if St=Ct=Rt=It=0. Additionally, if St=Ct=Rt=It=0, then V′t=0. Furthermore, using ab≤a2+b2/2, we have(40)V′t≤−bI∗N+2dS−2βτbI∗N−a−bS∗NS2t−a+2dC+β−3βτa+dC+β−1+βτbI∗N+bS∗NC2t−r+2dI−3βτbS∗N−2bS∗N−β−βτbI∗N−βτa+dC+βI2t−2dR−rR2t.Therefore, when(41)bI∗N+2dS−2βτbI∗N−a−bS∗N>0,a+2dC+β−3βτa+dC+β−1+βτbI∗N+bS∗N>0,r+2dI−3βτbS∗N−2bS∗N−β−βτbI∗N−βτa+dC+β>0,2dR−r>0,then V′ is negative definite. Hence, we have the following theorem.Theorem 6.
When the following conditions,(42)τ∈0,τ1τ1=minbI∗/N+2dS−a−bS∗/N2βbI∗/N,a+2dC+β−bI∗/N−bS∗/Nβ3a+dC+β+bI∗/N+bS∗/N,r+2dI−2bS∗/N−ββ3bS∗/N+bI∗/N+a+dC+βbS∗N<minbI∗N+2dS−a,a+2dC+β−bI∗N,r+2dI−β2,2dR−r>0,hold, then equilibrium 0,0,0,0 is locally asymptotically stable for system (32), that is, the positive equilibrium E∗S∗,C∗,I∗,R∗ is also locally asymptotically stable for system (8).Example 2.
Assumea=0.01,b=0.025,β=0.058,r=0.015,τ=4,N=1.0×107,λ1=1000,λ2=200,λ3=1,dS=0.0095,dC=0.015,dI=0.055,anddR=0.012; then, system (8) reads(43)S′t=1000+0.01Ct−0.025StIt1.0×107−0.0095St,C′t=200+0.025StIt1.0×107−0.01Ct−0.058Ct−τ−0.015Ct,I′t=1+0.058Ct−τ−0.015It−0.055It,R′t=0.015It−0.012Rt.
Through calculation,S∗,C∗,I∗,R∗=107749,2416,2016,2520, bI∗/N+2dS−a−bS∗/N/2βbI∗/N=14939,a+2dC+β−bI∗/N−bS∗/N/β3a+dC+β+bI∗/N+bS∗/N=6.7593,r+2dI−2bS∗/N−β/β3bS∗/N+bI∗/N+a+dC+β=13.6719,2dR−r=0.009>0, and τ1=min14939,6.7593,13.6719=6.7593. Because condition (42) for Theorem 6 is satisfied, S∗,C∗,I∗,R∗=107749,2416,2016,2520 is locally asymptotically stable, which is also obtained from Figure 2. Note that, in Figure 2, the time histories of St,Ct,It,andRt in system (43), for τ=4, with initial functions S−τ,C−τ,I−τ=107,30,5 and S0,C0,I0=107,2000,100, and t:0−4000, respectively, all converge to 107749,2416,2016,2520 as t tends to +∞. The time histories are presented as follows.Figure 2
The time histories ofSt,Ct,It,andRt in system (43), for τ=4, with initial functions S−τ,C−τ,I−τ,R−τ=107,30,5,0, S0,C0,I0,R0=107,2000,100,0, and t:0−4000, respectively.
(a)(b)(c)(d)


### 3.4. Hopf Bifurcation for Model 2

Simply, system (8) was nondimensionalized using the scaling St−S∗⟶St,Ct−C∗⟶Ct,It−I∗⟶It,Rt−R∗⟶Rt, and t⟶t/τ; then, we have(44)S′tC′tI′tR′t=τBStCtItRt+τCSt−1Ct−1It−1Rt−1+τFSt,Ct,It,Rt,St−1,Ct−1,It−1,Rt−1,with(45)B=−bI∗N−dSa−bS∗N0bI∗N−a−dCbS∗N000−r−dI000r−dR,C=00000−β000β000000,F=−bStItNbStItN00.Assumeτ=τ0+μ, then the Hopf bifurcation value of (44) is μ=0. When ϕ=ϕ1,ϕ2,ϕ3,ϕ4T∈C−1,0,R4, let(46)Lμϕ=τ0+μBϕ0+Cϕ−1.Using the Riesz representation theorem, there is a4×4 matrix ηθ,μ,θ∈−1,0, whose elements are functions of bounded variation such that(47)Lμϕ=∫−10dηθ,μϕθ,ϕ∈C−1,0,R4.Actually, we choose(48)ηθ,μ=τ0+μB,θ=0,0,θ∈−1,0,−τ0+μC,θ=−1,which satisfies (47).Whenϕ∈C−1,0,R4, we let the operator Aμ be(49)Aμϕθ=dϕθdθ,θ∈−1,0,∫−10dηξ,μϕξ,θ=0,and(50)Rμϕθ=0,θ∈−1,0,Fμ,ϕ,θ=0,where(51)Fμ,ϕ=τ0+μ12f20ϕ120+f11ϕ10ϕ20+12f02ϕ220+16f30ϕ13012l20ϕ120+l11ϕ10ϕ20+12l02ϕ220+16l30ϕ130+12f12ϕ10ϕ220+12f21ϕ120ϕ20+16f03ϕ230+⋯+12l12ϕ10ϕ220+12l21ϕ120ϕ20+16l03ϕ230+⋯.Therefore, system (8) is equal to the following operator equation:(52)ut′=Aμut+Rμut,with ut=St,Ct,It,RtT and ut=ut+θ for θ∈−1,0.Whenψ∈C10,1,R4∗, let(53)A∗ψs=−dψsds,s∈0,1,∫−10ϕ−ξdηξ,0,s=0.And a bilinear form is(54)ψs,ϕθ=ψ¯0ϕ0−∫−10∫ξ=0θψ¯ξ−θdηθϕξdξ,with ηθ=ηθ,0. Hence, A0 and A∗ are adjoint operators.We know the eigenvalues ofA0 are ±iω0τ0 and A∗ has the same eigenvalues. Furthermore, we can verify that the eigenvectors of A0 and A∗ corresponding to the eigenvalues iω0τ0 and −iω0τ0 are vectors qθ=q1,q2Teiω0τ0θθ∈−1,0 and q∗s=Dq1∗,q2∗eiω0τ0ss∈0,1, respectively, and(55)iω0τ0I−A0q1q2q3q4eiω0τ0θ=0000,−iω0τ0I−A∗Dq1∗q2∗q3∗q4∗eiω0τ0s=0000.Namely,(56)iω0τ0I−τ0B−τ0e−iω0τ0Cq1q2q3q4=0000,−iω0τ0I−τ0BT−τ0e−iω0τ0CTq1∗q2∗q3∗q4∗=0000.It means thatq1,q2,q3, and q4 follow as(57)iω0+bI∗N+dS−abS∗N0−bI∗Niω0+a+dC+βe−iω0τ0−bS∗N00−βe−iω0τ0iω0+r+dI000−rdRq1q2q3q4=0000.If(58)iω0+bI∗N+dS−abS∗N−bI∗Niω0+a+dC+βe−iω0τ0−bS∗N0−βe−iω0τ0iω0+r+dI=0,Equation (57) has not trivial solution. Assume q4=1, and we obtain(59)q3=dRr,q2=iω0+r+dIq3βe−iω0τ0,q1=iω0+a+dC+βe−iω0τ0q2−bS∗/Nq3bI∗/N.By the same method, we have(60)−iω0+bI∗N+dS−bI∗N0−a−iω0+a+dC+βe−iω0τ0−βe−iω0τ0bS∗N−bS∗N−iω0+r+dI=0.q1∗,q2∗,q3∗, and q4∗ have no trivial solution. Assume q1∗=1; we obtain(61)q2∗=−iω0+bI∗/N+dSbI∗/N,q3∗=−a+−iω0+a+dC+βe−iω0τ0q2∗βe−iω0τ0,q4∗=0.In order to satisfy the normalized conditionq∗s,qθ=1 and the orthogonal condition q∗s,q¯θ=0, we need to determine the value D. From (54), we know(62)q∗s,qθ=q∗¯0q0−∫−10∫ξ=0θq∗¯ξ−θdηθqξdξ=D¯q1∗¯q2∗¯q3∗¯q4∗¯q1q2q3q4−∫−10∫ξ=0θD¯q1∗¯q2∗¯q3∗¯q4∗¯e−iω0τ0ξ−θdηθq1q2q3q4eiω0τ0ξdξ=D¯q1∗¯q1+q2∗¯q2+q3∗¯q3+q4∗¯q4−D¯∫−10q1∗¯q2∗¯q3∗¯q4∗¯eiω0τ0θθdηθq1q2q3q4=D¯q1∗¯q1+q2∗¯q2+q3∗¯q3+q4∗¯q4+D¯τ0e−iω0τ0βq2q3∗¯−q2∗¯=1.So, we chooseD satisfying(63)1D¯=q1∗¯q1+q2∗¯q2+q3∗¯q3+q4∗¯q4+τ0e−iω0τ0βq2q3∗¯−q2∗¯.With the same notations in [31], we first calculate the coordinates describing the center manifold C0 at μ=0. Assume ut is the solution of system (8) for μ=0.Letzt=q∗,ut and(64)Wt,θ=utθ−2Reztqθ.On the center manifoldC0, we obtain(65)Wt,θ=Wzt,z¯t,θ.With(66)Wzt,z¯t,θ=W20θz22+W11θzz¯+W02θz¯22+W30θz36+⋯and for center manifold C0, the local coordinates in the direction of q∗ and q∗¯ are z and z¯, respectively. We know if ut is real, then W is real. Here, we only consider real solutions.To the solutionut∈C0 of system (8), owing to μ=0,(67)z′t=iω0τ0zt+q∗θ,τ0FWzt,z¯t,θ+2Reztqθ=iω0τ0z+q∗¯0τ0FWzt,z¯t,0+2Reztq0=iω0τ0z+q∗¯0τ0F0z,z¯,which can be rewritten as(68)z′t=iω0τ0zt+gz,z¯.With(69)gz,z¯=q∗¯0τ0FWzt,z¯t,0+2Reztq0=g20z22+g11zz¯+g02z¯22+g21z¯2z2+⋯.From (52), (68), and (64), we can obtain(70)W′=ut′−z′q−z¯′q¯=AW−2Req∗¯0τ0F0qθ,θ∈−1,0,AW−2Req∗¯0τ0F0q0+τ0F0,θ=0,≔AW+Hz,z¯,θ,with(71)Hz,z¯,θ=H20θz22+H11θzz¯+H02θz¯22+⋯.If we expand the above series and compare the corresponding coefficients, then we can obtain(72)−H20θ=A−2iω0τ0W20θ,−H11θ=AW11θ,−H02θ=A+2iω0τ0W02θ,….By (67), assume F0z,z¯=Fz2z2/2+Fzz¯zz¯+Fz¯2z¯2/2+Fz¯2zz2z¯/2+⋯ and (68); we obtain(73)g20=q∗¯0τ0Fz2,g11=q∗¯0τ0Fzz¯,g02=q∗¯0τ0Fz¯2,g21=q∗¯0τ0Fz2z¯.By (64), we can get utθ=Wt,θ+ztqθ+z¯tq¯θ, so(74)ut=Wt,0+ztq0+z¯tq¯0=q1q2q3q4z+q1¯q2¯q3¯q4¯z¯+W200z22+W110zz¯+W020z¯22+⋯,with(75)ut=StCtItRt,ut−1=St−1Ct−1It−1Rt−1,Wt,0=W1t,0W2t,0W3t,0W4t,0,Wt,−1=W1t,−1W2t,−1W3t,−1W4t,−1.In addition, we have(76)St=q1z+q1¯z¯+W1t,0Ct=q1z+q1¯z¯+W2t,0It=q1z+q1¯z¯+W3t,0Rt=q1z+q1¯z¯+W4t,0,St−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W1t,−1Ct−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W2t,−1It−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W3t,−1Rt−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W4t,−1,From(77)gz,z¯=q∗¯0τ0F0=D¯τ0q1∗¯q2∗¯q3∗¯q4∗¯−bNSIbNSI00=bτ0D¯Nq2∗¯−q1∗¯SI,and (76), we compare the coefficients with (69) and have(78)g20=2bτ0D¯Nq2∗¯−q1∗¯q1q3,g11=bτ0D¯Nq2∗¯−q1∗¯q1q3¯+q1¯q3,g02=2bτ0D¯Nq2∗¯−q1∗¯q1¯q3¯,g21=2bτ0D¯Nq2∗¯−q1∗¯q1W1130+q3W1110+q1¯W20302+q3¯W20102+.g20,g11, and g02 can be directly computed, but we need to compute W20θ and W11θ for θ∈−1,0 in the following in order to compute g21.Firstly, we study the situation ofθ∈−1,0 of (70). By (70), we obtain(79)W′=AW−2Req¯∗0τ0F0qθ=AW−q¯∗0τ0F0qθ−q∗0τ0F¯0q¯θ=AW−gz,z¯qθ−g¯z,z¯q¯θ,and(80)W′=AW+Hz,z¯,θ.Hence,(81)Hz,z¯,θ=−gz,z¯qθ−g¯z,z¯q¯θ=−g20z22+g11zz¯+g02z¯22qθ−g¯20z¯22+g¯11zz¯+g¯02z22q¯θ=H20θz22+H11θzz¯+H02θz¯22.Therefore,(82)H20θ=−g20qθ−g¯20q¯θ=−A−2iω0τ0W20θ,and(83)H11θ=−g11qθ−g¯11q¯θ=−AW11θ.By (82), we obtain(84)AW20θ=2iω0τ0W20θ+g20qθ+g¯02q¯θ.From (76), hence,(85)W20′θ=2iω0τ0W20θ+g20q0eiω0τ0θ+g¯02q¯0e−iω0τ0θ.If we solve the above equation forW20θ, then we have(86)W20θ=K1e2iω0τ0θ−g20q0iω0τ0eiω0τ0θ−g¯02q¯03iω0τ0e−iω0τ0θ.By the same method, from (83), we can obtain(87)AW11θ=g11qθ+g¯11q¯θ.From (76), hence,(88)W11′θ=g11q0eiω0τ0θ+g¯11q¯0e−iω0τ0θ.If we solve the above equation forW11θ, then we have(89)W11θ=K2+g11q0iω0τ0eiω0τ0θ−g¯02q¯0iω0τ0e−iω0τ0θ.Secondly, to obtainK1 and K2, we study the situation of θ=0 of (70). By (70), we obtain(90)W′=AW+τ0F0−2Req¯∗0τ0F0q0=AW+Hz,z¯,0.So,(91)Hz,z¯,0=−2Req¯∗0τ0F0q0+τ0F0=−gz,z¯q0−g¯z,z¯q¯0+τ0F0=−g20z22+g11zz¯+g02z¯22q0−g¯20z¯22+g¯11zz¯+g¯02z22q¯0+−bNSIbNSI00=H200z22+H110zz¯+H020z¯22.If we put (76) into the above equation and compare the coefficient of z2/2 and zz¯, then we obtain(92)H200=−g20q0−g¯02q¯0+−2bNq1q32bNq1q300,H110=−g11q0−g¯11q¯0+−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.ByAW200=2iω0τ0W200−H200, we obtain(93)τ0BW200+τ0CW20−1=2iω0τ0W200−H200.If (86) and (92) are put into (93), then we can obtain(94)τ0BK1−g20iω0τ0τ0Bq0−g¯023iω0τ0τ0Bq¯0+e−2iω0τ0τ0CK1−g20e−iω0τ0θiω0τ0τ0Cq0−g¯02eiω0τ03iω0τ0τ0Cq¯0−2iω0τ0K1−g20iω0τ0q0−g¯023iω0τ0q¯0.=g20q0+g¯02q¯0−−2bNq1q32bNq1q300.By(95)τ0B+e−iω0τ0τ0Cq0=iω0τ0q0,τ0B+eiω0τ0τ0Cq¯0=−iω0τ0q¯0.We obtain(96)2iω0τ0I−τ0B−e−2iω0τ0τ0CK1=−2bNq1q32bNq1q300.So,(97)K1=τ02iω0I−B−e−2iω0τ0C−1−2bNq1q32bNq1q300.ByAW110=−H110, we obtain(98)τ0BW110+τ0CW11−1=−H110.If (89) and (92) are put into (98), then we can obtain(99)τ0BK2+g11iω0τ0q0−g¯11iω0τ0q¯0+τ0CK2+e−iω0τ0g11iω0τ0q0−eiω0τ0g¯11iω0τ0q¯0=g11q0+g¯11q¯0−−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.By (95), we obtain(100)τ0B+τ0CK2=−−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.So,(101)K2=−τ0B+C−1−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.FromK1 and K2, we have W11θ and W20θ. Hence, we can determine g21. On the basis of the above analysis, we know each gij in (78) is decided by the parameters and delays in system (8). Therefore, we can have the following quantities:(102)C10=i2ω0τ0g20g11−2g112−13g022+g212,μ2=−ReC10Reλ′τ0,T2=−ImC10+μ2Imλ′τ0ω0τ0,β2=2ReC10.It determines the properties of bifurcating periodic solutions near the critical valueτ0, and we get the following theorem.Theorem 7.
The following results hold for the expression of (102).(i)
The direction of the Hopf bifurcation is determined by the sign ofμ2, that is, Hopf bifurcations are forward (backward) for μ2>0<0, respectively(ii)
The stability of bifurcating periodic solutions is determined by the sign ofβ2, that is, the bifurcation periodic solutions on the center manifold are stable (unstable) for β2<0>0, respectively(iii)
The periods of the bifurcating periodic solutions is determined by the sign ofT2, that is, the periods increase (decrease) for T2>0<0, respectivelyExample 3.
Assumea=0.1,b=0.025,β=0.058,r=0.015,N=1.0×107,Λ1=5000,Λ2=300,Λ3=2,dS=1.81×10−5,dC=2.33×10−5,dI=5.48×10−5,anddR=1.95×10−5, then system (8) reads(103)S′t=5000+0.1Ct−0.025StIt1.0×107−1.81×10−5St,C′t=300+0.025StIt1.0×107−0.1Ct−0.058Ct−τ−2.33×10−5Ct,I′t=2+0.058Ct−τ−0.015It−5.48×10−5It,R′t=0.015It−1.95×10−5Rt.
Through calculation, system (103) has a positive equilibrium S∗,C∗,I∗,R∗=16040222,86339,332760,2.5597×108. a1=0.058,a2=0.1159,a3=−0.0014,a4=0.0015,a5=7.0011×10−7,a6=2.7547×10−8,Δ1=0.1739,Δ2=1.9640×10−5,andΔ3=1.4291×10−11, satisfying condition (19). Furthermore, p=0.0070,q=4.1721×10−7,andg=−4.8939×10−13, so the Hopf critical value τc=311 by the calculation process in Section 3.2.
From Theorems3 and 4, this equilibrium 16040222,86339,332760,2.5597×108 of system (103) is asymptotically stable for τ∈0,311 and unstable for τ>311, which can also be found from Figures 3–6 for the values of interest, τ=4,10,15,25, respectively, since to the best knowledge of the authors, the most incubation delay is 24 days. Note that, in Figures 3–6, four orbits St,Ct,It,andRt start from initial point 1.6×107,8.6×104,3.3×105,2.5×108 for τ=4,10,15,25, respectively, and all converge to 16040222,86339,332760,2.5597×108.
Ast tends to +∞, for the greater values τ, we do not simulate it any more.Figure 3
The time histories ofSt,Ct,It,andRt in system (103), for τ=4, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)Figure 4
The time histories ofSt,Ct,It,andRt in system (103), for τ=10, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)Figure 5
The time histories ofSt,Ct,It,andRt in system (103), for τ=15, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)Figure 6
The time histories ofSt,Ct,It,andRt in system (103), for τ=25, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)


## 3.1. Construction of Model 2

With the joint efforts of the whole country, the epidemic in some countries has been brought under control. However, due to the immigration of imported cases and the virus on imported cold-chain food, the epidemic has occurred many times. For the convenience of discussion, we will collectively refer to imported personnel and imported cold-chain food as imported from abroad.The corresponding model (Model 2) is(8)S′t=Λ1+aCt−bStItN−dSSt,C′t=Λ2+bStItN−aCt−βCt−τ−dCCt,I′t=Λ3+βCt−τ−rIt−dIIt,R′t=rIt−dRRt,where Λ1,Λ2,andΛ3 present the immigration rates of susceptible class, closely contacted class, and infected class from the foreign countries, respectively. dS,dC,dI,anddR present the death rates of susceptible individuals, closely contacted individuals, infected, and removed, respectively. Moreover, N=S+C+I+R is a variable, which keeps increasing with the input of Λ1,Λ2, and Λ3 and holds decreasing with output of the population which we do not consider. Therefore, combining the two cases into consideration, we can assume that the population is nearly a constant, the total population of a country or region.System (8) has no disease-free equilibrium 0,0,0,0 or other boundary equilibrium; there exists only the endemic equilibrium E∗S∗,C∗,I∗,R∗. From (8), it follows that(9)C∗=Λ1+Λ2−dSS∗β+dC,I∗=Λ3r+dI+βΛ1+Λ2−dSS∗r+dIβ+dC,R∗=rI∗dR,where S∗ satisfies(10)m1S∗2−m2S∗+m3=0,m1=bβdS>0,m2=bβΛ1+Λ2+bΛ3β+dC+NdSr+dIβ+dC+a>0,m3=Nr+dI.Λ1β+dC+aΛ1+Λ2>0. We analyze the root of (10), by calculation:(11)Δ=m22−4m1m3=aNdSr+dI−bβΛ1+Λ22+NdSr+dIβ+dC2+bdCΛ32+2b2βdCΛ3Λ1+Λ2+Λ3+2abβNdSΛ3r+dI+2abNdSdCΛ3r+dI+2NdSr+dIβ+dCa+bβΛ2+Λ3+bdCΛ3−bβΛ1.WhenΔ>0, equation (10) has two positive roots S1 and S2. Furthermore, if Λ1+Λ2−dSSi>0,i=1,2, then system (8) has two positive equilibria E1∗S1∗,C1∗,I1∗,R1∗ and E2∗S2∗,C2∗,I2∗,R2∗.WhenΔ>0, equation (10) has two positive roots S1 and S2. However, if Λ1+Λ2−dSS1<0, or Λ1+Λ2−dSS2<0, then system (8) has the unique equilibrium, i.e., the positive equilibrium E2∗S2∗,C2∗,I2∗,R2∗ or E1∗S1∗,C1∗,I1∗,R1∗.WhenΔ=0, equation (10) has one positive root S∗. Furthermore, if Λ1+Λ2−dSS∗>0, then system (8) has the unique equilibrium, i.e., the positive equilibrium E∗S∗,C∗,I∗,R∗.


## 3.2. Locally Asymptotically Stability for Model 2

The Jacobian matrix [29] on E∗ of sysem (8) is(12)JE∗=−bI∗N−dSa−bS∗N0bI∗N−a−dC−βe−λτbS∗N00βe−λτ−r−dI000r−dR.Hence, we can get the characteristic equation as(13)λ+dRλ+bI∗N+dSλ+a+dC+βe−λτλ+r+dR+βe−λτbI∗NbS∗N−βe−λτbS∗Nλ+bI∗N+dS−abI∗Nλ+r+dR=0.So, we have the characteristic rootλ=−dR, and the other characteristic roots satisfying(14)λ3+a1λ2e−λτ+a2λ2+a3λe−λτ+a4λ+a5e−λτ+a6=0,with(15)a1=β,a2=a+dC+r+dI+bI∗N+dS,a3=βr+dI+bI∗N+dS−bS∗N,a4=bI∗NdC+r+dI+dSa+dC+r+dI+r+dIa+dC,a5=βr+dIbI∗N+dS−dSbS∗N,a6=dSr+dIa+dC+dCbI∗Nr+dI.By (14), when τ=0, we have(16)λ3+a1+a2λ2+a3+a4λ+a5+a6=0,with(17)a1+a2=β+a+dC+r+dI+bI∗N+dS>0,a3+a4=βr+dI+bI∗N+dS−bS∗N+bI∗NdC+r+dI+dSa+dC+r+dI+r+dIa+dC,a5+a6=βr+dIbI∗N+dS−dSbS∗N+dSr+dIa+dC+dCbI∗Nr+dI.Assume(18)Δ1=a1+a2,Δ2=a1+a21a5+a6a3+a4=a1+a2a3+a4−a5+a6,Δ3=a1+a210a5+a6a3+a4a1+a200a5+a6=a5+a6Δ2.Using Routh–Hurwotz criterion, when(19)Δ1>0,Δ2>0,Δ3>0.Then, all roots of (16) have negative real parts. So, we can obtain the following result.Theorem 2.
If condition (19) holds, all roots of (16) have negative real parts. The endemic equilibrium E∗ of system (8) is asymptotically stable for τ=0 if condition (19) holds.Next, we suppose condition (19) is established. In order to determine if the real parts of some roots of (9) reaches to zero and ultimately become positive as τ varies continuously, let λ=iωω≠0 be the root; we have(20)a5−a1ω2cosωτ+a3ωsinωτ=a2ω2−a6,a3ωcosωτ+a1ω2−a5sinωτ=ω3−a4ω.Thus,a5−a1ω22+a3ω2=a2ω2−a62+ω2ω2−a42, in other words,(21)ω6+a22−a12−2a4ω4+a42−a32+2a1a5−2a2a6ω2+a62−a52=0.Assumez=ω2, equation (21) can be converted to(22)z3+pz2+qz+g=0,with p=a22−a12−2a4,q=a42−a32+2a1a5−2a2a6,andg=a62−a52.Assumehz=z3+pz2+qz+gandΓ=p2−3q. We have the following result.Lemma 1 (see [30]).
(i) Wheng<0, equation (22) has at least one positive root.(i)
Wheng≥0 and Γ=p2−3q<0, equation (22) has no positive root.\(ii)
Wheng≥0 and Γ=p2−3q≥0, equation (22) has positive root if and only if z1=−p+Γ/3>0 and hz1≤0Therefore, the following theorem is obtained by Lemma1.Theorem 3.
When(23)g<0.Equation (22) has at least one positive root.Moreover, we suppose condition (23) is established. Equation (14) has at least one positive root, i.e., equation (22) has at least one positive root. And equation (22) is a cubic, so it may have one, two, or three positive roots. Without loss of generality, we suppose that it has three positive roots, expressed as z1, z2, and z3 Then, equation (21) has three positive roots, that is, ω1=z1,ω2=z2, and ω3=z3.Assume(24)τkj=1ωkθk+2jπ,k=1,2,3,j=0,1,2,….Withθk∈0,2π, meeting(25)sinθk=ωka3a2ωk2−a6+ωk2−a4a1ωk2−a5a3ωk2+a1ωk2−a52,cosθk=a2ωk2−a6a5−a1ωk2+a3ωk2ωk2−a4a3ωk2+a1ωk2−a52,so ±iωk are a pair of purely imaginary roots of equation (14) with τ=τkj,k=1,2,3,j=0,1,2,…; obviously,(26)limj⟶∞τkj=∞,k=1,2,3.So, we can assume(27)τ0=τk00=min1≤k≤3τk0,ω0=ωk0.When condition (23) is established, owing to Theorem 2, equation (14) always has positive root and has roots having zero real part for τ=τkj, and τ0 is the minimum of τ. Therefore, for all τ∈0,τ0, the roots of equation (14) have negative real parts. We can have the following theorem.Theorem 4.
Whenτ∈0,τ0, all roots of equation (14) have negative real parts. That is, when τ∈0,τ0, the endemic equilibrium E∗of system (8) is locally asymptotically stable.Assumeλτ=ατ+iωτ is the root of equation (14) satisfying ατ0=0andωτ0=ω0. We suppose h′z0≠0 in order to ensure that ±iω0 are simple pure imaginary roots of equation (14) for τ=τ0 and λτ satisfies the transversality condition.The two sides of Equation (14) are differentiated with respect toτ; we can obtain(28)dλτdτ=a1λ2+a3λ+a5λe−λτ3λ2+2a2λ+a4−e−λτa1λ2τ+a3λτ−2a1λ+a5τ−a3=Aτ+BτiDτ+Eτi,with(29)Aτ=e−ατa1α2−ω2+a3α+a5αcosωτ+ωsinωτ−e−ατ2a1αω+a3ωωcosωτ−αsinωτ,Bτ=e−ατa1α2−ω2+a3α+a5ωcosωτ−αsinωτ+e−ατ2a1αω+a3ωαcosωτ+ωsinωτ,Dτ=3α2−ω2+2a2α−a4−e−ατcosωτa1τα2−ω2+a3ατ−2a1α+a5τ−a3+e−ατsinωτ2a1ταω+a3τω−2a1ω,Eτ=6αω+2a2ω−e−ατcosωτ2a1ταω+a3τω−2a1ω−e−ατsinωτa1τα2−ω2+a3ατ−2a1α+a5τ−a3.Puttingτ=τ0 into the above equation, separating the real and imaginary parts, and by Equation (20) and the definition of hz, we obtain(30)dReλτdτ|τ=τ0=Aτ0Dτ0+Bτ0Eτ0D2τ0+E2τ0=ω023ω04+2a22−a12−2a4ω02+a42−a32+2a1a5D2τ0+E2τ0=ω02h′ω02D2τ0+E2τ0.WhendReλτ/dτ<0 for τ<τ0 and close to τ0, we obtain equation (14) has a root λτ=αλ+iωλ with αλ. This contradicts with Theorem 4, so we have the following conclusion.Theorem 5.
Whenτ=τ0, a pair of simple purely imaginary roots of equation (14) are ±iω0 and all other roots have negative real parts. In addition, dReλτ0/dτ>0.Hence,E∗ of system (8) is asymptotically stable (unstable) when τ<τ0τ>τ0, respectively, and τ0 satisfying Equation (27).


## 3.3. Local Asymptotic Stability of Model 2

AssumeSt=S∗+S˜Ct=C∗+C˜It=I∗+I˜Rt=R∗+R˜; then, system (8) is linearized as(31)S′t=−bI∗N+dSS˜t+aC˜t−bS∗NI˜t,C˜′t=bI∗NS˜t−a+dCC˜t−βC˜t−τ+bS∗NI˜t,I˜′t=βC˜t−τ−r+dII˜t,R˜′t=rI˜t−dRR˜t.For the simplicity of the symbol, we omit the tildes; system (31) can be rewritten as(32)S′t=−bI∗N+dSSt+aCt−bS∗NIt,C′t=bI∗NSt−a+dCCt−βCt−τ+bS∗NIt,I′t=βCt−τ−r+dIIt,R′t=rIt−dRRt.We rewrite the second equation of system (32) as(33)C′t=bI∗NSt−a+dC+βCt+bS∗NIt+β∫t−τtC′xdx=bI∗NSt−a+dC+βCt+bS∗NIt+β∫t−τtbI∗NSx−a+dCCx−βCx−τ+bS∗NIxdx.By the same method, from the third equation of the linearized system (32), we have(34)C′t=βCt−r+dIIt−β∫t−τtbI∗NSx−a+dCCx−βCx−τ+bS∗NIxdx.In order to prove the local asymptotic stability of the equilibrium0,0,0,0 of system (32), it is sufficient to consider the existence of a strict Lyapunov functional.LetV1t=S2t; the full-time derivative of V1t is(35)V1′t=−2bI∗N+dSS2+2aSC−2bS∗NSI.ConsideringV21t=C2t, from (33) and the inequality ab≤a2+b2/2, we obtain(36)V21′t=2bI∗NStCt−2a+dC+βC2t+2bS∗NCtIt+2βCt∫t−τtbI∗NSx−a+dCCx−βCx−τ+bS∗NIxdx≤2bI∗NStCt−2a+dC+βC2t+2bS∗NCtIt+βbI∗N+a+dC+β+bS∗NτC2t+β∫t−τtbI∗NS2x+a+dCC2x+βC2x−τ+bS∗NI2xdx.ConsideringV22t=β∫t−τt∫ztbI∗/NS2x+a+dCC2x+βC2x−τ+bS∗/NI2xdxdz, we have(37)V22′t=βτbI∗NS2t+a+dCC2t+βC2t−τ+bS∗NI2t−β∫t−τtbI∗NS2x+a+dCC2x+βC2x−τ+bS∗NI2xdx.In addition, consideringV23t=β2τ∫t−τtC2xdx and V2t=V21t+V22t+V23t, we obtain(38)V2′t≤2bI∗NStCt−2a+dC+βC2t+2bS∗NCtIt+βbI∗N+a+dC+β+bS∗NτC2t+βτbI∗NS2t+a+dCC2t+βC2t+bS∗NI2t.Similarly, from (34) and ab≤a2+b2/2, let V31t=I2t, V32t=V22t, V33t=V23t, and V3t=V31t+V32t+V33t; we obtain(39)V3′t≤2βCtIt−2r+dII2t+βbI∗N+a+dC+β+bS∗NτI2t+βτbI∗NS2t+a+dC+βC2t+bS∗NI2t.AssumeV4t=R2t and Vt=V1t+V2t+V3t+V4t. Obviously, Vt≥0 and Vt=0 if and only if St=Ct=Rt=It=0. Additionally, if St=Ct=Rt=It=0, then V′t=0. Furthermore, using ab≤a2+b2/2, we have(40)V′t≤−bI∗N+2dS−2βτbI∗N−a−bS∗NS2t−a+2dC+β−3βτa+dC+β−1+βτbI∗N+bS∗NC2t−r+2dI−3βτbS∗N−2bS∗N−β−βτbI∗N−βτa+dC+βI2t−2dR−rR2t.Therefore, when(41)bI∗N+2dS−2βτbI∗N−a−bS∗N>0,a+2dC+β−3βτa+dC+β−1+βτbI∗N+bS∗N>0,r+2dI−3βτbS∗N−2bS∗N−β−βτbI∗N−βτa+dC+β>0,2dR−r>0,then V′ is negative definite. Hence, we have the following theorem.Theorem 6.
When the following conditions,(42)τ∈0,τ1τ1=minbI∗/N+2dS−a−bS∗/N2βbI∗/N,a+2dC+β−bI∗/N−bS∗/Nβ3a+dC+β+bI∗/N+bS∗/N,r+2dI−2bS∗/N−ββ3bS∗/N+bI∗/N+a+dC+βbS∗N<minbI∗N+2dS−a,a+2dC+β−bI∗N,r+2dI−β2,2dR−r>0,hold, then equilibrium 0,0,0,0 is locally asymptotically stable for system (32), that is, the positive equilibrium E∗S∗,C∗,I∗,R∗ is also locally asymptotically stable for system (8).Example 2.
Assumea=0.01,b=0.025,β=0.058,r=0.015,τ=4,N=1.0×107,λ1=1000,λ2=200,λ3=1,dS=0.0095,dC=0.015,dI=0.055,anddR=0.012; then, system (8) reads(43)S′t=1000+0.01Ct−0.025StIt1.0×107−0.0095St,C′t=200+0.025StIt1.0×107−0.01Ct−0.058Ct−τ−0.015Ct,I′t=1+0.058Ct−τ−0.015It−0.055It,R′t=0.015It−0.012Rt.
Through calculation,S∗,C∗,I∗,R∗=107749,2416,2016,2520, bI∗/N+2dS−a−bS∗/N/2βbI∗/N=14939,a+2dC+β−bI∗/N−bS∗/N/β3a+dC+β+bI∗/N+bS∗/N=6.7593,r+2dI−2bS∗/N−β/β3bS∗/N+bI∗/N+a+dC+β=13.6719,2dR−r=0.009>0, and τ1=min14939,6.7593,13.6719=6.7593. Because condition (42) for Theorem 6 is satisfied, S∗,C∗,I∗,R∗=107749,2416,2016,2520 is locally asymptotically stable, which is also obtained from Figure 2. Note that, in Figure 2, the time histories of St,Ct,It,andRt in system (43), for τ=4, with initial functions S−τ,C−τ,I−τ=107,30,5 and S0,C0,I0=107,2000,100, and t:0−4000, respectively, all converge to 107749,2416,2016,2520 as t tends to +∞. The time histories are presented as follows.Figure 2
The time histories ofSt,Ct,It,andRt in system (43), for τ=4, with initial functions S−τ,C−τ,I−τ,R−τ=107,30,5,0, S0,C0,I0,R0=107,2000,100,0, and t:0−4000, respectively.
(a)(b)(c)(d)


## 3.4. Hopf Bifurcation for Model 2

Simply, system (8) was nondimensionalized using the scaling St−S∗⟶St,Ct−C∗⟶Ct,It−I∗⟶It,Rt−R∗⟶Rt, and t⟶t/τ; then, we have(44)S′tC′tI′tR′t=τBStCtItRt+τCSt−1Ct−1It−1Rt−1+τFSt,Ct,It,Rt,St−1,Ct−1,It−1,Rt−1,with(45)B=−bI∗N−dSa−bS∗N0bI∗N−a−dCbS∗N000−r−dI000r−dR,C=00000−β000β000000,F=−bStItNbStItN00.Assumeτ=τ0+μ, then the Hopf bifurcation value of (44) is μ=0. When ϕ=ϕ1,ϕ2,ϕ3,ϕ4T∈C−1,0,R4, let(46)Lμϕ=τ0+μBϕ0+Cϕ−1.Using the Riesz representation theorem, there is a4×4 matrix ηθ,μ,θ∈−1,0, whose elements are functions of bounded variation such that(47)Lμϕ=∫−10dηθ,μϕθ,ϕ∈C−1,0,R4.Actually, we choose(48)ηθ,μ=τ0+μB,θ=0,0,θ∈−1,0,−τ0+μC,θ=−1,which satisfies (47).Whenϕ∈C−1,0,R4, we let the operator Aμ be(49)Aμϕθ=dϕθdθ,θ∈−1,0,∫−10dηξ,μϕξ,θ=0,and(50)Rμϕθ=0,θ∈−1,0,Fμ,ϕ,θ=0,where(51)Fμ,ϕ=τ0+μ12f20ϕ120+f11ϕ10ϕ20+12f02ϕ220+16f30ϕ13012l20ϕ120+l11ϕ10ϕ20+12l02ϕ220+16l30ϕ130+12f12ϕ10ϕ220+12f21ϕ120ϕ20+16f03ϕ230+⋯+12l12ϕ10ϕ220+12l21ϕ120ϕ20+16l03ϕ230+⋯.Therefore, system (8) is equal to the following operator equation:(52)ut′=Aμut+Rμut,with ut=St,Ct,It,RtT and ut=ut+θ for θ∈−1,0.Whenψ∈C10,1,R4∗, let(53)A∗ψs=−dψsds,s∈0,1,∫−10ϕ−ξdηξ,0,s=0.And a bilinear form is(54)ψs,ϕθ=ψ¯0ϕ0−∫−10∫ξ=0θψ¯ξ−θdηθϕξdξ,with ηθ=ηθ,0. Hence, A0 and A∗ are adjoint operators.We know the eigenvalues ofA0 are ±iω0τ0 and A∗ has the same eigenvalues. Furthermore, we can verify that the eigenvectors of A0 and A∗ corresponding to the eigenvalues iω0τ0 and −iω0τ0 are vectors qθ=q1,q2Teiω0τ0θθ∈−1,0 and q∗s=Dq1∗,q2∗eiω0τ0ss∈0,1, respectively, and(55)iω0τ0I−A0q1q2q3q4eiω0τ0θ=0000,−iω0τ0I−A∗Dq1∗q2∗q3∗q4∗eiω0τ0s=0000.Namely,(56)iω0τ0I−τ0B−τ0e−iω0τ0Cq1q2q3q4=0000,−iω0τ0I−τ0BT−τ0e−iω0τ0CTq1∗q2∗q3∗q4∗=0000.It means thatq1,q2,q3, and q4 follow as(57)iω0+bI∗N+dS−abS∗N0−bI∗Niω0+a+dC+βe−iω0τ0−bS∗N00−βe−iω0τ0iω0+r+dI000−rdRq1q2q3q4=0000.If(58)iω0+bI∗N+dS−abS∗N−bI∗Niω0+a+dC+βe−iω0τ0−bS∗N0−βe−iω0τ0iω0+r+dI=0,Equation (57) has not trivial solution. Assume q4=1, and we obtain(59)q3=dRr,q2=iω0+r+dIq3βe−iω0τ0,q1=iω0+a+dC+βe−iω0τ0q2−bS∗/Nq3bI∗/N.By the same method, we have(60)−iω0+bI∗N+dS−bI∗N0−a−iω0+a+dC+βe−iω0τ0−βe−iω0τ0bS∗N−bS∗N−iω0+r+dI=0.q1∗,q2∗,q3∗, and q4∗ have no trivial solution. Assume q1∗=1; we obtain(61)q2∗=−iω0+bI∗/N+dSbI∗/N,q3∗=−a+−iω0+a+dC+βe−iω0τ0q2∗βe−iω0τ0,q4∗=0.In order to satisfy the normalized conditionq∗s,qθ=1 and the orthogonal condition q∗s,q¯θ=0, we need to determine the value D. From (54), we know(62)q∗s,qθ=q∗¯0q0−∫−10∫ξ=0θq∗¯ξ−θdηθqξdξ=D¯q1∗¯q2∗¯q3∗¯q4∗¯q1q2q3q4−∫−10∫ξ=0θD¯q1∗¯q2∗¯q3∗¯q4∗¯e−iω0τ0ξ−θdηθq1q2q3q4eiω0τ0ξdξ=D¯q1∗¯q1+q2∗¯q2+q3∗¯q3+q4∗¯q4−D¯∫−10q1∗¯q2∗¯q3∗¯q4∗¯eiω0τ0θθdηθq1q2q3q4=D¯q1∗¯q1+q2∗¯q2+q3∗¯q3+q4∗¯q4+D¯τ0e−iω0τ0βq2q3∗¯−q2∗¯=1.So, we chooseD satisfying(63)1D¯=q1∗¯q1+q2∗¯q2+q3∗¯q3+q4∗¯q4+τ0e−iω0τ0βq2q3∗¯−q2∗¯.With the same notations in [31], we first calculate the coordinates describing the center manifold C0 at μ=0. Assume ut is the solution of system (8) for μ=0.Letzt=q∗,ut and(64)Wt,θ=utθ−2Reztqθ.On the center manifoldC0, we obtain(65)Wt,θ=Wzt,z¯t,θ.With(66)Wzt,z¯t,θ=W20θz22+W11θzz¯+W02θz¯22+W30θz36+⋯and for center manifold C0, the local coordinates in the direction of q∗ and q∗¯ are z and z¯, respectively. We know if ut is real, then W is real. Here, we only consider real solutions.To the solutionut∈C0 of system (8), owing to μ=0,(67)z′t=iω0τ0zt+q∗θ,τ0FWzt,z¯t,θ+2Reztqθ=iω0τ0z+q∗¯0τ0FWzt,z¯t,0+2Reztq0=iω0τ0z+q∗¯0τ0F0z,z¯,which can be rewritten as(68)z′t=iω0τ0zt+gz,z¯.With(69)gz,z¯=q∗¯0τ0FWzt,z¯t,0+2Reztq0=g20z22+g11zz¯+g02z¯22+g21z¯2z2+⋯.From (52), (68), and (64), we can obtain(70)W′=ut′−z′q−z¯′q¯=AW−2Req∗¯0τ0F0qθ,θ∈−1,0,AW−2Req∗¯0τ0F0q0+τ0F0,θ=0,≔AW+Hz,z¯,θ,with(71)Hz,z¯,θ=H20θz22+H11θzz¯+H02θz¯22+⋯.If we expand the above series and compare the corresponding coefficients, then we can obtain(72)−H20θ=A−2iω0τ0W20θ,−H11θ=AW11θ,−H02θ=A+2iω0τ0W02θ,….By (67), assume F0z,z¯=Fz2z2/2+Fzz¯zz¯+Fz¯2z¯2/2+Fz¯2zz2z¯/2+⋯ and (68); we obtain(73)g20=q∗¯0τ0Fz2,g11=q∗¯0τ0Fzz¯,g02=q∗¯0τ0Fz¯2,g21=q∗¯0τ0Fz2z¯.By (64), we can get utθ=Wt,θ+ztqθ+z¯tq¯θ, so(74)ut=Wt,0+ztq0+z¯tq¯0=q1q2q3q4z+q1¯q2¯q3¯q4¯z¯+W200z22+W110zz¯+W020z¯22+⋯,with(75)ut=StCtItRt,ut−1=St−1Ct−1It−1Rt−1,Wt,0=W1t,0W2t,0W3t,0W4t,0,Wt,−1=W1t,−1W2t,−1W3t,−1W4t,−1.In addition, we have(76)St=q1z+q1¯z¯+W1t,0Ct=q1z+q1¯z¯+W2t,0It=q1z+q1¯z¯+W3t,0Rt=q1z+q1¯z¯+W4t,0,St−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W1t,−1Ct−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W2t,−1It−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W3t,−1Rt−1=e−iω0τ0q1z+eiω0τ0q1¯z¯+W4t,−1,From(77)gz,z¯=q∗¯0τ0F0=D¯τ0q1∗¯q2∗¯q3∗¯q4∗¯−bNSIbNSI00=bτ0D¯Nq2∗¯−q1∗¯SI,and (76), we compare the coefficients with (69) and have(78)g20=2bτ0D¯Nq2∗¯−q1∗¯q1q3,g11=bτ0D¯Nq2∗¯−q1∗¯q1q3¯+q1¯q3,g02=2bτ0D¯Nq2∗¯−q1∗¯q1¯q3¯,g21=2bτ0D¯Nq2∗¯−q1∗¯q1W1130+q3W1110+q1¯W20302+q3¯W20102+.g20,g11, and g02 can be directly computed, but we need to compute W20θ and W11θ for θ∈−1,0 in the following in order to compute g21.Firstly, we study the situation ofθ∈−1,0 of (70). By (70), we obtain(79)W′=AW−2Req¯∗0τ0F0qθ=AW−q¯∗0τ0F0qθ−q∗0τ0F¯0q¯θ=AW−gz,z¯qθ−g¯z,z¯q¯θ,and(80)W′=AW+Hz,z¯,θ.Hence,(81)Hz,z¯,θ=−gz,z¯qθ−g¯z,z¯q¯θ=−g20z22+g11zz¯+g02z¯22qθ−g¯20z¯22+g¯11zz¯+g¯02z22q¯θ=H20θz22+H11θzz¯+H02θz¯22.Therefore,(82)H20θ=−g20qθ−g¯20q¯θ=−A−2iω0τ0W20θ,and(83)H11θ=−g11qθ−g¯11q¯θ=−AW11θ.By (82), we obtain(84)AW20θ=2iω0τ0W20θ+g20qθ+g¯02q¯θ.From (76), hence,(85)W20′θ=2iω0τ0W20θ+g20q0eiω0τ0θ+g¯02q¯0e−iω0τ0θ.If we solve the above equation forW20θ, then we have(86)W20θ=K1e2iω0τ0θ−g20q0iω0τ0eiω0τ0θ−g¯02q¯03iω0τ0e−iω0τ0θ.By the same method, from (83), we can obtain(87)AW11θ=g11qθ+g¯11q¯θ.From (76), hence,(88)W11′θ=g11q0eiω0τ0θ+g¯11q¯0e−iω0τ0θ.If we solve the above equation forW11θ, then we have(89)W11θ=K2+g11q0iω0τ0eiω0τ0θ−g¯02q¯0iω0τ0e−iω0τ0θ.Secondly, to obtainK1 and K2, we study the situation of θ=0 of (70). By (70), we obtain(90)W′=AW+τ0F0−2Req¯∗0τ0F0q0=AW+Hz,z¯,0.So,(91)Hz,z¯,0=−2Req¯∗0τ0F0q0+τ0F0=−gz,z¯q0−g¯z,z¯q¯0+τ0F0=−g20z22+g11zz¯+g02z¯22q0−g¯20z¯22+g¯11zz¯+g¯02z22q¯0+−bNSIbNSI00=H200z22+H110zz¯+H020z¯22.If we put (76) into the above equation and compare the coefficient of z2/2 and zz¯, then we obtain(92)H200=−g20q0−g¯02q¯0+−2bNq1q32bNq1q300,H110=−g11q0−g¯11q¯0+−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.ByAW200=2iω0τ0W200−H200, we obtain(93)τ0BW200+τ0CW20−1=2iω0τ0W200−H200.If (86) and (92) are put into (93), then we can obtain(94)τ0BK1−g20iω0τ0τ0Bq0−g¯023iω0τ0τ0Bq¯0+e−2iω0τ0τ0CK1−g20e−iω0τ0θiω0τ0τ0Cq0−g¯02eiω0τ03iω0τ0τ0Cq¯0−2iω0τ0K1−g20iω0τ0q0−g¯023iω0τ0q¯0.=g20q0+g¯02q¯0−−2bNq1q32bNq1q300.By(95)τ0B+e−iω0τ0τ0Cq0=iω0τ0q0,τ0B+eiω0τ0τ0Cq¯0=−iω0τ0q¯0.We obtain(96)2iω0τ0I−τ0B−e−2iω0τ0τ0CK1=−2bNq1q32bNq1q300.So,(97)K1=τ02iω0I−B−e−2iω0τ0C−1−2bNq1q32bNq1q300.ByAW110=−H110, we obtain(98)τ0BW110+τ0CW11−1=−H110.If (89) and (92) are put into (98), then we can obtain(99)τ0BK2+g11iω0τ0q0−g¯11iω0τ0q¯0+τ0CK2+e−iω0τ0g11iω0τ0q0−eiω0τ0g¯11iω0τ0q¯0=g11q0+g¯11q¯0−−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.By (95), we obtain(100)τ0B+τ0CK2=−−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.So,(101)K2=−τ0B+C−1−bNq1q3¯+q1¯q3bNq1q3¯+q1¯q300.FromK1 and K2, we have W11θ and W20θ. Hence, we can determine g21. On the basis of the above analysis, we know each gij in (78) is decided by the parameters and delays in system (8). Therefore, we can have the following quantities:(102)C10=i2ω0τ0g20g11−2g112−13g022+g212,μ2=−ReC10Reλ′τ0,T2=−ImC10+μ2Imλ′τ0ω0τ0,β2=2ReC10.It determines the properties of bifurcating periodic solutions near the critical valueτ0, and we get the following theorem.Theorem 7.
The following results hold for the expression of (102).(i)
The direction of the Hopf bifurcation is determined by the sign ofμ2, that is, Hopf bifurcations are forward (backward) for μ2>0<0, respectively(ii)
The stability of bifurcating periodic solutions is determined by the sign ofβ2, that is, the bifurcation periodic solutions on the center manifold are stable (unstable) for β2<0>0, respectively(iii)
The periods of the bifurcating periodic solutions is determined by the sign ofT2, that is, the periods increase (decrease) for T2>0<0, respectivelyExample 3.
Assumea=0.1,b=0.025,β=0.058,r=0.015,N=1.0×107,Λ1=5000,Λ2=300,Λ3=2,dS=1.81×10−5,dC=2.33×10−5,dI=5.48×10−5,anddR=1.95×10−5, then system (8) reads(103)S′t=5000+0.1Ct−0.025StIt1.0×107−1.81×10−5St,C′t=300+0.025StIt1.0×107−0.1Ct−0.058Ct−τ−2.33×10−5Ct,I′t=2+0.058Ct−τ−0.015It−5.48×10−5It,R′t=0.015It−1.95×10−5Rt.
Through calculation, system (103) has a positive equilibrium S∗,C∗,I∗,R∗=16040222,86339,332760,2.5597×108. a1=0.058,a2=0.1159,a3=−0.0014,a4=0.0015,a5=7.0011×10−7,a6=2.7547×10−8,Δ1=0.1739,Δ2=1.9640×10−5,andΔ3=1.4291×10−11, satisfying condition (19). Furthermore, p=0.0070,q=4.1721×10−7,andg=−4.8939×10−13, so the Hopf critical value τc=311 by the calculation process in Section 3.2.
From Theorems3 and 4, this equilibrium 16040222,86339,332760,2.5597×108 of system (103) is asymptotically stable for τ∈0,311 and unstable for τ>311, which can also be found from Figures 3–6 for the values of interest, τ=4,10,15,25, respectively, since to the best knowledge of the authors, the most incubation delay is 24 days. Note that, in Figures 3–6, four orbits St,Ct,It,andRt start from initial point 1.6×107,8.6×104,3.3×105,2.5×108 for τ=4,10,15,25, respectively, and all converge to 16040222,86339,332760,2.5597×108.
Ast tends to +∞, for the greater values τ, we do not simulate it any more.Figure 3
The time histories ofSt,Ct,It,andRt in system (103), for τ=4, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)Figure 4
The time histories ofSt,Ct,It,andRt in system (103), for τ=10, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)Figure 5
The time histories ofSt,Ct,It,andRt in system (103), for τ=15, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)Figure 6
The time histories ofSt,Ct,It,andRt in system (103), for τ=25, with initial functions S−τ,C−τ,I−τ,R−τ=S0,C0,I0,R0=1.6×107,8.6×104,3.3×105,2.5×108 and t:68000−70000, respectively.
(a)(b)(c)(d)


## 4. Conclusions

In this study, we divided people into four categories:St susceptible, Ct closely contacted, It infective, and Rt removed according to the current situation in the middle and late stages of COVID-19, investigating the SCIR models with immigration and without immigration, and obtained the following main results:(1)
Ifb<r/2, then the disease-free equilibrium E0S0,0,0 of Model 1 is globally asymptotically stable(2)
IfΔ1>0,Δ2>0,Δ3>0, the endemic equilibrium E∗ of Model 2 is asymptotically stable for τ∈0,τ0(3)
If condition (42) holds, then the positive equilibrium E∗S∗,C∗,I∗,R∗ of Model 2 is locally asymptotically stable by constructing Lyapunov function(4)
By the bifurcation formulae, in [31], we qualitatively analyze the properties for the occurring Hopf bifurcations of Model 2Moreover, considering neither the immigration nor the death rate during the outbreak, Model 1 has only the disease-free equilibrium. However, considering both the immigration and the death rate during late stage, Model 2 has only one or two endemic equilibria. Both indicate the importance of finding closely contacted and overseas imports on epidemic control, respectively.



---

*Source: 9182830-2022-05-11.xml*