# Complex Behaviors of Epidemic Model with Nonlinear Rewiring Rate

Authors: Ding Fang; Yongxin Zhang; Wendi Wang
Journal: Complexity
 (2020)

---

## Abstract

An SIS propagation model with the nonlinear rewiring rate on an adaptive network is considered. It is found by bifurcation analysis that the model has the complex behaviors which include the transcritical bifurcation, saddle-node bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcation. Especially, a bifurcation curve with “S” shape emerges due to the nonlinear rewiring rate, which leads to multiple equilibria and twice saddle-node bifurcations. Numerical simulations show that the model admits a homoclinic bifurcation and a saddle-node bifurcation of the limit cycle.


---

## Body

## 1. Introduction

Mathematical models of epidemic propagations on networks have been proposed and studied by many papers (see [1–3] and the references cited therein) due to the importance in applications. Normile [4] described the network of SARS transmissions through hospitals and transportation. Zhang et al. [5, 6] modeled the transmission of H7N9 through transportation networks in China, respectively. Wang et al. [7–9] studied the impact of virus transmission on the network of mobile devices, respectively.In recent years, there has been much interest in studying epidemic propagation on adaptive networks [10–15]. This is because behavior changes during the course of epidemic propagation are quite common and have the strong influence on the progression of disease. Indeed, an individual who links to an infective node tends to break the linkage and connect to a safer neighbor. For these reasons, Gross et al. [16] extended the classical SIS model in a static network to the adaptive network whose structure changes with the progression of epidemic disease. This model is important because it simulates the spread of epidemic disease in the network with the behavior changes and captures the key features of disease evolutions. There have been a number of studies about epidemic propagation on the adaptive network in this direction (see [11–13] and the references cited therein).The following model, which originates from [16, 17], is presented in [18]:(1)S˙=−τSI+γI,I˙=τSI−γI,SI.=−γSI+γII−τSI−τISI+τSSI−ωSI,II.=−2γII+2τSI+2τISI,SS.=2γSI−2τSSI+2ωSI,where S and I denote the expected numbers of susceptible and infectious nodes, respectively. SI represents the expected number of SI edges, which means a susceptible node contacts with an infected node; SS and II have similar meanings. SSI and ISI denote the expected numbers of SSI and ISI triples [17]. An infected node can infect any susceptible neighbor at rate τ and recover to susceptible state at rate γ. The rewiring term is ωSI, which means that the susceptible individual disconnects with an infected neighbor at rate ω called as rewiring rate and randomly selects a susceptible individual to reconnect immediately. Juher et al. [15] analyzed the dynamics of model (1). Bodó and Simon [14] gave further mathematical analysis about the stability of equilibria and Hopf bifurcation for model (1). By numerical simulations, they found that the model has the complex behaviors such as the Bogdanov–Takens bifurcation and double limit cycles. A very important work was given by Zhang et al. [19] recently. They obtained the complete analytical conditions for bifurcations of the model that reveals the mechanism of various complex phenomena in model (1).Notice that humans may take different behaviors as the risk of disease transmission varies [20, 21]. For example, during the SARS outbreak in 2003, the government and public opinion urged people to take protection measures by reducing infectious contacts, and people adopt the stronger responses when the reported cases were increased [22, 23]. Thus, social contacts tend to disconnect infectious nodes and turn to safe connections. These activities result in the rise of rewiring rate when the number of infectious individuals increases. For this reason, we modify the rewiring rate in model (1) from ωSI to ωSI+β˜ISI, where β˜ is a nonnegative coefficient that describes the influence of infective number on the rewiring rate. Then, model (1) is extended to(2)S˙=−τSI+γI,I˙=τSI−γI,SI.=−γ+τSI+γII−τISI+τSSI−ω+β˜ISI,II.=−2γII+2τSI+2τISI,SS.=2γSI−2τSSI+2ω+β˜ISI.It is clear that model (2) extends model (1) by introducing the infection-dependent response coefficient β˜. This is important for a larger epidemic outbreak where human behaviors are crucially affected by the number of infected individuals, which is an index of infection risk.The objective of this paper is to present the analytical analysis for bifurcations of model (2). We find that the model exhibits the more complicated behaviors due to the nonlinear rewiring rate. Besides the transcritical bifurcation, saddle-node bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcation in model (1), there exists a bifurcation curve with “S” shape, which is attributed to the nonlinear rewiring rate and leads to multiple stable equilibria and twice saddle-node bifurcations. More importantly, the disease can persist even though the basic reproduction number is less than unity when the bifurcation at the disease-free equilibrium is forward. Moreover, we find a scenario where a stable periodic solution coexists with a stable equilibrium.This paper is organized as follows. In Section2, we simplify the model by the technique of moment closure and study the distribution of equilibrium points of the reduced model. The detailed analysis of stability of the steady states and bifurcations of the system is presented in Section 3. Numerical simulations are presented in Section 4 to reveal more interesting global behaviors of the model and illustrate the effect of adaptive rewiring on epidemic spread. The paper ends with discussions in Section 5.


## 2. Model Reduction and Steady States

