Complex Behaviors of Epidemic Model with Nonlinear Rewiring Rate

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.


Introduction
Mathematical models of epidemic propagations on networks have been proposed and studied by many papers (see [1][2][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][8][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][11][12][13][14][15]. is 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. is 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. ere have been a number of studies about epidemic propagation on the adaptive network in this direction (see [11][12][13] and the references cited therein). e following model, which originates from [16,17], is presented in [18]:  [17]. An infected node can infect any susceptible neighbor at rate τ and recover to susceptible state at rate c. e 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. ey 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].
us, social contacts tend to disconnect infectious nodes and turn to safe connections. ese 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) It is clear that model (2) extends model (1) by introducing the infection-dependent response coefficient β. is 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. e 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.
is paper is organized as follows. In Section 2, we simplify the model by the technique of moment closure and study the distribution of equilibrium points of the reduced model. e 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. e paper ends with discussions in Section 5.

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]: e 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. e second equality in (3) has the similar meaning.
Let N be the number of vertices of the network and n be the average degree of the graph.
to obtain in which, the identities of P S + P I � 1 and P SI + P SS + P II � 1 are used, and β � βN. In view of biological background, we 2 Complexity assume that the parameters τ, c, and n are positive constants, and the parameters ω and β are nonnegative constants. It is easy to verify the domain Ω � P I , P SI , P SS 0 ≤ P SI , 0 ≤ P SS , 0 ≤ P I < 1, P SI + P SS ≤ 1.
is a positive invariant set of system (5). In the rest of this paper, we study the dynamical behaviors of the system in Ω. e disease-free steady state (DEF) of system (5) is E 0 � (0, 0, 1). With the help of the computation method in [25], we get the basic reproduction number: e endemic steady states are determined by Obviously, from equations (8) and (10), we have By substituting these expressions into equation (9), we can deduce P I which satisfies the following equation: where A � β, Set Motivated by [26], we define the discriminant of the cubic function F(P I ) by Notice that A ≥ 0 if Δ 0 � 0. Similar to [26], we see that the positive roots of equation (12) are classified by the signs of Δ 0 , B, C, and D. en, we state the distribution of equilibria of equation (12) in three cases. C1: R 0 > 1: (1) If C ≤ 0, then F(P I ) has a unique positive root.
(2) Suppose that C > 0 and B ≥ 0. en, F(P I ) has a unique positive root. (1) Suppose that C < 0. en, F(P I ) has a unique positive root.
(2) Suppose that C � 0 and B < 0. en, F(P I ) has a unique positive root. 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.

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.

Backward Bifurcation.
For the stability of the diseasefree steady state (DEF), we have the following proposition.

Proposition 1.
e DEF is asymptotically stable if R 0 < 1 and is unstable if R 0 > 1.

Complexity
Proof. It is easy to find that − c is an eigenvalue of the Jacobian matrix at DEF. e other two eigenvalues λ 1 and λ 2 are the roots of If R 0 < 1, then τn < c + ω, which leads to 3c − τn + ω + τ > 0. us, both λ 1 and λ 2 have negative real parts, and therefore, the DEF is asymptotically stable. If R 0 > 1, then one root of (16) is positive, which means that the DEF is unstable.
ere exist two endemic steady states when R 0 < 1. is implies the possibility that the system has a backward bifurcation. Let us now figure out the conditions for the backward bifurcation to appear. Set By the technique of central manifold [27], we obtain the following theorem. 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 R 0 . From Table 1, we find that the system may undergo a saddle-node bifurcation when the two equilibria coalesce into one equilibrium. us, we derive the conditions for the saddlenode bifurcation. Interestingly, we show that the system has a forward bifurcation curve that looks like "S" shape.

Saddle-Node and Hopf Bifurcations.
In this section, we first present the conditions that system (5) undergoes a saddle-node bifurcation. e Jacobian matrix of the system at an endemic state E * � (P * I , P * SI , P * SS ) is Range of R 0 Other conditions Equilibria of system (5) One positive equilibrium Two positive equilibria 1,2 mean that there are a simple equilibrium and another equilibrium of multiplicity 2. One positive equilibrium 2 means that there is an equilibrium of multiplicity 2. One positive equilibrium 3 means that there is an equilibrium of multiplicity 3. e remaining cases represent a simple equilibrium.
4 Complexity e eigenvalues of J(E * ) are the roots of where us, the matrix J(E * ) admits a zero eigenvalue if and only if G 0 (P * I ) � 0, i.e., F ′ (P * I ) � 0. Furthermore, zero is a simple eigenvalue if G 1 (P * I ) ≠ 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. en, system (5) has an endemic equilibrium E � (P I , P SI , P SS ) of multiplicity 2 if one of the following is satisfied:

