Bifurcation analysis of a pair-wise epidemic model on adaptive networks

1.   Introduction

Modelling has become an important approach in understanding and control of infectious diseases over recent years. Classical infectious disease model has a long history which can trace back to the study of smallpox in $1760$ by Bernoulli ^[1]. Since then the classical compartmental epidemic models are used more and more to study the epidemic spreading in human population. The homogeneous mixing model in early stage assume that the population are well-mixing, meaning that each pair of individuals has an equal probability of contacting which is idealized obviously. With the rise of complex network, it is found that the complex networks can characterize the contact pattern of people more practical. And increasingly, more and more researchers realize that modeling the spreading of epidemic incorporating complex social network structures is more reasonable. Thus Pastor-Satorras et al. established and studied the susceptible-infective-susceptible(SIS) model based on the node dynamics on static homogeneous and scale-free networks (so called node-based epidemic model), respectively ^[2]. They proved that there has no the epidemic threshold in scale-free network when the parameters are proper. Furthermore, Pastor-Satorras and Vespignani studied the transmission threshold of epidemic on scale-free networks with the limited scale ^[3]. Then Boguna et al. gave the critical value of SIS epidemic model with correlations on complex networks ^[4], and they found that the epidemic threshold inversely proportional to the largest eigenvalue of the connectivity of the average number of contaction with the type of a node with connectivity $k$ linking to the other node with degree $k^{'}$.

Following this idea, a large body of investigations on network epidemiology emerged based on the modeling approach of node-based epidemic model. Mieghem et al. ^[5] introduced a node-level epidemic model by remoulding a traditional SIS model. Yang et al. ^[6,7,8] studied the spreading of virus by taking into account heterogeneous, the influence of patch forwarding and bi-virus competing, respectively. Zhang et al. ^[9] established an SIS model for epidemic spreading on semi-directed networks, which contains the outbreak threshold on undirected networks and directed networks as special cases. Yang et al. ^[10] incorporated the infection age into on epidemic model on complex networks and studied the effect of which on epidemic spreading. Peng et al. ^[11] considered the dual influences of dynamical contact networks' structures and heterogeneous spatial constraints. Sun et al. ^[12] investigated the effect of feedback regulation on vegetation patterns in semi-arid environments on homogeneous networks, and gave a conclusion and compare between the homogeneous mixing models and network-based epidemic models^[13]. However, the node-based epidemic models on complex networks neglects the correlation of connections, Keeling et al. ^[14,15,16] proposed the moment closure epidemic model by taking the different types of edges into the dynamic models (so called edge-based epidemic model). Wang et al. ^[17,18,19] put forward several different edge-based epidemic models investigating the effect of multiple routes of transmission, infectious force in latent period on random networks and two-stage contact process on complex networks, respectively. Yang et al. ^[20] proposed a new coumputational approach for the basic reproduction number of epidemic models on complex networks. Li et al. ^[21,22] studied the spread of several specific infectious diseases through modeling approach of epidemic models on homogeneous complex networks.

Furthermore, modeling epidemic on dynamical networks has attracted tremendous interest recently. With the study of infections disease, people find that adaptive behavior is frequently happening. For example, if the susceptible discover they have infected neighbors, the edges between the susceptibles and infected neighbors can be automatically interrupted to avoid being infected. The edge after disconnecting can sometimes rewire with other non-infected node or don't rewire sometimes. Then most of the researches either mainly focus on the studies of structures of complex networks^[23,24] or dynamics of node-based epidemic models on networks ^{[25,26,27,28]}. Moreover, Zhang et al. ^[28] established the epidemic models on complex networks to study the spreading of seasonal-influenza diseases, and they found that that preferential imitation can invalidate targeted subsidy policies epidemic spreading. Gross et al. ^[29,30]studied a SIS edge-based model in dynamical networks. It's worth mentioning that Zhang et al. ^[31] gave a rigorously proof that the adaptive behaviors can lead to rich dynamics like degenerate Hopf bifurcation, homoclinic bifurcation and Bogdanov-Takens bifurcation analytical, etc. All of these studies suggested that a simple rewiring rule for the network has a profound impart on the network structure and dynamics behavior. In addition, Shaw et al. studied susceptible-infective-recovered-susceptible(SIRS) model in dynamics network and employ Monte Carlo for stochastic simulation. In this paper, it is assumed that once susceptible node has lost the connection with an infected neighbour, it will re-established connections only with other noninfected node, and there exist threshold and different steady state solutions when rewiring rate is zero. When rewiring rate isn't zero, there have rich dynamics which contains bistable phenomenon and Hopf bifurcation, transcritical bifurcation and saddle-node bifurcation ^[32]. By contrast, Zanette and Risau-Gusman ^[33] assumed that the new connections are allowed connected with other nodes randomly whatever their states are. This sort of rewiring is appropriate for asymptomatic diseases. Besides that, Schwarzkopf et al. ^[34] considered a epidemic model with an discriminate rewiring mechanism, depending on the nodes' degree instead of the infection process. More recently, Rogers et al. study a SIRS pair-wise proxy model by combining random activation-deletion links with 'smart' rewiring, they develop a method to predict the effects of stochasticity in adaptive network ^[35]. Moreover, Juher et al. ^[36] studied the epidemic spreading with two extinction scenarios by bifurcation analysis, and they obtained an analytical condition for the occurrence of bistability based on^[29]. Zhou et al. revealed the phenomenon of epidemic reemergence in growing network by the dynamics of bistability and oscillation^[37].

Make a survey of above studies, the triples approximation formula used to close the moment closure models is based on the Poisson distribution, in which the average degree are often used to stand for the nodes' degrees. However, there are very large difference among degrees of nodes of real world complex networks since which are usually heterogeneous. Thus the triples approximation formula overlooked the degree of center node of triples which has dramatic effects on the accuracy of approximation and the analysis of system. By incorporating the center node's degree as a variable of triples approximation formula, Kiss et al. ^[38] proposed a more precise approximation formula which can consider the effect of heterogeneous of the complex networks. Above all, in this paper, we study how the topological structures of adaptive networks influences epidemic dynamics using the new triples approximation formula. We firstly close the model and reduce which to a three dimensional system by some reasonable assumptions. Owing to the more nonlinearity of the new formula, the analysis of the corresponding closed model is more difficult. By the qualitative theory and the stability theory of ordinary differential equations, the basic reproduction number is derived, and the rich dynamics of backward bifurcation, saddle-node bifurcation and Hopf bifurcation are analyzed.

In our paper, we develop a new triple approximation formula ^[38] and make further theoretical analysis on the basis of the previous researches. The new triple approximation formula can involve the dramatic effects of center nodes' degree of triples on the structures, which makes the approximation more accurate. We close the epidemic model on adaptive networks on the basis of new triple approximation formula, and the analytical conditions for rich dynamics like backward bifurcation and saddle-node bifurcation are given. The outline of this paper is as follows. In Section 2, we give the SIS epidemic model based on the adaptive networks, then an equivalent nondimensionalized three dimensional model is derived under the new triple approximation formula. In Section 3, the existence and stabilities of equilibria are studied, and the occurrence of backward bifurcation and saddle-node bifurcation are proven analytically. In Section 4, the normal form of Hopf bifurcation and the first Lyapunov coefficient are derived by a numerical example. What is more, numerical simulations (bifurcation diagrams and time series plots) holding for these results are also provided. Finally a brief conclusion is given in Section 5.

2.   The model

We consider an adaptive network with the constant number of nodes and connections are $N$ and $K$, respectively. Thus the average degree of the network can be denoted as $n = \frac{2K}{N}$. Since there are no isolated nodes in the network and the network is connected, the degree of every node is larger that $1$, i.e. $n > 1$.

The individuals (nodes) of the network are divided into two groups: the susceptibles (S), and the infected (I). In the model, transmission will occur from infected nodes to susceptible nodes along an edge connecting an infected node with a susceptible individual (so called SI link). For every SI link, the susceptible is infected by his or her infective neighbor with fixed probability of $\tau$ per unit time. In additional, each infected individual recovers to susceptible at rate $\gamma$ independently. Furthermore, since the susceptible individuals are likely to protect themselves, so the connections between susceptible and infected individuals are also being rewired adaptively during the spreading of the epidemic. Each SI link is rewired with probability $w$ to connect the susceptible node to another randomly selected susceptible one. Assume that there are no self-loops and multiple connections between nodes to form in the rewiring process.

Let $[S]$ and $[I]$ denote the total number of susceptible and infected individuals in the adaptive network $\mathcal{G}$. Moreover, $[AB]$ and $[ABC]$ represent the number of pairs of type $A-B$ (with one member in state $A$, and the other member in state $B$) and the number of triples of type $A-B-C$ (with one edge member in state $A$, the middle member in state $B$, and the other edge member in state $C$), respectively, where $A, B, C\in{S, I}$. Then we have the balanced relationships, $[S]+[I] = N$ and $[SS]+2[SI]+[II] = 2K$, with respect to the nodes and edges of the network, respectively. Hence the spreading process on an adaptive network can be described by the following $SIS$ pairwise epidemic model

One can verify that $\dot{[S]}+\dot{[I]} = 0$ and $\dot{[SS]}+2\dot{[SI]}+\dot{[II]} = 0$, which corresponds to the balance conditions of nodes and edges respectively. However, these equations (2.1) are unclosed and cannot be analyzed. The proceeds of the computation and analysis need to seek approximations that allow us to close these equations at a lower dimension. That is to say, in order to close the equations (2.1) it is necessary to approximate the triples $[SSI]$ and $[ISI]$ that appear. Keeling et. al. ^[15] proposed an approximate formula in terms of pairs and singles

where $n$ is the average degree of the network, and $A, B, C\in\{S, I\}$. However, since all triples are formed from two pairs sharing a common central node, the approximation of triples depends on the degree of central node. Therefore Kiss et. al. ^[38] modified the equation (2.2) by replacing the average degree of the whole network with the group of central individuals

where $k_{B}$ denotes the expected average degree of central nodes with type $B$. In fact, all types of triples $[SSI]$ and $[ISI]$ are susceptible, therefore the average degree of susceptible individuals is only needed, which given by

Henceforth, we use equation (2.3) to approximate the triple structures in the network. Furthermore, utilizing the balance conditions of constant number of nodes and connections, the system (2.1) is closed at the level of pairs and singles

in which the derivations are taken with respect to $t'$. Therefore for biological feasibility, we just consider the system (2.5) in the following domain

It is not difficult to verify that all of the solutions in the first octant will stay or enter inside the set $\Omega_0$.

Moreover, for simplicity of expression and analysis, we will rescale the model (2.5) into a dimensionless system. Let

therefore system becomes into the following nondimensionalized model

where the dots represent the derivatives with respect to time $t$.

Since the average degree $n > 1$ holds, and according to (2.7), then one can define the parameter space as following

Further, due to $\Omega_0$ is positive invariant of system (2.5), then we have

is positive invariant of system (2.8).

In the rest of the paper, we will study the dynamics and bifurcations of the dimensionless system (2.8) in $\Omega$ with parameters in $\Lambda$.

3.   Dynamics analysis of the model

3.1.   Existence of equilibria

We find equilibria by setting the right hand side of system (2.8) equal to zero. One trivial equilibrium point is the disease-free equilibrium $(DFE)$, which is denoted by $E_{0}(0, 0, 1)$. The stability of the $E_0$ is determined by the Jacobian,

Evaluating the $J$ on the disease-free equilibrium ${\rm E_0}$ yields the characteristic polynomial,

Then we define the reproduction number as

Then, ${\rm E_0}$ is stable (unstable) for $R_0 \! < \! 1$ ($R_0 \! > \! 1$), as expected.

Next, we discuss the existence of endemic equilibrium. As mentioned above, the solution according to the endemic equilibrium, $x > 0$, is founded by letting the right side of the equations of the system (2.5) to zero. Then the coordinates of any endemic equilibrium $E(x, y, z)$ satisfy

where

it should be noted that $\Delta_{1} > 0$ always holds since whose discriminant is $-16A_1(A_2+1) < 0$, and we discard the other branch root of $z(x) = \frac{(A_{1}+A_{2}+1)(1-x)-x-\sqrt{\Delta_{1}}}{2A_{1}n}$ since it is negative. Define

and

where

Then the coordinate $x$ is a root of the following cubic polynomial equation

By a direct calculation, we analyze the roots of following equivalent equation $f(x)$ of (3.5) within the interval $(0, 1)$

where

In order to obtain the exact solutions of equation (3.6) by utilizing the ShengJin's formula ^[39], then we have

Because of the uncertainty of the sign of $\Delta$, we list all the possible real roots formulas below:

Cases 1. $\Delta < 0$ (implies $A > 0$), equation (3.6) has $3$ different real roots

where $\theta = \arccos\frac{-y_{N}}{h}$, $\alpha = \arccos\frac{y_{N}}{h}$ and $x_{11} > x_{12} > x_{13}$.

Cases 2. $\Delta = 0$, equation (3.6) has $3$ real roots,

(ⅰ) If $A > 0$, two of which are repeated roots

(ⅱ) If $A = 0$, there are $3$ equal roots at $x_{1, 2, 3} = x_N$.

Cases 3. $\Delta > 0$, the unique real root of equation (3.6) is

where $Y_{1, 2} = \frac{-y_{N}\pm\sqrt{\Delta}}{2 b_{3}}$.

We denote the possible endemic equilibria by $E_{ij}\left(x_{ij}, y(x_{ij}), z(x_{ij})\right)$ ($i, j = 1, 2, 3$ and $i\geqslant j$).

For $(A_{1}, A_{2}, n)\in\Lambda$, $f(1) = A_{1}(n-1)+A_{2} > 0$, $f''(1) = 2A_{2}(A_{1}(2n-1)+A_{2}) > 0$, we have the following three cases about the existence of endemic equilibria.

(1) $R_0 > 1$: In this case $f(0) = b_{0} < 0$. It is not difficult to derive that $b_{2} > A_{2}(A_{1}+2) > 0$, therefore $f''(0) = 2b_{2} > 0$, then $f''(x) > 0$ always holds when $x\in(0, 1)$, that is to say, $f'(x)$ is monotone increasing in the interval $(0, 1)$. We discuss the roots of $f(x)$ in the following situations:

(a) If $b_{1}\geqslant0$, i.e. $f'(0)\geqslant0$, thus $f'(x) > 0$ always holds when $x\in(0, 1)$, implying that $f(x)$ is monotone increasing in the interval $(0, 1)$.

(b) If $b_{1} < 0$, i.e. $f'(0) < 0$, on account of

which guarantees $x_{1}^{*}\in(0, 1)$. Therefore it is not difficult to derive that $f(x)$ is decreasing and increasing monotonically in the interval $(0, x_{1}^{*})$ and $(x_{1}^{*}, 1)$, respectively.

Hence, system (2.8) has a unique endemic equilibrium whose $x$-coordinate is the largest real root of $f(x)$.

(ⅰ) If $\Delta < 0$, the unique endemic equilibrium is $E_{11}\left(x_{11}, y(x_{11}), z(x_{11})\right)$ since $x_{11} > x_{12} > x_{13}$;

(ⅱ) If $\Delta = 0$, equation (3.6) has $3$ real roots, since $b_2$ and $b_0$ have opposite sign, then $A = 0$ and $\Delta = 0$ cannot be simultaneously true, namely it is nonexistent of $3$ equal roots. Hence, we just consider the case of $A > 0$, i.e. there has $3$ real roots, $x_{21}$, $x_{22} = x_{23}$, two of which are repeated roots. According to the Vieta theorem

then we can obtain $x_{21} > 0$, $x_{22} = x_{23} < 0$, therefore the unique endemic equilibrium is $E_{21}\left(x_{21}, y(x_{21}), z(x_{21})\right)$.

(ⅲ) If $\Delta > 0$, it is obviously that the unique endemic equilibrium is $E_{31}\left(x_{31}, y(x_{31}), z(x_{31})\right)$.

(2) $R_0 = 1$: In this case $f(0) = b_{0} = 0$, $b_{2} = A_{2}(A_{1}+2) > 0$ and $f(x) = x(b_{3}x^2+b_{2}x+b_{1})$, recall that $f(1) > 0$, therefore system (2.8) has a unique endemic equilibrium if and only if $b_{1} < 0$. Actually it is not difficult to prove that $\Delta < 0$ holds in this case, and the $3$ different real roots are

respectively. Thus the unique endemic equilibrium is $E_{11}\left(x_{11}, y(x_{11}), z(x_{11})\right)$.

(3) $R_0 < 1$: In this case $f(0) = b_{0} > 0$, thus $f(x) = 0$ has at least one negative real root, namely system (2.8) has up to two endemic equilibria. It is obvious that system (2.8) has no endemic equilibrium if $b_{1}\geqslant0$ and $b_{2}\geqslant0$ due to all of the coefficients of equation (3.6) are non-negative in this case. Therefore we just need to consider the existence of positive roots of $f(x) = 0$ in $(0, 1)$ when $b_{1} < 0$ or $b_{2} < 0$. It is not difficult to derive that $x_{1}^\ast > 0$.

(a) $\Delta < 0$ implies $A > 0$, $f(x) = 0$ has 3 different roots which satisfy $x_{13} < x_{12} < x_{11}$, then $x_{13} < 0$. Since $f(x)$ is monotonically decrease and increase at the interval $(x_{2}^\ast, x_{1}^\ast)$ and $(x_{1}^\ast, +\infty)$ respectively, thus $x_{12} < x_{1}^\ast < x_{11}$. According to the Vieta theorem

then $x_{11}$ and $x_{12}$ have the same sign, we have $0 < x_{12} < x_{11}$. On the other hand, $1\in[x_{12}, x_{11}]$ is false due to $f(1) > 0$, that is to say either $1 < x_{12} < x_{1}^\ast < x_{11}$ or $x_{12} < x_{1}^\ast < x_{11} < 1$, then it is not difficult to derive that the latter case holds if and only if $x_{1}^\ast < 1$ which is equivalent to $f'(1) > 0$ since $f''(1) > 0$. Above all, $f(x) = 0$ has 2 different roots satisfying $0 < x_{12} < x_{11} < 1$ if and only if $x_{1}^\ast < 1$ (i.e. $f'(1) > 0$), otherwise there are no roots in the interval $(0, 1)$. System (2.8) has 2 endemic equilibria $E_{11}\left(x_{11}, y(x_{11}), z(x_{11})\right)$ and $E_{12}\left(x_{12}, y(x_{12}), z(x_{12})\right)$.