Notice that the time evolutions of pairs in equation (2) depend on triples, which makes the analysis of model (2) impossible. A common method is to use the moment closure technique to obtain the closed form of these equations. Here, we follow Zhang et al. [19] to assume the degree distribution of the network obeys the Poisson distribution. In this case, the triples can be approximated by pairs and singles according to Keeling [24]:(3)SSI=SSSIS,ISI=SISIS.The first equality in (3) means that SSI links are determined by the product of number of SS links and the probability that a given susceptible node has an SI link. The second equality in (3) has the similar meaning.LetN be the number of vertices of the network and n be the average degree of the graph. Since the derivative 2SI+SS+II equals zero, we have 2SI+SS+II=nN. Then, let us nondimensionalize the system by the transformations(4)PS=SN,PI=IN,PSI=2SInN,PSS=SSnN,PII=IInN,to obtain(5)P˙I=τn2PSI−γPI,P˙SI=−γ+τ+ωPSI+2γ1−PSI−PSS+τnPSS−PSI2PSI1−PI−βPIPSI,P˙SS=γPSI−τnPSSPSI1−PI+ωPSI+βPIPSI,in which, the identities of PS+PI=1 and PSI+PSS+PII=1 are used, and β=β˜N. In view of biological background, we assume that the parameters τ,γ, and n are positive constants, and the parameters ω and β are nonnegative constants. It is easy to verify the domain(6)Ω=PI,PSI,PSS0≤PSI,0≤PSS,0≤PI<1,PSI+PSS≤1.is a positive invariant set of system (5). In the rest of this paper, we study the dynamical behaviors of the system in Ω.The disease-free steady state (DEF) of system (5) is E0=0,0,1. With the help of the computation method in [25], we get the basic reproduction number:(7)R0=τnγ+ω.The endemic steady states are determined by(8)0=τn2PSI−γPI,(9)0=−γPSI+2γ1−PSI−PSS−τPSI−τn2PSIPSI1−PI+τnPSSPSI1−PI−ωPSI−βPIPSI,(10)0=γPSI−τnPSSPSI1−PI+ωPSI+βPIPSI.Obviously, from equations (8) and (10), we have(11)PSI=2γPIτn,PSS=−βτnPI2+−γ−ω+βτnPI+γ+ωτn.By substituting these expressions into equation (9), we can deduce PI which satisfies the following equation:(12)FPI≔API3+ℬPI2+CPI+D=0,where(13)A=β,ℬ=ω−2β−τ,C=−2ω+β+τ+τn,D=γ+ω1−R0.Set(14)A=ℬ2−3Ac,B=ℬC−9AD,C=C2−3ℬD.Motivated by [26], we define the discriminant of the cubic function FPI by(15)Δ0≔B2−4AC.Notice thatA≥0 if Δ0=0. Similar to [26], we see that the positive roots of equation (12) are classified by the signs of Δ0,ℬ,C, and D. Then, we state the distribution of equilibria of equation (12) in three cases.C1:R0>1:(1)
IfC≤0, then FPI has a unique positive root.(2)
Suppose thatC>0 and ℬ≥0. Then, FPI has a unique positive root.(3)
Suppose thatC>0 and ℬ<0. Then, FPI has a unique positive root if Δ0>0, has a positive root of multiplicity 2 and another positive root if Δ0=0 and A>0, has a positive root of multiplicity 3 if Δ0=A=0, and has three distinct positive roots if Δ0<0.C2:R0<1:(1)
Suppose thatC≥0 and ℬ≥0. Then, FPI has no positive root.(2)
Suppose thatC≥0 and ℬ<0. Then, FPI has no positive root if Δ0>0, has a positive root of multiplicity 2 if Δ0=0, and has two positive roots if Δ0<0.(3)
Suppose thatC<0. Then, FPI has no positive root if Δ0>0, has a positive root of multiplicity 2 if Δ0=0, and has two positive roots if Δ0<0.C3:R0=1:(1)
Suppose thatC<0. Then, FPI has a unique positive root.(2)
Suppose thatC=0 and ℬ<0. Then, FPI has a unique positive root.(3)
Suppose thatC≥0 and ℬ≥0. Then, FPI has no positive root.(4)
Suppose thatC>0 and ℬ<0. Then, FPI has no positive root if Δ0>0, has a positive root of multiplicity 2 if Δ0=0, and has two positive roots if Δ0<0.Based on the above analysis about the roots of equation (12), we obtain the distributions of endemic equilibria for system (5), which are shown in Table 1. It is obvious that the distributions of the endemic equilibria are complex, which means that there exist rich bifurcations in system (5). We will elucidate these bifurcations in Section 3.Table. 1
Distribution of the endemic equilibria.
Range ofR0
Other conditions
Equilibria of system (5)
R
0
<
1
C
≥
0
ℬ
≥
0
No positive equilibrium
C
≥
0
ℬ
<
0
Δ
0
<
0
Two positive equilibria
C
≥
0
ℬ
<
0
Δ
0
=
0
One positiveequilibrium2
C
≥
0
ℬ
<
0
Δ
0
>
0
No positive equilibrium
C
<
0
Δ
0
<
0
Two positive equilibria
C
<
0
Δ
0
=
0
One positiveequilibrium2
C
<
0
Δ
0
>
0
No positive equilibrium
R
0
>
1
C
≤
0
One positive equilibrium
C
>
0
ℬ
≥
0
One positive equilibrium
C
>
0
ℬ
<
0
Δ
0
<
0
Three positive equilibria
C
>
0
ℬ
<
0
Δ
0
>
0
One positive equilibrium
C
>
0
ℬ
<
0
Δ
0
=
0 A>0
Two positiveequilibria1,2
C
>
0
ℬ
<
0
Δ
0
=
0 A=0
One positiveequilibrium3
R
0
=
1
C
<
0
One positive equilibrium
C
>
0
ℬ
<
0
Δ
0
>
0
No positive equilibrium
C
>
0
ℬ
<
0
Δ
0
=
0
One positiveequilibrium2
C
>
0
ℬ
<
0
Δ
0
<
0
Two positive equilibria
C
≥
0
ℬ
≥
0
No positive equilibrium
C
=
0
ℬ
<
0
One positive equilibrium
Two positiveequilibria1,2 mean that there are a simple equilibrium and another equilibrium of multiplicity 2. One positive equilibrium2 means that there is an equilibrium of multiplicity 2. One positive equilibrium3 means that there is an equilibrium of multiplicity 3. The remaining cases represent a simple equilibrium.


## 3. Stability and Bifurcation Analysis

In this section, we study the stability of steady states to get the sufficient and necessary conditions of backward bifurcation and Hopf bifurcation. Besides, we show that the system undergoes the Bogdanov–Takens bifurcation.


### 3.1. Backward Bifurcation

