key: cord-010739-28qfmj9x
authors: Sherborne, N.; Blyuss, K. B.; Kiss, I. Z.
title: Bursting endemic bubbles in an adaptive network
date: 2018-04-09
journal: nan
DOI: 10.1103/physreve.97.042306
sha: 
doc_id: 10739
cord_uid: 28qfmj9x

The spread of an infectious disease is known to change people's behavior, which in turn affects the spread of disease. Adaptive network models that account for both epidemic and behavioral change have found oscillations, but in an extremely narrow region of the parameter space, which contrasts with intuition and available data. In this paper we propose a simple susceptible-infected-susceptible epidemic model on an adaptive network with time-delayed rewiring, and show that oscillatory solutions are now present in a wide region of the parameter space. Altering the transmission or rewiring rates reveals the presence of an endemic bubble—an enclosed region of the parameter space where oscillations are observed.

The spread of an infectious disease changes the behavior of individuals, and this, in turn, affects the spread of the disease [1] . Broadly speaking, responses to an epidemic fall into two categories: coordinated and uncoordinated. Coordinated responses include vaccination and quarantine schemes, travel restrictions, and information spread through mass media. Uncoordinated responses cover individuals adapting their behavior based on their own perceived risk; this includes improved hygiene regimens and avoiding crowded places and public transport during outbreaks. Surveys consistently identify such precautionary measures taken by individuals during epidemic outbreaks [2, 3] . Fear of becoming infected during the 2003 SARS epidemic in Hong Kong caused huge behavioral shifts; air travel into Hong Kong dropped by as much as 80% [4] . Responses to a large study covering numerous European and Asian regions revealed that, in the event of an influenza pandemic, 75% of people would avoid public transport, and 20-30% would try to stay indoors [5] . These behavioral shifts change the potential routes for transmission and can alter the size and time scale of an epidemic [6] .

In the context of epidemic models on networks, perhaps, the most widespread approach to couple epidemics and behavior is by using adaptive networks, where behavioral changes are captured by link rewiring based on the disease status of nodes [6, 7] . Gross et al. [8, 9] considered a simple susceptible-infectedsusceptible (SIS) model with rewiring, in which susceptible nodes disconnect from infected neighbors at rate ω, and immediately reconnect to a randomly chosen susceptible node. This simple model led to bistability and to oscillatory solutions, albeit with oscillations limited to an extremely narrow region of the parameter space. This rewiring procedure has since been extended to consider scenarios where both the susceptible and * k.blyuss@sussex.ac.uk infected nodes can rewire, and diseases with a latent period [10] . Zhang et al. [11] presented a further alternative, where news about past prevalence influences whether nodes choose to disconnect edges. The authors found an estimate of the critical delay that induces a Hopf bifurcation, thus causing periodicity. Tunc et al. [12] studied a network model with temporary deactivation of edges between susceptible and infected individuals. On a growing network, Zhou et al. [13] showed that cutting links between susceptible and infected individuals can lead to epidemic reemergence, with long periods of low disease prevalence punctuated by large outbreaks.

Periodic cycles and disease reemergence are evident in realworld data. Many diseases are subject to seasonal peaks, which have been studied extensively [14, 15] . Often a sinusoidal or other form of time-varying transmission parameter is used to imitate seasonality, which can lead to multiennial peaks [16] . A number of models have identified other possible causes of periodicity in epidemic dynamics. To give one example, Hethcote et al. [17] showed that in a well-mixed population temporary immunity in susceptible-infectedrecovered-susceptible (SIRS)-or susceptible-exposedinfected-recovered-susceptible (SEIRS)-type models as represented by a time delay can result in the emergence of periodic solutions when the immunity period exceeds some critical value.

One should note that seasonality alone cannot explain all cases of oscillations. In both the United Kingdom and the United States, the 2009 H1N1 pandemic occurred in two distinct waves separated by a few months [18, 19] . Other diseases have shown more long-term trends. Incidence reports of mycoplasma pneumonia have found evidence of epidemic cycles in many different countries, with periodicity of three to five years [20, 21] . Recently, it has been suggested that syphilis exhibits periodic cycling [22] , although these findings have been subsequently questioned [23] . While it is difficult to pinpoint the specific causes of periodicity in the dynamics of these diseases, if syphilis epidemics are indeed cyclical, then changes in human behavior have been proposed as the likely explanation [24] .