Furthermore, system (5) undergoes a saddle-node bifurcation at
The proof of Theorem 2 is omitted because it is similar to that in [19,26,28].
From Theorems 1 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 c � 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 a 0 � 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 SN 1 and SN 2 for β > 2.3. It is notable that the equilibrium SN 2 coincides with DEF at β � 2.3 if R 0 � 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 eorems 1 and 2, when β > 2.3, the bifurcation curve will remain Sshaped. More importantly, the disease can persist even though R 0 < 1 when the bifurcation at DEF is forward (see the dotted dash bifurcation curve of Figure 1).
By (21), we know G 2 (P * I ) > 0. Using the Routh-Hurwitz criterion and invoking the method to examine the Hopf bifurcation in [29], we see that the equilibrium point Moreover, the Hopf bifurcation occurs when

Bogdanov-Takens Bifurcation.
When Δ 0 � 0 and one of the conditions (a), (b), (c), and (d) in eorem 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 If G 1 (P I ) ≠ 0, by eorem 2, E is a saddle node. If G 1 (P I ) � 0, it is obvious that the eigenvalues of the Jacobian matrix at E are 0, 0, − G 2 (P I ) < 0. We show in the following E is a Bogdanov-Takens point.

then E is a Bogdanov-Takens point, and system (5) localized at E is topologically equivalent to
Proof. We just outline the main steps and leave the cumbersome details in Appendix B. First, let us bring E to the origin by substituting x � P I − P I , y � P SI − P SI , z � P SS − P SS . en, system (5) becomes Because E is an equilibrium of multiplicity of 2, we can get G 0 (P I ) � 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 � − G 2 (P I ) < 0 if G 1 (P I ) � 0. e generalized eigenvectors corresponding to 0 are and the eigenvectors corresponding to λ 3 are

Complexity
Here, the symbol ′ means the transpose of a vector. en, the matrix V � (V 1 , V 2 , V 3 ) is used for the nonsingular linear transformation: x As a result, we get e expressions of the coefficients are shown in Appendix B. By the center manifold theorem [30], system (31) restricted to the center manifold becomes We use the near-identity transformation: and then rewrite X, Y into μ and ] to get It is proved in Appendix B that M 20 > 0 and M 11 + 2L 20 > 0. If we rescale the time and variables by then system (34) is topologically equivalent to normal form (26). From eorem 3, we need two parameters to unfold the codimension two singularity of E. Firstly, let (τ, ω, n, c 1 , β 1 ) satisfy the conditions in eorem 3. en, we have In the following, we consider system (5) for parameters (β, c) in a neighborhood (β 1 , c 1 ). Substituting into system (5), we obtain a perturbed system where (ε 1 , ε 2 ) is in the neighborhood of (0, 0). 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 � P I − P I , y � P SI − P SI , z � P SS − P SS . en, using the method in Appendix B, we get Next, we apply transformation (30) to system (38) to obtain According to the method in [19,26], the local center manifold has the following form: en, system (40) restricted to this center manifold is We use the near-identity transformation and rewrite X, Y into μ and ] to get 8 Complexity Now, we can rescale the time and variables by formula (35) and such that system (44) becomes Here, we omit the υ(ε) expansions because of the complexity and only verify the nondegeneracy conditions. After a series of calculations, we get where Since � 6 √ /3 < P I < 1, it is sufficient to prove that ABC ≠ 0. If A � 0, we get which leads to contradiction. Similarly, it can be proved that BC ≠ 0. erefore, the nondegeneracy property holds, and the proof is completed.

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][32][33]. First, we choose β and c as bifurcation parameters to obtain the bifurcation diagrams in Figure 2. e plane is split into seven regions by the saddle-node bifurcation curve, the Hopf bifurcation curve, and the line of R 0 � 1. e 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.
It is notable that the Bogdanov-Takens point implies the existence of a homoclinic orbit.
us, one can select the parameters in region II in Figure 2 such that a homoclinic orbit coexists with a stable endemic equilibrium. is is the case if β � 7.89 with τ � 0.5, n � 6, ω � 2.9, c � 0.13, as shown in Figure 3(a).
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 E 1 becomes unstable, and a stable periodic solution appears due to the supercritical Hopf bifurcation Complexity 9 when β � 7.99 (see Figure 3(c)). us, the system has two limit cycles. e 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 P I , 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. erefore, the final size of P I settles down to the higher state. As β passes through β 1 and stays in (β 1 , β h ), the final size of P I depends on the initial size of P I . If the initial size of P I is smaller, the final size of P I 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. Finally, we compare model (5) with stochastic simulations. e 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 P I in the model and simulation are the same. is 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. erefore, approximation (3) is very suitable for the moment closure.

Discussion
In this paper, we have extended the network model in [16] by incorporating the adaptive rewiring rate, that is, the constant   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. e parameter values are c � 0.1, τ � 1.5, n � 10, ω � 2.9, and β � 0.05 (green), 0.1 (blue), or 0.15 (red). e initial value is P I � 0.2, P SI � 0.19, and P SS � 0.76.
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.
rough rigorous mathematical analysis, we have obtained the conditions for the backward bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation. ese 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 saddlenode 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 R 0 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.