For the stability of the disease-free steady state (DEF), we have the following proposition.Proposition 1.
The DEF is asymptotically stable ifR0<1 and is unstable if R0>1.Proof.
It is easy to find that−γ is an eigenvalue of the Jacobian matrix at DEF. The other two eigenvalues λ1 and λ2 are the roots of(16)λ2+λ3γ−τn+ω+τ+2γγ−τn+ω=0.
IfR0<1, then τn<γ+ω, which leads to 3γ−τn+ω+τ>0. Thus, both λ1 and λ2 have negative real parts, and therefore, the DEF is asymptotically stable. If R0>1, then one root of (16) is positive, which means that the DEF is unstable.
There exist two endemic steady states whenR0<1. This implies the possibility that the system has a backward bifurcation. Let us now figure out the conditions for the backward bifurcation to appear. Set(17)a0=2γnγ+ω−β+γ+ω2n+1γ+ωn−1,b0=2γn22γn+γ+ω.
By the technique of central manifold [27], we obtain the following theorem.Theorem 1.
System (5) exhibits a backward bifurcation at DEF when R0 crosses unity if a0>0 and a forward bifurcation if a0<0.
The proof of Theorem 1 is deferred to Appendix A. It is noting that the coefficientβ has a large impact on the backward bifurcation although it has no effect on R0. From Table 1, we find that the system may undergo a saddle-node bifurcation when the two equilibria coalesce into one equilibrium. Thus, we derive the conditions for the saddle-node bifurcation. Interestingly, we show that the system has a forward bifurcation curve that looks like “S” shape.


### 3.2. Saddle-Node and Hopf Bifurcations

In this section, we first present the conditions that system (5) undergoes a saddle-node bifurcation. The Jacobian matrix of the system at an endemic state E∗=PI∗,PSI∗,PSS∗ is(18)JE∗=−γτn202γPI∗−2βPI∗2+3β−ω−2γPI∗−β+ω+γτn1−PI∗2−2γPI∗1−PI∗−τ−2γ2γPI∗1−PI∗−12γPI∗−2βPI∗−γ−ω+βτn1−PI∗0−2γPI∗1−PI∗.The eigenvalues ofJE∗ are the roots of(19)λ3+G2PI∗λ2+G1PI∗λ+G0PI∗=0,where(20)G2PI∗=−γPI∗−τPI∗+3γ+τ−1+PI∗=γPI∗+τ1−PI∗+3γ1−PI∗>0,(21)G1PI∗=γ2βPI∗3+−3β+ω−τPI∗2+β+3γ−ωPI∗+2γ+τ−1+PI∗2,(22)G0PI∗=2γ2PI∗F′PI∗1−PI∗2.Thus, the matrixJE∗ admits a zero eigenvalue if and only if G0PI∗=0, i.e., F′PI∗=0. Furthermore, zero is a simple eigenvalue if G1PI∗≠0, and Table 1 shows that E∗ is an endemic equilibrium of multiplicity 2 or 3 when Δ0=0. Consequently, we can state the following theorem.Theorem 2.
AssumeΔ0=0. Then, system (5) has an endemic equilibrium E^=P^I,P^SI,P^SS of multiplicity 2 if one of the following is satisfied:(a)
R0<1, C≥0, and ℬ<0(b)
R0<1 and C<0(c)
R0>1, C>0, ℬ<0, and A>0(d)
R0=1, C>0, and ℬ<0
Furthermore, system (5) undergoes a saddle-node bifurcation at E^ if G1PI^≠0.
The proof of Theorem 2 is omitted because it is similar to that in [19, 26, 28].
From Theorems1 and 2, we find that the rewiring rate β has remarkable impact on behaviors of system (5). To illustrate this, we depict the bifurcation diagrams of system (5) with different β in Figure 1, where the basal parameters are γ=0.1,n=6,ω=2.9. In Figure 1, when the rewiring rate β increases from zero, the type of bifurcation changes. For 0<β<2.3, the bifurcation at DEF is backward, and there exists a saddle-node bifurcation SN (see the solid bifurcation curve of Figure 1). At β=2.3, which leads the critical condition a0=0, the DEF becomes an equilibrium of multiplicity 3 (see the dashed bifurcation curve of Figure 1). As β increases further, the bifurcation at DEF becomes forward, and there exist two saddle-node bifurcations SN1 and SN2 for β>2.3. It is notable that the equilibrium SN2 coincides with DEF at β=2.3 if R0=1. In a word, we find that the coefficient β can induce the bifurcation curve which looks like “S” shape (see the dotted dash bifurcation curve and dotted bifurcation curve of Figure 1). By Theorems 1 and 2, when β>2.3, the bifurcation curve will remain S-shaped. More importantly, the disease can persist even though R0<1 when the bifurcation at DEF is forward (see the dotted dash bifurcation curve of Figure 1).
By (21), we know G2PI∗>0. Using the Routh–Hurwitz criterion and invoking the method to examine the Hopf bifurcation in [29], we see that the equilibrium point PI∗,PSI∗,PSS∗ is stable if(23)G0PI∗>0,G2PI∗G1PI∗−G0PI∗>0.
Moreover, the Hopf bifurcation occurs when(24)G2PI∗G1PI∗−G0PI∗=0,G1PI∗>0.Figure 1
Bifurcation diagrams of system (5) where SN, SN1, and SN2 represent saddle-node points. The parameter values are γ=0.1, n=6, ω=2.9, and β=0 (solid), 2.3 (dashed), 3.5 (dotted dash), or 15 (dotted).


### 3.3. Bogdanov–Takens Bifurcation