Intuitively, and as shown by empirical observations, one would expect oscillations to appear in epidemic models where behavior is considered. If an individual is aware of the state of their neighbors and responds accordingly, then times of high prevalence will be associated with greater caution, curbing further spread. Conversely, without advance warning, behavior will return to normal as prevalence wanes, enabling a second wave of the epidemic. Despite this intuition, adaptive network models have so far not been able to show such robust oscillations over reasonable regions of the parameter space. To tackle this problem, we introduce a simple SIS model on an adaptive network with N nodes. Infected nodes transmit the disease to susceptible neighbors at rate β across links, and recover and become susceptible again at rate γ , independently of the network. Susceptible nodes cut links that connect them to infected neighbours at rate ω and, after a fixed time delay of length τ , reconnect to susceptible nodes chosen uniformly at random from all such available nodes. The delay between cutting and reconnecting is crucial. It is unrealistic to expect that alternative contacts can be identified and established arbitrarily quickly. The delay represents both people's hesitance to make new contacts and also the potential lack of availability of such new contacts when an epidemic is spreading thorough a population [5] .

To construct the mean-field model, we use the pairwise approximation method [25] . The number of nodes in the susceptible or infected state at time t is denoted by [S] and [I ], respectively; [SS], [SI ], and [I I ] denote the number of connected pairs of nodes in the respective states, with all pairs being doubly counted. The explicit dependence on time is dropped for simplicity. For the moment closure approximation we use the assumption that once a node is fixed, typically a susceptible node, then the states of the neighbors are Poisson distributed [26] . This leads to

to express the number of connected triples [8, 25] . The delay before an S-I edge is rewired to an S-S edge introduces a complication, as not all newly formed edges will be between two susceptible nodes. To see this, consider an example of a susceptible node with two or more infected neighbors. At some time t 1 it disconnects from one of these neighbors. Then, in the interval (t 1 ,t 1 + τ ) another infected neighbor transmits the disease to it. If it then remains infected until time t 1 + τ , the new edge will be of an I -S type rather than S-S. To deal with this issue we use a technique similar to that used by Kiss et al. [27] for a pairwise model with an infectious period of fixed length. Consider y p (t) to be the cohort of susceptible nodes that have cut a link at time t − τ and are waiting to reconnect. The expected number of infected neighbors a susceptible node has is approximated by [SI ]/[S]. Therefore, the rate at which nodes in the cohort become infected over the interval (t − τ,t) iṡ

The solution to this ordinary differential equation is

A member of the cohort infected at some time u ∈ (t − τ,t) may recover before time t. To ensure that we only consider nodes which remain infected, we must include the probability that a node infected at time u remains infected until time t in the integral term of (2). This is the survival probability of the recovery process, and it is given by e −γ (t−u) . Therefore, the rate at which new S-S edges are formed is

If the exponential term in (3) is denoted by x(t), the rate at which new I -S edges are formed is

With this in mind, the mean-field model iṡ Comparison between the solution of (4) and numerical simulation. Three sets of results are shown, ω = 0 (top), ω = 1 (middle), and ω = 1.4 (bottom). Other parameters are β = 0.6, γ = 1, τ = 6, and k = 10. Simulation results are averaged across 100 iterations on random networks of 1000 nodes. All simulations begin by randomly selecting a node to infect at time t = 0. Simulation runs which die out are discarded and performed again.

When τ = 0, the dynamics of (4) are equivalent to the wellknown model of Gross et al. [8] . Figure 1 shows a comparison between the solution of the new model (4) and numerical simulation. The agreement is excellent despite the simplicity of the model and the fact that the moment closures do not reflect the changing network structure. In particular, both the solution and simulation results exhibit similar oscillatory behavior for the same parameter values. These results validate the model and allow us to analyze its behavior.

First, consider the basic reproductive ratio, R 0 , defined as the expected number of secondary infections caused by a single typical infectious node in an otherwise wholly susceptible population. One can find R 0 for the delayed rewiring model (4) N,0, k N,0,0,1) . Performing this analysis gives