(b) $\Delta = 0$ implies $A > 0$ (otherwise one can get a contradiction $b_{1} > 0$ and $b_{2} > 0$). $f(x) = 0$ has 3 roots $x_{21} < 0$ and $x_{22} = x_{23} = x_{1}^\ast > 0$. To sum up, $f(x) = 0$ has 2 roots satisfying $0 < x_{22} = x_{23} < 1$ if and only if $x_{1}^\ast < 1$ (i.e. $f'(1) > 0$), otherwise there are no roots in the interval $(0, 1)$. System (2.8) has 2 identical endemic equilibria $E_{22}\left(x_{22}, y(x_{22}), z(x_{22})\right)$ which $x_{22} = x_{23} = x_{1}^\ast > 0$.

(c) $\Delta > 0$. It is obviously that $f(x) = 0$ has a unique real root $x_{31} < 0$, there are no roots in the interval $(0, 1)$.

Note that if $b_{1} < 0$ or $b_{2} < 0$ and $A > 0$, one have $f(x_{2}^\ast)\geqslant f(0) > 0$ since $x_{1}^\ast > 0$ and $f(x)$ is monotonically increase and decrease at the interval $(-\infty, x_{2}^\ast)$ and $(x_{2}^\ast, x_{1}^\ast)$ respectively. Therefore the sign of $\Delta$ is completely consistent with $f(x_{1}^\ast)$ because of $\Delta = f(x_{2}^\ast)f(x_{1}^\ast)$.

We summarize the result above in the following theorem.

Theorem 3.1. For system (2.8) with $(A_{1}, A_{2}, n)\in\Lambda$ and $(x, y, z)\in\Omega$, $f(1) = A_{1}(n-1)+A_{2} > 0$, $f''(1) = 2A_{2}(A_{1}(2n-1)+A_{2}) > 0$,

$(1)$ The disease-free equilibrium $E_{0}$ always exists.

${\rm (2)}$ When $R_0 > 1$, there exists a unique endemic equilibrium which is (a) $E_{11}$ when $\Delta < 0$; (b) $E_{21}$ when $\Delta = 0$; (c) $E_{31}$ when $\Delta > 0$.

${\rm (3)}$ When $R_0 = 1$, there exists a unique endemic equilibrium $E_{11}$ if and only if $b_{1} < 0$.

${\rm (4)}$ When $R_0 < 1$, and

${\rm (a)}$ If $b_{1}\geqslant0$ and $b_{2}\geqslant0$ or $f'(1)\leqslant0$, there is no endemic equilibrium.

${\rm (b)}$ If $b_{1} < 0$ or $b_{2} < 0$ and $f'(1) > 0$, system (2.8) has two endemic equilibria $E_{11}$ and $E_{12}$ if and only if $\Delta < 0$(implies $A > 0$), and these two equilibria coalesce into one endemic equilibrium $E_{11} = E_{12} = E_{22}$ if and only if $\Delta = 0$(implies $A > 0$). Otherwise, system (2.8) has no endemic equilibrium.

3.2.   Stability of disease-free equilibrium

We primarily study the stability of disease-free equilibrium $E_0$ in this subsection. By considering the distribution of the roots of the characteristic equation (3.2) of Jacobian at disease-free equilibrium $E_0$, one can obtain ${\rm E_0}$ is stable (unstable) for $A_1+A_2+1 > A_1n$ ($A_1+A_2+1 < A_1n$). Then the disease-free equilibrium $E_0$ is locally asymptotically stable (resp. unstable) for $R_{0} < 1$ (resp. $R_{0} > 1$).

We determine the types of disease-free equilibrium $E_0(1, 0, 1)$ in the following.

Theorem 3.2. For the system (2.8), the disease free equilibrium $E_0(1, 0, 1)$ is

${\rm (a)}$ $R_0 < 1$: an attracting node;

${\rm (b)}$ $R_0 > 1$: an hyperbolic saddle;

${\rm (c)}$ $R_0 = 1$ and

(ⅰ) $A_2(A_1+A_2) > 1$: a saddle-node of codimension 1;

(ⅱ) $A_2(A_1+A_2) < 1$: a saddle-node of codimension 1;

(ⅲ) $A_2(A_1+A_2) = 1$: an attracting semi-hyperbolic node of codimension 2.

Proof. For the system (2.8), form the characteristic equation (3.2), it is not difficult to obtain that $E_0$ is an attracting node if $R_0 < 1$, and $E_0$ is an hyperbolic saddle if $R_0 > 1$.

If $R_0 = 1$, the eigenvalues of Jacobian at $E_{0}$ are $-1$, $-(A_1+2)$ and $0$. To determine the type of $E_{0}$, we linearize the system (2.8) at $E_{0}$ and diagonalize the linear part, then we obtain the following system

Therefore if $A_2(A_1+A_2)\neq1$, it is unnecessary to calculate the center manifold, and $E_{0}$ is a saddle-node of codimension 1.

If $A_2(A_1+A_2) = 1$, according to the center manifold theorem, there exists a center manifold for system (3.7), which have the following form

then one can obtain

so $E_{0}$ is an attracting semi-hyperbolic node of codimension 2.

4.   The analysis of bifurcation and numerical simulations

Owing to the complexity of the expressions of the equilibria, in this section we launch into a detailed discussion of the effect of adaptive behavior on the epidemic spreading through numerical simulations, including stabilities, backward bifurcation, saddle-node bifurcation and Hopf bifurcation, etc. In order to study the effect of topology of adaptive networks aroused by rewiring rate $w$ and average degree $n$, infection rate $\tau$ and recovery rate $\gamma$ on disease transmission respectively, we perform the numerical simulation of the system (2.9) based on the original parameter instead of $A_1$ and $A_2$. From the analysis (2.7) and (3.3) in Sections 2 and 3, the basic reproduction number can be represented as

First, we simulate the time series of total densities of infected individuals on adaptive networks based on the different values of rewiring rate $w$ by fixing the other three parameters $\tau = 0.008$, $\gamma = 0.002$ and $n = 10$. In Figure 1(a), (d) and (g), we choose $w = 0.004$, which corresponding to $A_1 = 4$ and $A_2 = 2$, in this case the basic reproduction number $R_0 = 5.714$; In Figure 1(b), (e) and (h), we choose $w = 0.04$, which corresponding to $A_1 = 4$ and $A_2 = 20$, in this case the basic reproduction number $R_0 = 1.600$; In Figure 1(c), (f) and (i), we choose $w = 0.4$, which corresponding to $A_1 = 4$ and $A_2 = 200$, in this case the basic reproduction number $R_0 = 0.195$; We can see that if $R_0 < 1$, the disease will disappear (see Figure 1(c)), otherwise $R_0 > 1$ the endemic equilibrium will stable (see Figures 1(a), (b)). Furthermore from Figures 1(a) and (b), we can see that the infected individuals tends to the same level even though the different initial values $x(0) = 0.96$ (red line), $x(0) = 0.55$ (green line) and $x(0) = 0.20$ (blue line) are selected, from which the infected individuals tends to zero, namely the disease extinction is determined by the basic reproduction number $R_0$ in this case. In the other hand, the disease will persist or extinct when the rewiring rate $w$ is lower or higher, respectively.

Next we show the backward bifurcation, saddle-node bifurcation and Hopf bifurcation can really happen by plotting the relations of infected individuals versus parameters $\tau$ and $w$ respectively, as shown in Figure 2. One can see that if $R_0 > 1$, there is an unique endemic equilibrium, however $2$ endemic equilibria appear at $R_0 = 1$ and disappear when coalesce into one endemic equilibrium with multiplicity $2$, namely backward bifurcation occurs when $R_0 = 1$ and disappears at the critical point $LP$ of saddle-node bifurcation. Since $R_0$ is monotonically decrease (or increase) with the increasing of $w$ (or $\tau$), the backward bifurcation curves are opening left and right in Figures 2(a) and 2(b) respectively. Comparing to Figure 2(a) and 2(b), it is different that there is a unique endemic equilibrium when the rewiring rate $w$ is small or the infection rate $\tau$ is large. In addition $2$ endemic equilibria appear at $R_0 = 1$ in common, the larger one is stable initially (red solid line)and the lower one is always a saddle node(blue dashed line).

We continue to explore the the bistability and oscillations of the system. In Figure 2(a), we can also see that there have a critical point HP across which the stability of endemic equilibrium with larger epidemicity is changed, then we can observe Hopf bifurcation bifurcating from critical point HP and disappears when $w$ is large (Figure 2(a)) or $\tau$ is small (Figure 2(b)). Since the disease-free equilibrium $E_0$ is stable when $R_0 < 1$, then the bistability occurs from the $R_0 = 1$ and disappears at the critical point HP, namely both the larger endemic equilibrium and disease-free equilibrium are stable. Furthermore, limit cycle will arise near the critical point HP, i.e. the periodic oscillations can occurs for the system (2.8). To conform above conjectures, we fix two parameters $\gamma = 0.002$ and $\tau = 0.008$, and choose $w = 0.3044$ close to the HP point in Figure 2(a), the phase portrait of system is shown in Figure 3. In this case the disease-free equilibrium $E_0$ is stable node and the larger endemic equilibrium is an unstable focus. Both the orbits of blue and green curve with initial conditions $P_1(0.8, 0.04, 0.15)$ and $P_3(0.7, 0.12, 0.15)$ spiral outward (or inward) to a same stable limit cycle. The orbit (red curve) with initial condition $P_2(0.85, 0.05, 0.10)$ spirals outward to the disease-free equilibrium $E_0$. There this example have proved that the system has a stable limit cycle, see Figure 3.

4.1.   The analysis of bifurcation

Furthermore we show the Hopf bifurcation will really happen through a concrete example in this Subsection. In additional, the type of Hopf bifurcation(supercritical or subcritical) is also determined by calculating the first Liapunov coefficient. By selecting the parameters $\gamma = 0.002$, $\tau = 0.008$ and $n = 10$, by a direct calculation one can obtain the critical value of rewiring rate at HP point is $w = 0.3008852$. In this case the larger endemic equilibrium is $E_1(0.9687103, 0.04843551, 0.1210720)$.

In order to derive the norm form of Hopf bifurcation, we draw into the following affine transformation,

where the columns of the nonsingular transformation matrix are the eigenvectors of Jacobian matrix of system (2.8) at $E_1$ with selected parameters, corresponding to the eigenvectors of $\lambda_{1, 2} = \pm0.0404113i$ and $\lambda_{3} = -0.2576751$, respectively. We obtain the following normal form dynamical system,

where

A direct calculation to obtain the center manifold which is

then the system restricted to the center manifold is given by

where

Thus from the Ref ^[40,41] we can calculate the first Liapunov coefficient of system $L_1$ which is

Therefore supercritical Hopf bifurcation occurs near the critical point HB which is consistent with the numerical simulation in Figure 3.

5.   Conclusion

For the purpose of assessing the impact of adaptive behavior on the epidemic spreading, a new triples approximation formula incorporating the effect of center nodes' degree of triples is proposed, by which the SIS epidemic model on adaptive complex networks is closed. We first transform the model into an equivalent nondimensionalized three dimensional system, and the existence of equilibria is derived. Then the basic reproduction number $R_{0}$ is obtained through the analysis of the stabilities of disease-free equilibrium $E_0$. In addition, the conditions of the existence of endemic equilibria are analyzed strictly. By the qualitative theory and the stability theory of ordinary differential equations, the local stabilities of the endemic equilibria are proved. Moreover, using the methods of bifurcation theory, we show that the model exhibits rich and complicated dynamics such as transcritical forward bifurcation, backward bifurcation, saddle-node bifurcation. Besides, owing to the complexity, we calculate the normal form of Hopf bifurcation numerically, and the first Lyapunov coefficient is derived to determine the type of Hopf bifurcation (supercritical or subcritical). Then we show that a stable limit cycle can bifurcate from the endemic equilibrium with larger epidemicity, and the numerical simulations which is consistent with the analytical results are given.

The topological structure of adaptive network has obvious significance on epidemic spreading in theory, i.e. adaptive behavior can lead to rich dynamics. At the point of the application, on account of that backward bifurcation implies that the the basic reproduction number $R_0 < 1$ is not enough to control the spread of infectious diseases any more, and Hopf bifurcation means the oscillations will occur, therefore adaptive behavior raises the difficulty of prevention and control of infectious diseases. Keeling et al. ^[14,16] studied the spread of the childhood epidemics and sexually transmitted diseases by moment closure epidemic model without adaptive behaviors, and they show that the models have no rich dynamics. Therefore, this paper show that adaptive topology is an important reason of periodic outbreaks of epidemic as well as the nonlinear incidence rate and recovery rate in well-mixing epidemic models, which expands people's understanding on the epidemic oscillation in theory. And in application, by controlling the parameters in a 'safe' region to break the occurrence the bifurcations, the research on this paper can provide a new guidance to prevention and control of diseases.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant Nos.11671241, 11601294, 61803242, 11501340, 11801340, the Youth Science Fund of Shanxi Province under Grant Nos. 201801D221001, 201801D221011, 201601D021012, 201801D221007 and Shanxi Scholarship Council of China under Grant No. 2016-011.

Conflict of interest

All authors declare no conflicts of interest in this paper.