WhenΔ0=0 and one of the conditions (a), (b), (c), and (d) in Theorem 2 is satisfied, system (5) has an endemic equilibrium E^=P^I,P^SI,P^SS of multiplicity 2. It is easy to verify A≠0 under these conditions. Furthermore, we see that the components of E^ are given by(25)P^I=−B2A,P^SI=−γBτnA,P^SS=−βB24τnA2+γ+ω−β2τnAB+γ+ωτn.IfG1PI^≠0, by Theorem 2, E^ is a saddle node. If G1PI^=0, it is obvious that the eigenvalues of the Jacobian matrix at E^ are 0, 0, −G2PI^<0. We show in the following E^ is a Bogdanov–Takens point.Theorem 3.
Assume one of the conditions (a), (b), (c), and (d) in Theorem2 is satisfied. If Δ0=0 and G1PI^=0, then E^ is a Bogdanov–Takens point, and system (5) localized at E^ is topologically equivalent to(26)μ˙=ν,ν˙=μ2+μν+Oμ,ν3.Proof.
We just outline the main steps and leave the cumbersome details in Appendix B. First, let us bringE^ to the origin by substituting x=PI−P^I,y=PSI−P^SI,z=PSS−P^SS. Then, system (5) becomes(27)x˙=−γx+τn2y,y˙=2γP^IτnβP^I+γ+ω1−P^I−γP^I1−P^I2−βx+−2γP^I1−P^I−τ−2γy+2γP^I1−P^I−1z−τnP^SIP^SI−2P^SS21−P^I3x2−βP^I2+τnP^SI−τnP^SS−2βP^I+β1−P^I2xy+τnP^SI1−P^I2xz−τn21−P^Iy2+τn1−P^Iyz+Ox,y,z3,z˙=2γP^Iτn−βP^I+γ+ω1−P^I+βx−2γP^I1−P^Iz−τnP^SSP^SI1−P^I3x2+βP^I2−τnP^SS−2βP^I+β1−P^I2xy−τnP^SI1−P^I2xz−τn1−P^Iyz+Ox,y,z3.
BecauseE^ is an equilibrium of multiplicity of 2, we can get G0PI^=0 from (21). From (19), we see that the Jacobian matrix at E^=P^I,P^SI,P^SS has three eigenvalues: λ1=λ2=0 and λ3=−G2P^I<0 if G1P^I=0. The generalized eigenvectors corresponding to 0 are(28)V1=−τnγ+2P^I−1β+ω,−2γγ+2P^I−1β+ω,1′,V2=τn2βP^I−β+γ+ωγ,0,−P^I+12γP^I′,and the eigenvectors corresponding to λ3 are(29)V3=−τnτ+γP^I−τ−3γ2γP^I2βP^I−β+γ+ω,τ+γP^I−τ−3γτ−2γP^I−τ−2γγP^I2βP^I−β+γ+ωP^I−1,1′.
Here, the symbol′ means the transpose of a vector. Then, the matrix V=V1,V2,V3 is used for the nonsingular linear transformation:(30)xyz=Vμνξ.
As a result, we get(31)μ˙=ν+L20μ2+L11μν+L02ν2+ξ⋅Oμ,ν,ξ+Oμ,ν,ξ3,ν˙=M20μ2+M11μν+M02ν2+ξ⋅Oμ,ν,ξ+Oμ,ν,ξ3,ξ˙=−G2P^Iξ+Oμ,ν,ξ2.
The expressions of the coefficients are shown in Appendix B. By the center manifold theorem [30], system (31) restricted to the center manifold becomes(32)μ˙=ν+L20μ2+L11μν+L02ν2+Oμ,ν3,ν˙=M20μ2+M11μν+M02ν2+Oμ,ν3.
We use the near-identity transformation:(33)μ=X+12L11+M02X2+L02XY+OX,Y3,ν=Y−L20X2+M02XY+OX,Y3,and then rewrite X, Y into μ and ν to get(34)μ˙=ν,ν˙=M20μ2+M11+2L20μν+Oμ,ν3.
It is proved in Appendix B thatM20>0 and M11+2L20>0. If we rescale the time and variables by(35)μ,ν,t⟶M20M11+2L202μ,M202M11+2L203ν,M11+2L20M20t,then system (34) is topologically equivalent to normal form (26).
From Theorem3, we need two parameters to unfold the codimension two singularity of E^. Firstly, let τ,ω,n,γ1,β1 satisfy the conditions in Theorem 3. Then, we have(36)β1=−3ω+3τP^I2+2ω−τ2P^I−13P^I2−2,γ1=τP^I−123P^I2−2.
In the following, we consider system (5) for parameters β,γ in a neighborhood β1,γ1. Substituting(37)β=β1+ε1,γ=γ1+ε2,into system (5), we obtain a perturbed system(38)P˙I=τn2PSI−γ1+ε2PI,P˙SI=−γ1+ϵ2+τ+ωPSI+2γ1+ϵ21−PSI−PSS+τnPSS−PSI2PSI1−PI−β1+ε1PIPSI,P˙SS=γ1+ε2PSI−τnPSSPSI1−PI+ωPSI+β1+ε1PIPSI,where ε1,ε2 is in the neighborhood of 0,0.Theorem 4.
If parametersε=ε1,ε2 are sufficiently small and τ,ω,n,γ1,β1 satisfy the conditions in Theorem 3, then system (44) admits a generic unfolding of the cusp singularity of codimension 2.Proof.
We present only the main steps to illustrate the proof and omit the details of calculation because they are similar to those in [19, 26]. Shift E^ to the origin by substituting x=PI−P^I,y=PSI−P^SI,z=PSS−P^SS. Then, using the method in Appendix B, we get(39)n=−3ω+3τP^I4+6ωP^I3−ω+5τP^I2+4τ−4ωP^I+2ω−ττ2P^I−13P^I2−2.
Next, we apply transformation (30) to system (38) to obtain(40)μ˙=ν+l00010ε1+l00001ε2+L20μ2+L11μν+L02ν2+∑i+j+k+p+q≥2i+j≠2ifi+j+k+p+q=2lijkpqμiνjξkε1pε2q,ν˙=m00010ε1+m00001ε2+M20μ2+M11μν+M02ν2+∑i+j+k+p+q≥2i+j≠2ifi+j+k+p+q=2mijkpqμiνjξkε1pε2q,ξ˙=−G2P^Iξ+g00010ε1+g00001ε2+∑i+j+k+p+q≥2gijkpqμiνjξkε1pε2q.
According to the method in [19, 26], the local center manifold has the following form:(41)ξ=g00010ε1+g00001ε2G2P^I+∑i+j+p+q≥2hijpqμiνjε1pε2q.
Then, system (40) restricted to this center manifold is(42)μ˙=ν+l00010ε1+l00001ε2+L20μ2+L11μν+L02ν2+∑i+j+p+q≥2i+j≠2ifi+j+p+q=2lijpqμiνjε1pε2q,ν˙=m00010ε1+m00001ε2+M20μ2+M11μν+M02ν2+∑i+j+p+q≥2i+j≠2ifi+j+p+q=2mijpqμiνjε1pε2q.
We use the near-identity transformation(43)μ=1+L1010ε1+L1001ε2+Oε2X+L0010ε1+L0001ε2+Oε2+Oε2Y+12L11+M02+OεX2+L02+OεXY+OεY2+OX,Y3,ν=1+M0110ε1+M0101ε2+Oε2Y−l00010ε1+l00001ε2+Oε2+M1010ε1+M1001ε2+Oε2X−L20+OεX2+M02+OεXY+OεY2+OX,Y3,and rewrite X, Y into μ and ν to get(44)μ˙=ν,ν˙=υ¯1ε+υ¯2εν+M20+Oε2μ2+2L20+M11+Oεμν+νOμ,ν2.
Now, we can rescale the time and variables by formula (35) and(45)υ¯1ε,υ¯2ε⟶M20+Oε232L20+M11+Oε4υ1ε,M20+Oε22L20+M11+Oευ2ε,such that system (44) becomes(46)μ˙=ν,ν˙=υ1ε+υ2εν+μ2+μν+νOμ,ν2.
Here, we omit theυε expansions because of the complexity and only verify the nondegeneracy conditions. After a series of calculations, we get(47)∂υ1,υ2∂ε1,ε2ε=0=−43P^I2−222P^I−15AB5τ3P^I32P^I2−2P^I+18C5,where(48)A=3ω−τP^I6+15τ−ωP^I5+10ω−13τP^I4+10ω+4τP^I3−8ω+9τP^I2+9τP^I−2τ,B=3ω−τP^I3+6τP^I2−2ω+5τP^I+2τ,C=3ω−τP^I3−3ω−τP^I2−2ω+τP^I+2ω.
Since6/3<P^I<1, it is sufficient to prove that ABC≠0. If A=0, we get(49)ω=τ3P^I6−15P^I5+13P^I4−4P^I3+9P^I2−9P^I+2P^I2P^I−1P^I−43P^I2−2<0,which leads to contradiction. Similarly, it can be proved that BC≠0. Therefore, the nondegeneracy property holds, and the proof is completed.