Note that increasing the rewiring rate decreases the epidemic threshold R 0 , but the length of the delay, τ , has no effect on the threshold. However, as we will show later, it does affect the final outcome of the epidemic. System (4) also has an endemic steady state, but its value is determined by a transcendental equation which can only be solved numerically. Using this result in the numerical linear stability analysis of (4) allows us to analyze the stability of the endemic equilibrium. As shown in Fig. 2(a) , changes to both τ and ω are capable of destabilizing the endemic equilibrium. Regardless of the value of τ , eventually high values of the rewiring rate make the DFE stable again. For most values of τ this coincides with the point where the endemic steady state becomes biologically infeasible (less than or equal to zero), leaving the DFE as the only plausible steady state for the system. However, for sufficiently small values of τ , the endemic steady state remains feasible, and there is a small region of bistability. Qualitatively, this behavior is the same for any choice of the other parameters, as long as the endemic steady state remains biologically feasible, as illustrated for Other parameters are the same as in Fig. 2(a) . The endemic equilibrium is unstable in the red/yellow region, stable in the green/blue region, and biologically infeasible in the white region. different values of β in Fig. 3 . This figure shows that increasing the disease transmission rate allows the endemic steady state to be feasible for a wider range of link-cutting rate ω, and it also lowers the critical time delay τ , at which this steady state becomes unstable. Figure 2 (b) shows the endemic equilibrium, as well as the minima and maxima of oscillations for a range of β and ω values, with oscillations being observed in a significant part of the parameter space. One can clearly see the formation of an endemic bubble that has been discovered earlier in other epidemic models [28, 29] . Interestingly, both ω and β appear to play similar roles in the formation of the endemic bubble, namely, they open it through a supercritical Hopf bifurcation of the endemic equilibrium and then close it through a subcritical Hopf bifurcation.

Increasing the length of the delay can only induce a supercritical Hopf bifurcation, resulting in the emergence of stable oscillations, beyond which point larger values of τ only increase the amplitude of oscillations until it settles on some steady level, as shown in Fig. 2 (c). One should note that the minima of oscillations get closer to zero for larger τ , suggesting that for large rewiring times there are periods of time with negligible disease prevalence, followed by major outbreaks, as illustrated in Fig. 2(d) . In the limit τ → ∞, disconnected edges are never redrawn and the epidemic dies out, partially due to the network becoming sparser.

For the case without time delay, Gross et al. [8] found bistability in a large region of the parameter space, and periodic oscillations in a much smaller region. By contrast, results shown in Fig. 2 demonstrate a large region in the parameter space with oscillatory behavior. Delay differential equation (DDEs) are known to often produce oscillatory dynamics, and bubbles similar to those shown in Fig. 2(b) have been reported in other biological and epidemic models [28, 29] .

Let us now discuss the origins of oscillatory behavior in our model. The delay between disconnecting an edge and drawing a new one means that the total number of edges, and thus also the mean degree, is not constant. Whenever a susceptible node chooses to rewire, the total number of edges in the network decreases by two (since all edges are bidirectional) until time τ passes, and the edge is redrawn. The mean degree k(t) at any time t can be calculated directly from this argument as follows:

Figure 2(d) shows that oscillations are driven by the dynamics of k(t). During the early stages of an outbreak with a high rewiring rate, k(t) falls rapidly, as susceptible nodes cut links in response to the propagation of the disease. If the value of τ is large enough, then after a certain time the number of edges in the network is small enough to effectively starve the disease of transmission routes, and prevalence falls. These edges are then redrawn at the same rate as they were cut τ time ago, and k(t) grows, which allows the disease to spread again. Figure 2 (d) illustrates this behavior both in simulation and in the meanfield model (4) , showing how after the initial outbreak each new wave of infection is preceded by the recovery of network connectivity. The effect of oscillatory interactions between network connectivity and the propagating disease may be more pronounced in network simulations. Gross et al. [8] found that adaptive rewiring without delay can lead to the formation of highly connected clusters of susceptible nodes that are vulnerable to disease once any one node becomes infected. Since the model (4) does not account for changes in network structure, i.e., the closure is the same for all times and it does not depend on the average degree or degree distribution, this can potentially explain the small discrepancy between the solution of the deterministic and simulation models observed in Fig. 1 .