## 3.1. Backward Bifurcation

For the stability of the disease-free steady state (DEF), we have the following proposition.Proposition 1.
The DEF is asymptotically stable ifR0<1 and is unstable if R0>1.Proof.
It is easy to find that−γ is an eigenvalue of the Jacobian matrix at DEF. The other two eigenvalues λ1 and λ2 are the roots of(16)λ2+λ3γ−τn+ω+τ+2γγ−τn+ω=0.
IfR0<1, then τn<γ+ω, which leads to 3γ−τn+ω+τ>0. Thus, both λ1 and λ2 have negative real parts, and therefore, the DEF is asymptotically stable. If R0>1, then one root of (16) is positive, which means that the DEF is unstable.
There exist two endemic steady states whenR0<1. This implies the possibility that the system has a backward bifurcation. Let us now figure out the conditions for the backward bifurcation to appear. Set(17)a0=2γnγ+ω−β+γ+ω2n+1γ+ωn−1,b0=2γn22γn+γ+ω.
By the technique of central manifold [27], we obtain the following theorem.Theorem 1.
System (5) exhibits a backward bifurcation at DEF when R0 crosses unity if a0>0 and a forward bifurcation if a0<0.
The proof of Theorem 1 is deferred to Appendix A. It is noting that the coefficientβ has a large impact on the backward bifurcation although it has no effect on R0. From Table 1, we find that the system may undergo a saddle-node bifurcation when the two equilibria coalesce into one equilibrium. Thus, we derive the conditions for the saddle-node bifurcation. Interestingly, we show that the system has a forward bifurcation curve that looks like “S” shape.


## 3.2. Saddle-Node and Hopf Bifurcations

In this section, we first present the conditions that system (5) undergoes a saddle-node bifurcation. The Jacobian matrix of the system at an endemic state E∗=PI∗,PSI∗,PSS∗ is(18)JE∗=−γτn202γPI∗−2βPI∗2+3β−ω−2γPI∗−β+ω+γτn1−PI∗2−2γPI∗1−PI∗−τ−2γ2γPI∗1−PI∗−12γPI∗−2βPI∗−γ−ω+βτn1−PI∗0−2γPI∗1−PI∗.The eigenvalues ofJE∗ are the roots of(19)λ3+G2PI∗λ2+G1PI∗λ+G0PI∗=0,where(20)G2PI∗=−γPI∗−τPI∗+3γ+τ−1+PI∗=γPI∗+τ1−PI∗+3γ1−PI∗>0,(21)G1PI∗=γ2βPI∗3+−3β+ω−τPI∗2+β+3γ−ωPI∗+2γ+τ−1+PI∗2,(22)G0PI∗=2γ2PI∗F′PI∗1−PI∗2.Thus, the matrixJE∗ admits a zero eigenvalue if and only if G0PI∗=0, i.e., F′PI∗=0. Furthermore, zero is a simple eigenvalue if G1PI∗≠0, and Table 1 shows that E∗ is an endemic equilibrium of multiplicity 2 or 3 when Δ0=0. Consequently, we can state the following theorem.Theorem 2.
AssumeΔ0=0. Then, system (5) has an endemic equilibrium E^=P^I,P^SI,P^SS of multiplicity 2 if one of the following is satisfied:(a)
R0<1, C≥0, and ℬ<0(b)
R0<1 and C<0(c)
R0>1, C>0, ℬ<0, and A>0(d)
R0=1, C>0, and ℬ<0
Furthermore, system (5) undergoes a saddle-node bifurcation at E^ if G1PI^≠0.
The proof of Theorem 2 is omitted because it is similar to that in [19, 26, 28].
From Theorems1 and 2, we find that the rewiring rate β has remarkable impact on behaviors of system (5). To illustrate this, we depict the bifurcation diagrams of system (5) with different β in Figure 1, where the basal parameters are γ=0.1,n=6,ω=2.9. In Figure 1, when the rewiring rate β increases from zero, the type of bifurcation changes. For 0<β<2.3, the bifurcation at DEF is backward, and there exists a saddle-node bifurcation SN (see the solid bifurcation curve of Figure 1). At β=2.3, which leads the critical condition a0=0, the DEF becomes an equilibrium of multiplicity 3 (see the dashed bifurcation curve of Figure 1). As β increases further, the bifurcation at DEF becomes forward, and there exist two saddle-node bifurcations SN1 and SN2 for β>2.3. It is notable that the equilibrium SN2 coincides with DEF at β=2.3 if R0=1. In a word, we find that the coefficient β can induce the bifurcation curve which looks like “S” shape (see the dotted dash bifurcation curve and dotted bifurcation curve of Figure 1). By Theorems 1 and 2, when β>2.3, the bifurcation curve will remain S-shaped. More importantly, the disease can persist even though R0<1 when the bifurcation at DEF is forward (see the dotted dash bifurcation curve of Figure 1).
By (21), we know G2PI∗>0. Using the Routh–Hurwitz criterion and invoking the method to examine the Hopf bifurcation in [29], we see that the equilibrium point PI∗,PSI∗,PSS∗ is stable if(23)G0PI∗>0,G2PI∗G1PI∗−G0PI∗>0.
Moreover, the Hopf bifurcation occurs when(24)G2PI∗G1PI∗−G0PI∗=0,G1PI∗>0.Figure 1
Bifurcation diagrams of system (5) where SN, SN1, and SN2 represent saddle-node points. The parameter values are γ=0.1, n=6, ω=2.9, and β=0 (solid), 2.3 (dashed), 3.5 (dotted dash), or 15 (dotted).


## 3.3. Bogdanov–Takens Bifurcation