To get a better understanding of the interplay between network topology and dynamics, it is worth looking at how delayed rewiring alters degree distribution. Time snapshots of several large networks in Fig. 4 show the evolution of the degree distribution at various key points of an epidemic in an oscillatory regime. The initial network topology (black solid lines) is quickly reorganized to a peaked distribution. The oscillations in prevalence cause slight but repeated changes in the degree distribution. Unsurprisingly, when prevalence is at or near its peak, nodes with a lower degree are more common. When the prevalence falls, the distribution curves shift to the right, and the shape of the distribution flattens slightly. When the endemic steady state is stable, the degree distribution stabilizes to a peaked distribution between the two extremes of the oscillatory regime. A very important observation is that, irrespective of the initial network topology, due to rewiring different networks eventually settle on a very similar skewed degree distribution. This implies that earlier conclusions derived for the specific closure (1) appropriate for Erdős-Rényi graphs are actually applicable to modeling long-term dynamics of different types of networks, for which the influence of the initial topology is low since a significant amount of rewiring has already taken place.

The particular strength of this model lies in its ability to exhibit rich behavior from a simple system of DDEs. Time delay captures the fact that finding alternative contacts takes time, and also during an epidemic many people try to temporarily reduce the number of their contacts. Such behavior can be modeled using this delayed rewiring process. Previous work separated the processes of edge destruction and creation, and with edge creation occurring at a fixed rate the number of edges in the network was bounded only by the network size [24, 30] . In the model presented above, edge creation is reduced to replenishing global network connectivity towards its original level. Therefore, this model is fundamentally different from earlier models, even when parameters are matched.

During the initial growth phase it is the rate at which potential transmission is avoided by cutting a link, not the delay before drawing a new edge, that determines whether a major outbreak will occur. Although the delay does not affect the basic reproductive ratio R 0 , it does impact the outcome of the epidemic [see Fig. 2(c) ]. The result of introducing the delay is that oscillations occur in a large region of the parameter space. This happens due to the interplay between the spread of the disease and the behavioral changes in response to the epidemic. When the length of the delay is significant, the network becomes more sparse, healthy individuals are at lower risk of infection, and over time the prevalence falls. When the new edges are then formed, the disease is once again able to spread, and the cycle repeats.

Understanding the nature and cause of oscillations may provide opportunities to eradicate the disease. For example, if public awareness campaigns can lead to an increase in the length of the delay, the prevalence of the disease will naturally fall close to zero, at which time a relatively minor intervention, such as quarantining those who remain infected, may be enough to eradicate the disease from the population entirely. Currently, the model assumes that only susceptible nodes rewire. However, in reality, infected nodes are also likely to change their behavior. Risau-Gusman and Zanette [10] considered a model of rewiring where infected nodes rewire with a given probability. It would be of great value to examine a similar situation under delayed rewiring, with time delay representing the time for which infected nodes partially isolate themselves before rewiring, in accordance with advice given by public health authorities. This would alter the nature of the variable x(t) in the model. For example, if only infected nodes rewire, x(t) ≈ e −γ τ . Preliminary tests of this rewiring scheme show behavior similar to the present model.

Numerical simulations have shown that a similar oscillatory behavior is observed for other initial network topologies, including scale-free networks. Furthermore, since rewiring nodes choose their new neighbors uniformly at random from all available susceptible nodes, the initial network topology itself is transient, as shown in Fig. 4 , and, as a result, over time our model becomes more relevant. Future work will look at how the degree distribution and oscillations are affected in the case when the network links are rewired not randomly but according to a preferential attachment or some fitness-based rule. This could result in some interesting new dynamics due to the competition between the increased probability of highly connected nodes receiving new links, and the increased probability of infection.

Proc. Natl. Acad. Sci. USA

Correlation equations and pair approximations for spatial ecologies

The authors are grateful to anonymous reviewers for their helpful comments and suggestions. N.S. acknowledges funding for his Ph.D. studies from the Engineering and Physical Sciences Research Council, Grant No. EP/M506667/1, and the University of Sussex.