WhenΔ0=0 and one of the conditions (a), (b), (c), and (d) in Theorem 2 is satisfied, system (5) has an endemic equilibrium E^=P^I,P^SI,P^SS of multiplicity 2. It is easy to verify A≠0 under these conditions. Furthermore, we see that the components of E^ are given by(25)P^I=−B2A,P^SI=−γBτnA,P^SS=−βB24τnA2+γ+ω−β2τnAB+γ+ωτn.IfG1PI^≠0, by Theorem 2, E^ is a saddle node. If G1PI^=0, it is obvious that the eigenvalues of the Jacobian matrix at E^ are 0, 0, −G2PI^<0. We show in the following E^ is a Bogdanov–Takens point.Theorem 3.
Assume one of the conditions (a), (b), (c), and (d) in Theorem2 is satisfied. If Δ0=0 and G1PI^=0, then E^ is a Bogdanov–Takens point, and system (5) localized at E^ is topologically equivalent to(26)μ˙=ν,ν˙=μ2+μν+Oμ,ν3.Proof.
We just outline the main steps and leave the cumbersome details in Appendix B. First, let us bringE^ to the origin by substituting x=PI−P^I,y=PSI−P^SI,z=PSS−P^SS. Then, system (5) becomes(27)x˙=−γx+τn2y,y˙=2γP^IτnβP^I+γ+ω1−P^I−γP^I1−P^I2−βx+−2γP^I1−P^I−τ−2γy+2γP^I1−P^I−1z−τnP^SIP^SI−2P^SS21−P^I3x2−βP^I2+τnP^SI−τnP^SS−2βP^I+β1−P^I2xy+τnP^SI1−P^I2xz−τn21−P^Iy2+τn1−P^Iyz+Ox,y,z3,z˙=2γP^Iτn−βP^I+γ+ω1−P^I+βx−2γP^I1−P^Iz−τnP^SSP^SI1−P^I3x2+βP^I2−τnP^SS−2βP^I+β1−P^I2xy−τnP^SI1−P^I2xz−τn1−P^Iyz+Ox,y,z3.
BecauseE^ is an equilibrium of multiplicity of 2, we can get G0PI^=0 from (21). From (19), we see that the Jacobian matrix at E^=P^I,P^SI,P^SS has three eigenvalues: λ1=λ2=0 and λ3=−G2P^I<0 if G1P^I=0. The generalized eigenvectors corresponding to 0 are(28)V1=−τnγ+2P^I−1β+ω,−2γγ+2P^I−1β+ω,1′,V2=τn2βP^I−β+γ+ωγ,0,−P^I+12γP^I′,and the eigenvectors corresponding to λ3 are(29)V3=−τnτ+γP^I−τ−3γ2γP^I2βP^I−β+γ+ω,τ+γP^I−τ−3γτ−2γP^I−τ−2γγP^I2βP^I−β+γ+ωP^I−1,1′.
Here, the symbol′ means the transpose of a vector. Then, the matrix V=V1,V2,V3 is used for the nonsingular linear transformation:(30)xyz=Vμνξ.
As a result, we get(31)μ˙=ν+L20μ2+L11μν+L02ν2+ξ⋅Oμ,ν,ξ+Oμ,ν,ξ3,ν˙=M20μ2+M11μν+M02ν2+ξ⋅Oμ,ν,ξ+Oμ,ν,ξ3,ξ˙=−G2P^Iξ+Oμ,ν,ξ2.
The expressions of the coefficients are shown in Appendix B. By the center manifold theorem [30], system (31) restricted to the center manifold becomes(32)μ˙=ν+L20μ2+L11μν+L02ν2+Oμ,ν3,ν˙=M20μ2+M11μν+M02ν2+Oμ,ν3.
We use the near-identity transformation:(33)μ=X+12L11+M02X2+L02XY+OX,Y3,ν=Y−L20X2+M02XY+OX,Y3,and then rewrite X, Y into μ and ν to get(34)μ˙=ν,ν˙=M20μ2+M11+2L20μν+Oμ,ν3.
It is proved in Appendix B thatM20>0 and M11+2L20>0. If we rescale the time and variables by(35)μ,ν,t⟶M20M11+2L202μ,M202M11+2L203ν,M11+2L20M20t,then system (34) is topologically equivalent to normal form (26).
From Theorem3, we need two parameters to unfold the codimension two singularity of E^. Firstly, let τ,ω,n,γ1,β1 satisfy the conditions in Theorem 3. Then, we have(36)β1=−3ω+3τP^I2+2ω−τ2P^I−13P^I2−2,γ1=τP^I−123P^I2−2.
In the following, we consider system (5) for parameters β,γ in a neighborhood β1,γ1. Substituting(37)β=β1+ε1,γ=γ1+ε2,into system (5), we obtain a perturbed system(38)P˙I=τn2PSI−γ1+ε2PI,P˙SI=−γ1+ϵ2+τ+ωPSI+2γ1+ϵ21−PSI−PSS+τnPSS−PSI2PSI1−PI−β1+ε1PIPSI,P˙SS=γ1+ε2PSI−τnPSSPSI1−PI+ωPSI+β1+ε1PIPSI,where ε1,ε2 is in the neighborhood of 0,0.Theorem 4.
If parametersε=ε1,ε2 are sufficiently small and τ,ω,n,γ1,β1 satisfy the conditions in Theorem 3, then system (44) admits a generic unfolding of the cusp singularity of codimension 2.Proof.
We present only the main steps to illustrate the proof and omit the details of calculation because they are similar to those in [19, 26]. Shift E^ to the origin by substituting x=PI−P^I,y=PSI−P^SI,z=PSS−P^SS. Then, using the method in Appendix B, we get(39)n=−3ω+3τP^I4+6ωP^I3−ω+5τP^I2+4τ−4ωP^I+2ω−ττ2P^I−13P^I2−2.
Next, we apply transformation (30) to system (38) to obtain(40)μ˙=ν+l00010ε1+l00001ε2+L20μ2+L11μν+L02ν2+∑i+j+k+p+q≥2i+j≠2ifi+j+k+p+q=2lijkpqμiνjξkε1pε2q,ν˙=m00010ε1+m00001ε2+M20μ2+M11μν+M02ν2+∑i+j+k+p+q≥2i+j≠2ifi+j+k+p+q=2mijkpqμiνjξkε1pε2q,ξ˙=−G2P^Iξ+g00010ε1+g00001ε2+∑i+j+k+p+q≥2gijkpqμiνjξkε1pε2q.
According to the method in [19, 26], the local center manifold has the following form:(41)ξ=g00010ε1+g00001ε2G2P^I+∑i+j+p+q≥2hijpqμiνjε1pε2q.
Then, system (40) restricted to this center manifold is(42)μ˙=ν+l00010ε1+l00001ε2+L20μ2+L11μν+L02ν2+∑i+j+p+q≥2i+j≠2ifi+j+p+q=2lijpqμiνjε1pε2q,ν˙=m00010ε1+m00001ε2+M20μ2+M11μν+M02ν2+∑i+j+p+q≥2i+j≠2ifi+j+p+q=2mijpqμiνjε1pε2q.
We use the near-identity transformation(43)μ=1+L1010ε1+L1001ε2+Oε2X+L0010ε1+L0001ε2+Oε2+Oε2Y+12L11+M02+OεX2+L02+OεXY+OεY2+OX,Y3,ν=1+M0110ε1+M0101ε2+Oε2Y−l00010ε1+l00001ε2+Oε2+M1010ε1+M1001ε2+Oε2X−L20+OεX2+M02+OεXY+OεY2+OX,Y3,and rewrite X, Y into μ and ν to get(44)μ˙=ν,ν˙=υ¯1ε+υ¯2εν+M20+Oε2μ2+2L20+M11+Oεμν+νOμ,ν2.
Now, we can rescale the time and variables by formula (35) and(45)υ¯1ε,υ¯2ε⟶M20+Oε232L20+M11+Oε4υ1ε,M20+Oε22L20+M11+Oευ2ε,such that system (44) becomes(46)μ˙=ν,ν˙=υ1ε+υ2εν+μ2+μν+νOμ,ν2.
Here, we omit theυε expansions because of the complexity and only verify the nondegeneracy conditions. After a series of calculations, we get(47)∂υ1,υ2∂ε1,ε2ε=0=−43P^I2−222P^I−15AB5τ3P^I32P^I2−2P^I+18C5,where(48)A=3ω−τP^I6+15τ−ωP^I5+10ω−13τP^I4+10ω+4τP^I3−8ω+9τP^I2+9τP^I−2τ,B=3ω−τP^I3+6τP^I2−2ω+5τP^I+2τ,C=3ω−τP^I3−3ω−τP^I2−2ω+τP^I+2ω.
Since6/3<P^I<1, it is sufficient to prove that ABC≠0. If A=0, we get(49)ω=τ3P^I6−15P^I5+13P^I4−4P^I3+9P^I2−9P^I+2P^I2P^I−1P^I−43P^I2−2<0,which leads to contradiction. Similarly, it can be proved that BC≠0. Therefore, the nondegeneracy property holds, and the proof is completed.


## 4. Numerical Simulations

In this section, we use numerical software MATCONT to demonstrate how the dynamical behaviors of model (5) are influenced by the parameters [31–33]. First, we choose β and γ as bifurcation parameters to obtain the bifurcation diagrams in Figure 2. The plane is split into seven regions by the saddle-node bifurcation curve, the Hopf bifurcation curve, and the line of R0=1. The distribution of endemic equilibria and their stability are described in Table 2. Furthermore, the Hopf bifurcation curve meets the saddle-node bifurcation curve at a Bogdanov–Takens point (BT) and is divided into two parts by the degenerate Hopf bifurcation point (GH), where the upper part is subcritical Hopf and the lower part is supercritical.Figure 2
Bifurcation diagrams in the planeβ,γ. The parameter values are fixed as τ=0.5,n=6,ω=2.9.Table 2
Stability of the steady states in different regions in Figure2.
DEF
E
1
E
2
E
3
I
Stable
Does not exist
Does not exist
Does not exist
II
Stable
Stable
Unstable
Does not exist
III
Stable
Unstable
Unstable
Does not exist
IV
Unstable
Stable
Does not exist
Does not exist
V
Unstable
Stable
Unstable
Stable
VI
Unstable
Unstable
Unstable
Stable
VII
Unstable
Does not exist
Does not exist
Stable
E
1
,
E
2
,
and
E
3 are three endemic steady states arranged according to sizes.It is notable that the Bogdanov–Takens point implies the existence of a homoclinic orbit. Thus, one can select the parameters in region II in Figure2 such that a homoclinic orbit coexists with a stable endemic equilibrium. This is the case if β=7.89 with τ=0.5,n=6,ω=2.9,γ=0.13, as shown in Figure 3(a).Figure 3
Global dynamics of system (5), where τ=0.5,n=6,ω=2.9,andγ=0.13. (a) Coexistence of a homoclinic orbit with a stable equilibrium in region II, where β=7.89. (b) Coexistence of an unstable cycle with a stable equilibrium in region II, where β=7.91. (c) Coexistence of an unstable cycle with a stable cycle in region III, where β=7.99. (d) E1 and E2 are unstable, and orbits tend to the disease-free equilibrium in region III, where β=8.1.
(a)
(b)
(c)
(d)Interestingly, an unstable periodical solution emerges asβ increases due to the homoclinic orbit is broken, which is shown in Figure 3(b), where β=7.91. As β increases further, the equilibrium E1 becomes unstable, and a stable periodic solution appears due to the supercritical Hopf bifurcation when β=7.99 (see Figure 3(c)). Thus, the system has two limit cycles. The unstable limit cycle shrinks, and the stable limit cycle expands as the rise of β until these two cycles coalesce with each other at β=8.01. Finally, system (5) has no periodic solution, and its orbits except for the endemic equilibria with the stable manifolds of the saddle tend to the disease-free equilibrium when β>8.01, as shown in Figure 3(d).In order to show the influence of rewiring coefficientβ on PI, we choose β as a bifurcation parameter to obtain the bifurcation diagrams in Figure 4(a). As β increases gradually from 0, system (5) has only one equilibrium which is globally asymptotically stable. Therefore, the final size of PI settles down to the higher state. As β passes through β1 and stays in β1,βh, the final size of PI depends on the initial size of PI. If the initial size of PI is smaller, the final size of PI stays in the lower branch of the bifurcation curve. Otherwise, the final size approaches the upper stable branch. When β crosses βh, the lager equilibrium point loses stability due to the occurrence of supercritical Hopf bifurcation, and a stable periodical solution emerges. As β>β2, the system has only one equilibrium which is globally stable and has very low value. It is clear that the rewiring coefficient β tends to decrease the endemic level.Figure 4
(a) Bifurcation diagram of equilibria.SN1 and SN2 represent saddle-node bifurcations, Hopf stands for a Hopf bifurcation, and solid (dotted) curves represent stable (unstable) equilibrium states. The parameter values are τ=0.5,n=6,ω=2.9,andγ=0.05. (b) Coexistence of a stable cycle with a stable equilibrium. A stable cycle is shown by the red curve, and a stable equilibrium is represented by E3, where β=9.60357 is close to βh shown in (a).
(a)
(b)Finally, we compare model (5) with stochastic simulations. The random network is generated by using the Gillespie algorithm [34]. In stochastic simulations, the initial network degree distribution of each sample is basically the same, which obeys the Poisson distribution, but the specific structure is different. It is shown in Figure 5 that model (5) and the simulations fit well when β is suitably large to cause a low infection rate. However, when β is smaller than a critical value, the error of model (5) to stochastic simulations increases, although the trend of value PI in the model and simulation are the same. This is because the strong adaptive response from the higher coefficient β decreases the infection to low level so that the degree distribution of the network is approximated well by the Poisson distribution. Therefore, approximation (3) is very suitable for the moment closure.Figure 5
Comparison of the time evolution of prevalence in pairwise model (5) and simulation models, with the continuous curves which are obtained from system (5), and the dotted curves are obtained from average of 50 simulations. The parameter values are γ=0.1,τ=1.5,n=10,ω=2.9,andβ=0.05 (green), 0.1 (blue), or 0.15 (red). The initial value is PI=0.2,PSI=0.19,andPSS=0.76.


## 5. Discussion

In this paper, we have extended the network model in [16] by incorporating the adaptive rewiring rate, that is, the constant rewiring rate of the model in [16] is allowed to be the rate which is correlated with the infective number to describe the self-protection of individuals to reduce infection risk. Through rigorous mathematical analysis, we have obtained the conditions for the backward bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcation. These extend the mathematical analysis in [19]. Interestingly, the coefficient β of adaptive rewiring can induce the complicated behaviors of model (5). As shown in Figure 1, the model exhibits the bifurcation curve with “S” shape due to the adaptive rewiring rate, which leads to multiple equilibria and twice saddle-node bifurcations. More importantly, the disease can persist even though the basic reproduction number is less than unity when the bifurcation at the disease-free equilibrium is forward. Moreover, we find a scenario where a stable periodic solution coexists with a stable equilibrium, which is shown in Figure 4(b). Finally, by Figure 4(a), we see that increasing β tends to decrease the infection level, although β has no direct contribution to R0 from Section 2.From the previous theoretical analysis and simulations, we can infer that if the magnitude of rewiring is positively related to the scale of infection, it produces remarkable help for the epidemic control. Especially if the adaptive response is strong enough, the epidemic outbreak can be contained.We have investigated the qualitative behaviors of system (5). It is very interesting to apply the theoretical results to the problems with real data, as done in the paper [35]. It would also be important to consider the case where the rewiring rate is a saturated function of infected nodes and consider the time delay of adaptive rewiring. We leave these as future projects.



---

*Source: 7310347-2020-05-08.xml*