# Large epidemic thresholds emerge in heterogeneous networks of heterogeneous nodes

###### Abstract

One of the famous results of network science states that networks with heterogeneous connectivity are more susceptible to epidemic spreading than their more homogeneous counterparts. In particular, in networks of identical nodes it has been shown that network heterogeneity, i.e. a broad degree distribution, can lower the epidemic threshold at which epidemics can invade the system. Network heterogeneity can thus allow diseases with lower transmission probabilities to persist and spread. However, it has been pointed out that networks in which the properties of nodes are intrinsically heterogeneous can be very resilient to disease spreading. Heterogeneity in structure can enhance or diminish the resilience of networks with heterogeneous nodes, depending on the correlations between the topological and intrinsic properties. Here, we consider a plausible scenario where people have intrinsic differences in susceptibility and adapt their social network structure to the presence of the disease. We show that the resilience of networks with heterogeneous connectivity can surpass those of networks with homogeneous connectivity. For epidemiology, this implies that network heterogeneity should not be studied in isolation, it is instead the heterogeneity of infection risk that determines the likelihood of outbreaks.

In the exploration of complex systems, epidemiology plays an important role as a source for toy models and case studies, but also an area where a real world impact can be made vespignani2012modelling (); barrat2008dynamical (); hufnagel2004forecast (); brockmann2005dynamics (); colizza2007predictability (); janssen2006toward (). It has been pointed out that new diseases have emerged whenever environmental change brought humans in contact with new pathogen or disease vectors, i.e. animal hosts of a given disease karlen1996man (). The past decades have brought rapid environmental change, a growing world population, and increasing long-range connectivity in relevant networks, due to human travel and livestock transports brockmann2013hidden (). Together with decreasing vaccination levels and misuse of antibiotics, this has led to both emergence of new diseases and the return of old ones, sometimes in the form of highly resistant strains. In the future antibiotics, vaccinations, and quarantine are bound to remain first line of defence against these threats. However, insights from physics that can improve the efficiency of these measures, even if only by a small measure, have the potential to save many lives and relieve the economic burden created by the disease.

Many current studies seek to determine the so-called epidemic threshold, the critical level of infectivity that a pathogen needs to surpass to spread and cause large outbreaks castellano2010thresholds (); boguna2013nature (). This threshold depends on many factors including the structure of the underlying network of contacts, the heterogeneity in the host population, and behavioral responses to the disease. Among these, the effect of network structure is perhaps best understood pastor2001epidemic (); colizza2007invasion (); parshani2010epidemic (); belik2011natural (). It can be shown that the ability of a disease to spread is generally related to the leading eigenvalue of the networks adjacency or non-backtracking matrices wang2003epidemic (); castellano2010thresholds (); hamilton2014tight (); karrer2014percolation (); rogers2015assessing (). Factors that increase this eigenvalue, lead to lower epidemic thresholds. Two well-known factors that facilitate the spreading of diseases are high network connectivity and network heterogeneity. Here, network heterogenity can be defined as the second moment of the degree distribution, the probability distribution of the number of links on a randomly chosen node. In extreme cases, for example in scale-free networks, the epidemic threshold may vanish in the thermodynamic limit such that diseases with arbitrarily low infectivity can still spread pastor2001epidemic (); demirel2012absence ().

With respect to individual heterogeneity in the host population, some questions remain open. Generally intra-individual and link heterogeneity yang2012epidemic (); wang2014epidemic (), such as different levels of resistance, times of contact, differences in infectitivity or recovery rates, reduce the size and risk of epidemics miller2007epidemic (); miller2008bounding (); neri2011heterogeneity (); neri2011effect (); katriel2012size (); smilkov2014beyond (). For instance, it was shown analytically and numerically that epidemics are most likely if infectivity is homogeneous and least likely if the variance of infectivity is maximized miller2007epidemic (); miller2008bounding (). Comparable results were found in lattice models and in a biological experiment neri2011heterogeneity (). However, heterogeneous susceptibility can make networks more vulnerable to the spread of diseases if the correlation between a node’s degree and susceptibility are positive smilkov2014beyond (). This is intuitive as the intra-individual heterogeniety, in this case, amplifies the negative effect of network heterogeniety. Finally, it was reported that heterogeneous susceptibility can in some scenarios cause a secondary epidemic after a primary outbreak katriel2012size (); rodrigues2009heterogeneity ().

The studies referenced above focussed on epidemics on static networks. Another active line or research concerns the effect of dynamic, adaptive volz2009epidemic (); taylor2012epidemic (); gorochowski2012evolving (); li2010epidemic () or temporal networks holme2012temporal (); cui2014efficient (). While these types of networks are closely related, dynamical networks emphasize the overall statistical effect of changing connectivity volz2009epidemic (); taylor2012epidemic (), adaptive networks emphasize the dynamical response of network structure to the disease stategross2008adaptive (); gross2009adaptive (), and temporal networks focus on the impact of the specific timings of events perra2012activity (); lentz2013unfolding ().

An area that is so far poorly understood is how the behavioral response of individuals to epidemics reshapes the network. Beside institutionalised responses, such as mandatory vaccinations and quarantine, humans react to the outbreak of a major epidemic in a variety of ways, for instance by increasing hygiene, using protective measures such as face masks, reducing contact with infected individuals zhang2014suppression (), or avoiding contacts with other humans when infected. The common element in these responses is a reduction in the frequency of contacts between infectious and susceptible individuals. This limits the disease propagation by decreasing the effective network connectivity. However, the ultimate effect of behavioral responses cannot be understood as a static reduction of connectivity. If individuals respond to the epidemic state of other individuals by altering their interactions then a complex dynamical feedback loop is formed in which the state of nodes affects the evolution of the topology, while the topology governs the dynamics of the nodes’ state. Thus an adaptive network is created.

In the context of epidemics it has been shown in a simple model gross2006epidemic () that trying to avoid contact with infected is highly effective against small outbreaks but less effective against an established epidemic. Instead of a single epidemic threshold, such models possess an invasion threshold for new diseases and a (lower) persistence threshold for established epidemics. Thus, a parameter region is formed where a disease cannot spread when it is newly introduced, but can persist when it is already established. Physically speaking, a hysteresis loop is formed and the percolation transition at the onset of the epidemic becomes discontinuous. Under certain conditions the behavioural feedback loop can also lead to the emergence of epidemic cycles gross2006epidemic (); gross2008adaptive ().

The analysis of wide variety of related models showed that these observations are robust over a wide class of models zanette2008infection (); risau2009contact (); janssen2006toward (); shaw2008fluctuating (); schwartz2010rewiring (); wang2011epidemic (). Furthermore, studies showed that behavioral response can increase the impact of targeted vaccinations shaw2010enhanced (), and that the timing of interventions can be more important than in static networks yang2012efficient ().

In this paper, we study the combined effect of heterogenity in intra-individual parameters and the behavioral response. We show that, starting from a well-mixed network, a heterogeneous connectivity is formed. It is known that heterogeneous networks of heterogeneous nodes can be very resistant to disease outbreaks, if certain correlations are presentsmilkov2014beyond (). Here we show that these correlations naturally arise in the adaptive network, and that the resulting network configuration is generally significantly more resistant to outbreaks than a network with homogeneous topology. Our analysis suggests that the decisive property governing disease invasion is not network heterogeneity but the heterogeneity of the effective disease risk of agents.

## Results

Heterogeneous adaptive SIS model. We consider a network of agents connected by bidirectional links. We distinguish two types of agents, which we denote as type A and type B, which differ by their resistance to the disease. The type is an internal property of the agent that does not change. Hence, the proportion of agents of type A, , and the proportion of agents of type B, , are constant in time. For conciseness we denote the total number of agents in a given state by . Moreover, nodes carry an additional internal variable that indicates their epidemic state. We consider a variant of the susceptible-infected-susceptible (SIS) model of epidemic diseases anderson1991infectious (), such that a given node is either susceptible to the disease, state S, or infected (and infectious), state I. The proportion of nodes in the S and I state is denoted as and , respectively.

The network is initialized as an Erdős-Renyí random graph, . These networks have a narrow, Poissonian degree distribution, such that the network connectivity is homogeneous in the initial state. Each node is randomly assigned a type and epidemic state such that desired values of and the desired initial values of and are realized. Time evolution of the network is then driven by three processes, namely the a) recovery of infected nodes, b) contact avoidance, and c) contagion. These are implemented as follows: a) Infected nodes recover at rate , returning to the susceptible state. b) A given link, connecting a susceptible agent to an infected agent is rewired at rate . In an rewiring event the original link is cut and a new link between the susceptible node and a randomly chosen other susceptible node is created. c) For every link connecting a susceptible to an infected node, the disease is transmitted along the link at a rate that is dependent on the type of the susceptible node.

In the following we assume that such that nodes of type A are more susceptible to the disease than nodes of type B. Our mathematical results hold for parameters in arbitrary units, but the rates can be thought of having the dimension of nodes-per-time (recovery) and links-per-time (contagion, rewiring) respectively. Throughout the paper we balance the parameters and such that . The variance of susceptibility is .

Most of our results below are found by analytical calculation or continuation of solution branches and are thus non-simulative in nature. However, we compare these results to large agent-based simulations. In these simulations we use an event-driven (Gillespie) algorithm to simulate the stochastic process described above. In simulation we use an Erdős-Rényi network with nodes and links with recovery rate , rewiring rate and average susceptibility , unless noted otherwise. We vary as a proxy for heterogeneity. For every choice of and , the parameter is set such that . All parameters used in simulation runs are stated in the caption of the respective figure.

Hysteresis in the heterogeneous model. To gain some basic intuition, let us first investigate the system by explicit agent-based simulation of the network model zschaler2013largenet2 (). For this purpose we evolve the system in time, according to the stochastic rules, until a stationary level of epidemic prevalence is approached. Repeating this procedure for different values of infectivity reveals the diagram shown in Fig. 1. The Figure is qualitatively similar to results from the homogeneous adaptive SIS model: Epidemics starting from a small proportion of initially infected agents go extinct deterministically if the infectivity is below a certain threshold, , which we identify as the invasion threshold. By contrast, established epidemics, i.e. simulation runs starting with a higher initial proportion of infected, can persist if the infectivity surpasses a different persistence threshold .

The persistence threshold is lower than the invasion threshold, such that a bistable region is created. In this region an epidemic that enters the system at low density goes extinct, but an established large epidemic, that perhaps entered the system earlier when parameters were different, can persist. The bistable region constitutes a hysteresis loop: If we slowly increase the infectivity the extinct state is stable until , where a jump to the endemic state occurs. If we lower again, the system persists in the endemic state up to , where it collapses back down to the extinct state.

In the following we explore the effect of heterogeneity on the thresholds and . To gain analytical insights, we consider a moment expansion of the system demirel2014moment (). We use symbols of the form and with and to respectively denote the proportion of agents and per capita density of links between agents of a given type. For instance is the proportion of agents that are infected and of type A, and is the per capita density of links between susceptible agents of type A and infected agents of type B. All of these variables are normalized with respect to the total number of nodes . Given the number of infected nodes of a given type we can thus find the number of susceptible nodes by using the conservation law .

The time evolution of the proportion of nodes that are infected and of type A and B can be respectively written as

(1) |

(2) |

For the link densities, using a pair-approximation leads to equations of the form

(3) |

where the terms on the right hand side describe the impact of the different processes on the motif considered, in this example. For instance the first term corresponds to the creation of --links due to recovery of the infected node in --links. In total the -nodes recover at the rate . Every such recovery event creates an expected number of --links that is identical to the average number of --links anchored on an -node, which is . In summary, the change in the density of --links due to recovery of -nodes is , which explains the first term in Eq. (3).

In total the moment expansion yields a system of 11 ordinary differential equations. For conciseness we show the remaining equations in the Methods.

We solve the moment equations by numerical continuation of solution branches using AUTO doedel2010auto (). This reveals branches of stable and unstable steady states (Fig. 1). As in homogeneous adaptive SIS model the limits of the hysteresis loop are marked by a fold bifurcation and a transcritical bifurcation point. In the fold bifurcation the endemic steady state collides with an unstable saddle and the two states annihilate. In the transcritical bifurcation the saddle state intersects the healthy steady state, which causes the healthy state to become unstable. The value of in the fold bifurcation point thus marks the persistence threshold and the critical value in the transcritical bifurcation point marks the invasion threshold .

The comparison between the continuation results and agent-based simulation, in Fig. 1, shows that both methods are in good agreement regarding the location of the solution branches. Also the values for the persistence threshold agree. However, the continuation predicts a much higher value of the invasion threshold than the agent-based simulation.

To understand the discrepancy in the thresholds let us again consider the plot in more detail. The diagram shows a projection of the full 11-dimensional space spanned by the moment equations. In general bifurcation diagrams of this type identify the steady states uniquely. However, this is not the case when the epidemic is extinct. If there are no infected left, then the dynamics freezes independently of the connectivity of the nodes. Thus, as long as the state under consideration is stationary regardless of the values of and . Hence the zero line in the bifurcation diagram is really a 2-dimensional plane of absorbing steady states.

We can now resolve the discrepancy between the continuation and the simulation results. Continuation of the unstable solution branch leads to a specific point on the plane of extinct states. This landing point is uniquely defined and marks the invasion threshold of the network with the corresponding values of and . However, these values are not identical to those considered in the numerical simulation.

We argue that both of the invasion thresholds have significance. The simulation is valid for the well-mixed initial system, where the network structure has not yet adapted to the presence of the disease. Thus this threshold is relevant in case of the arrival of a new disease. By contrast the threshold found by continuation corresponds to a case where the network structure has adapted to the disease, for instance due to repeated exposure to same or similar pathogens. In the following, we therefore refer to the two thresholds as the initial and adapted invasion threshold, respectively.

Invasion thresholds and heterogeneity. We emphasize that the observed discrepeancy between the initial and the adapted invasion threshold could not appear in networks of identical nodes. For identical nodes the extinct state is unique on the level of the pair approximation (), and thus both thresholds must coincide. The results in Fig. 2 show that is indeed the case, while different thresholds are observed in all networks with heterogeneous nodes.

We note that the adapted invasion threshold is always higher than the initial threshold. Thus adaptation increases the robustness of the system to disease invasion. Let us therefore explore the adaptation in more detail. The adaptation is driven by the rewiring which means that nodes that are frequently infected in average lose links, while nodes that are rarely infected gain links. On the population level this means that the average degree of nodes of type A, , and the average degree of nodes of type B, change dynamically in response to the average proportion of nodes of each type that are infected. Before trying to compute the ratio let us first point out that a lower bound is , in this case the degrees are equal, and thus nodes of type A would be infected more frequently, due to their higher susceptibility. Hence would decrease and the ratio would increase. An upper bound is provided by . In this case, implies the that both types of nodes get infected at equal rate. Because the degree of nodes of type B is higher than the degree of nodes of type A, an infected node of type B will lose links more rapidly and hence the ratio will decrease.

The numerical value of the degree ration is shown in Fig. 3. To gain also an analytical understanding we resort to a description of the system that is coarser-grained than the full moment expansion. First, note that, on a population level , the mean degree of nodes of type , obeys a differential equation of the form

(4) |

where the two terms denote rewiring losses and gains, and auxialliary variables and have been defined to contain all other factors which do not depend on the index . In the steady state we find

(5) |

Thus, when we compute the ratio

(6) |

the factors and vanish. Using a mean field approximation, the epidemic state variables follow equations of the form

(7) |

where the terms capture the effects of contagion and recovery, respectively, and again auxiliary variables and have been defined that contain all other factors that don’t explicitly depend on . We use the same trick as before and consider the steady state, where

(8) |

and hence

(9) |

Substituting this result into Eq. (6) we find

(10) |

and hence

(11) |

The result of this mean field argument is in good agreement with numerical results (Fig. 3). It implies that the rewiring mechanism considered drives the system to a state where the less susceptible nodes (type B) have a higher mean degree, which partly, but not fully, compensates for their lower susceptibility. Thus in the adapted network the nodes of type B get infected more often than they would in a network with homogeneous degree distribution, but still less often than the nodes of type A.

To verify that the self-organization of the link distribution explains the observed discrepancy between the initial and the adapted invasion threshold, we turn to the agent-based simulation again. However, in this case we start the simulation from an artificially created adapted state. To initialize this state we simulate the system with the same set of initial parameters until the system reaches the stationary state. Then we retain the self-organized link pattern, but reassign all epidemic states, such that all agents are susceptible except for 20 initially infected. Then we simulate the system again until it either reaches the endemic state or the epidemic goes extinct. We locate the epidemic threshold by running a series of such simulations and find the point where the probability to reach the disease-free state becomes zero. The epidemic threshold that is thus found coincides with the result from continuation of the equation-based model (Fig. 4).

Higher thresholds in heterogeneous networks. Let us emphasize that the adapted network has a more heterogeneous degree distribution (Fig. 5) but also a higher epidemic threshold. Thus comparing the adapted and the maximally random state, the more heterogeneous network is actually more robust to the invasion of the pathogen. This appears counter-intuitive in the light of many recent work in network epidemiology, which showed that more heterogeneous contact networks are generally more susceptible to disease invasion. However, it was already shown in smilkov2014beyond () that heterogeneous networks become less susceptible if there is a correlation such that the highly connected nodes have lower susceptibility to the disease. The present results show that a simple but plausible local adaptation rule can drive this process so far that the network with heterogeneous degree distribution is more robust against disease invasion than an Erdős-Renyi random graph of the same mean degree.

Link heterogeneity is clearly a double-edged sword. On the one hand it is intuitive that some amount of heterogeneity in the degree is advantageous if it means that more links lead to nodes with low susceptibility. On the other hand, epidemic theory suggests that in the limit of very heterogeneous networks, the robustness of the network should decline because of the high excess degreepastor2001epidemic (); moreno2002epidemic (). This suggests that there should be an optimum level of heterogeneity, where the epidemic invasion threshold is highest.

We can understand the interplay between the two effects by a link-centric percolation argument. We consider the limit of low disease prevalence and focus on the active links, i.e. links connecting an infected node to a susceptible node. Given a single focal active link we can estimate the expected number of secondary active that is created by transmission along the focal link by

(12) |

where is the constant mean degree of the system. Defining and substituting and we can express the link reproductive number as a function of the degree ratio . Although the resulting expression for is relatively complex, we can find it’s minimum, i.e. the most robust point, by differentiating, which yields

(13) |

This coincides with the upper bound for the degree ratio, or, in other words, the point where the nodes of types A and B are infected at an identical rate. Therefore increasing the degree heterogeneity by increasing the degree ratio is advantageous to the point where the more resistant agents become infected more often than their less resistant counterparts. Thus it appears that the decisive characteristic that determines its robustness to disease invasion is not the heterogeneity of the the network structure, but the heterogeneity of the disease risk to which the agents are exposed.

The independence of the result above from other parameters suggests that it is true in a wider class of systems, but this intuition will need to be validated in further investigations. For the adaptive system this result means that the network always operates in the regime where a higher degree ratio and therefore more heterogeneity has a stabilizing effect. We could in principle construct a system that self-organizes to the optimal point, by replacing the per-link rewiring rate by a rewiring rate per infected node. However this variant of the model is beyond the scope of the present paper.

## Discussion

In this paper we studied an adaptive heterogeneous SIS model numerically and analytically. The analysis revealed that heterogeneity in the intrinsic parameters of the nodes induces heterogeneity in the connectivity: Over time more resistant nodes gain more links until a steady state is reached, in which nodes with higher resistance are still less likely to contract the disease, but more highly connected than average nodes.

A well known result in network science is that more heterogeneous networks are less resistant to the invasion of diseases. However, in the self-organized networks studied here the opposite effect is observed. In comparison to random networks the self-organized networks are both more resistant to the disease and more heterogeneous in connectivity.

The evolved networks gain their resistance to the disease from correlations between intrinsic parameters and the node connectivity. The increased resistance thus arises directly from the heterogeneity. A comparable effect is not possible in networks of identical nodes. While the specific bifurcation structure of the present model may depend on modelling assumptions, the basic interplay between the connectivity and the invasion threshold arises from the fundamental physics of the spreading process and is thus likely to be a generic feature that is observed across many models. It thus appears plausible that also in real world epidemics some anti-correlation between the true susceptibility and the effective degree of agents will be induced. By concentrating more of the remaining links of the more resilient individually the network will generally become both more heterogeneous and more resistant to the disease.

It is possible that in real networks the self-organized heterogeneity is relatively minor compared to network heterogeneity that is unrelated to the epidemic in question. In that case the phenomenon described here would still occur but be of lesser importance. The extent to which self-organized heterogeneity plays a role in real world diseases will certainly depend on the disease in question, and can probably only be assessed through further empirical work.

For epidemiology the results reported here imply, that network heterogeneity should not be considered in isolation. If there is an underlying heterogeneity in the susceptibility to the disease then a heterogeneous network may be more resistant to disease invasion than its homogeneous counterpart. Moreover, a vaccination campaign that targets the most highly connected nodes may end up vaccinating the wrong people as these nodes may also have the strongest natural resistance against the disease.

Perhaps most importantly, our results suggest that the invasibility of the network is not governed by the heterogenity of the network alone, but by the heterogeneity of effective disease risk, which takes both node degree and susceptibility into account.

## Methods

Moment-closure approximation. The time evolution of the proportion of nodes that are infected and of type can then be written as

(14) |

The two terms of this equation capture the recovery of infected nodes and the infection of susceptible nodes, respectively. In the second term we used the symbol to denote the number of links between a node of state S and type and a node of state I and type , normalized with respect to the total number of nodes. This quantity therefore has the dimension of links per node. For determining these link densities, we can write additional evolution equations, which are in turn depend on the density of triplet chains of nodes of given type and state , this yields

(15) |

(16) |

(17) |

where

(18) |

In the equation for , the S-S-links, the first term accounts for the creation of --links by recovery of an infected node in an S-I-link. The factor needs to be included to avoid double counting in the case of identical indices. The second term accounts for the destruction of --links due to infection of one of the two S-nodes by a node that is external to the link. The number of such infected nodes outside the S-S-link is given by the number of triplets and . For a single --link the infection rate due to infected connected to the -node is , i.e. the expected number of triplet chains that run through one specific the S-S-link and include an infected node on the side, multiplied by the effective infection rate . To obtain the total rate we multiply by the number of those links, , which cancels the denominator, take the infected on the other side of the S-S-link into account, and multiply to avoid double-counting. The final term in the equation accounts for creation of S-S-links by rewiring. Links of the type - are rewired at rate . When rewiring the susceptible node cuts the link to the infected and connects to a randomly chosen susceptible. In such a rewiring event, a new --link is created if the newly-chosen partner is of type A. This is the case with probability . Such that the total rate of --link creation from --link rewiring is . Taking all possible combinations of indices and double-counting into account leads to the term in the equation.

In the equation for --links the first term accounts for the loss of these links due to recovery of one of the linked nodes. The second term, , accounts for the creation of these links due to transmission of infection inside an S-I-link. Similarly, the final term accounts for the creation of I-I-link by infection of the S-node in an S-I-link from an internal source. In this term triplet chains appear in analogy to the term for the loss of S-S-links due to infection from sources external to the link, discussed above. In the equation for the S-I-links, the first term accounts for the creation of these links due to recovery of one infected node in an I-I-link. The second term accounts for the loss of S-I-links, both due to recovery of the I-node and due to internal transmission of the infection, resulting in an I-I-link. The fourth term captures the creation of S-I-links due to infection of an S-node in an S-S-link, whereas the fourth term captures the loss of S-I-links due to infection of the S-node from a source external to the link. Finally, the last term accounts for the loss of S-I-links due to rewiring.

To cut the progression to ever larger network motifs one approximates the density of triplet chains by a moment closure approximation, here the pair approximation

(19) |

where is a factor arising from symmetries, such that if , if either or , and if .

Inherent in the moment-closure approximation is the assumption that long-ranged correlations vanish, such that the densities of motifs beyond the cut-off conform to statistical expectations. This assumption is the main source of inaccuracies in this type of approximation demirel2014moment (). The approximation can still be used to identify phase transitions in the adaptive SIS model as the correlations associated with these are still captured. However, the approximation performs poorly in fragmentation-type transitions, found for instance in the adaptive voter modelbohme2011analytical ().

The additional equations from the moment closure approximation that are used in the continuation are

(20) |

(21) |

(22) |

(23) |

(24) |

(25) |

(26) |

(27) |

## References

- (1) Vespignani, A. Modelling dynamical processes in complex socio-technical systems. Nat. Phys. 8, 32–39 (2012).
- (2) Barrat, A., Barthelemy, M., & Vespignani, A. Dynamical processes on complex networks, Vol. 1. (Cambridge University Press, Cambridge, 2008).
- (3) Hufnagel, L., Brockmann, D., & Geisel, T. Forecast and control of epidemics in a globalized world. Proc. Natl. Acad. Sci. USA 101, 15124–15129 (2004).
- (4) Brockmann, D. et al. Dynamics of modern epidemics. SARS: A case study in emerging infections, Ch. 11, 81-91, 2005.
- (5) Colizza, V., Barrat, A., Barthélemy, M., & Vespignani, A. Predictability and epidemic pathways in global outbreaks of infectious diseases: the sars case study. BMC medicine 5, 34 (2007).
- (6) Janssen, M. A. et al. Toward a network perspective of the study of resilience in social-ecological systems. Ecology and Society 11, 15 (2006).
- (7) Karlen, A. Man and microbes: Disease and plagues in history and modern times. (Simon and Schuster, 1996).
- (8) Brockmann, D. & Helbing, D. The hidden geometry of complex, network-driven contagion phenomena. Science 342, 1337-1342 (2013).
- (9) Castellano, C. & Pastor-Satorras, R. Thresholds for epidemic spreading in networks. Phys. Rev. Lett. 105, 218701 (2010).
- (10) Boguñá, M., Castellano, C. & Pastor-Satorras, R. Nature of the epidemic threshold for the susceptible-infected-susceptible dynamics in networks. Phys. Rev. Lett. 111, 068701 (2013).
- (11) Pastor-Satorras, R. & Vespignani, A. Epidemic spreading in scale-free networks. Phys. Rev. Lett. 86, 3200 (2001).
- (12) Colizza, V. & Vespignani, A. Invasion threshold in heterogeneous metapopulation networks. Phys. Rev. Lett. 99, 148701 (2007).
- (13) Parshani, R., Carmi, S. & Havlin, S. Epidemic threshold for the susceptible-infectious-susceptible model on random networks. Phys. Rev. Lett. 104, 258701 (2010).
- (14) Belik, V., Geisel, T. & Brockmann, D. Natural human mobility patterns and spatial spread of infectious diseases. Phys. Rev. X 1, 011001 (2011).
- (15) Wang, Y., Chakrabarti, D., Wang, C. X. & Faloutsos, C. Epidemic spreading in real networks: An eigenvalue viewpoint. In Reliable Distributed Systems, 2003. Proceedings. 22nd International Symposium on, IEEE 25-34 (2003).
- (16) Hamilton, K. E. & Pryadko, L. P. Tight lower bound for percolation threshold on an infinite graph. Phys. Rev. L 113, 208701 (2014).
- (17) Karrer, B., Newman, M. E. J. & ZdeborovÃ¡, L. Percolation on sparse networks. Phys. Rev. L 113, 208702 (2014).
- (18) Rogers, T. Assessing node risk and vulnerability in epidemics on networks. Europhys. Lett. 109, 28005 (2015)
- (19) Demirel, G. & Gross, T. Absence of epidemic thresholds in a growing adaptive network. arXiv preprint 1209.2541 (2012).
- (20) Yang, Z. M. & Zhou, T. Epidemic spreading in weighted networks: An edge-based mean-field solution. Phys. Rev. E 85, 056106 (2012).
- (21) Wang, W., Tang, M., Zhang, H. F., Gao, H., Do, Y. & Liu, Z. H. Epidemic spreading on complex networks with general degree and weight distributions. Phys. Rev. E 90, 042803 (2014).
- (22) Miller, J. C. Epidemic size and probability in populations with heterogeneous infectivity and susceptibility. Phys. Rev. E 76, 010101 (2007).
- (23) Miller, J. C. Bounding the size and probability of epidemics on networks. J. Appl. Prob. 45, 498-512 (2008).
- (24) Neri, F. M., Pérez-Reche, F. J., Taraskin, S. N. & Gilligan, C. A. Heterogeneity in susceptible–infected–removed (sir) epidemics on lattices. J. R. Soc. Interface 8, 201–209 (2011).
- (25) Neri, F. M. et al. The effect of heterogeneity on invasion in spatial epidemics: from theory to experimental evidence in a model system. PLoS Comput. Bio. 7, e1002174 (2011).
- (26) Katriel, G. The size of epidemics in populations with heterogeneous susceptibility. J. Math. Biol 65, 237-262 (2012).
- (27) Smilkov, D., Hidalgo, C. A. & Kocarev, L. Beyond network structure: How heterogeneous susceptibility modulates the spread of epidemics. Sci. Rep. 4, 4795 (2014).
- (28) Rodrigues, P., Margheri, A., Rebelo, C. & Gomes, M. G. M. Heterogeneity in susceptibility to infection can explain high reinfection rates. J. Theor. Biol 259, 280-290 (2009).
- (29) Volz, E & Meyers, L. A. Epidemic thresholds in dynamic contact networks. J. R. Soc. Interface 6, 233-241 (2009).
- (30) Taylor, M., Taylor, T. J. & Kiss, I. Z. Epidemic threshold and control in a dynamic network. Phys. Rev. E 85, 016103 (2012).
- (31) Gorochowski, T. E., Bernardo, M. D. & Grierson, C. S. Evolving dynamical networks: A formalism for describing complex systems. Complexity 17, 18-25 (2012).
- (32) Li, X., Cao, L. & Cao, G. F. Epidemic prevalence on random mobile dynamical networks: individual heterogeneity and correlation. Eur. Phys. J. B 75, 319-326 (2010).
- (33) Holme, P. & Saramäki, J. Temporal networks. Phys. rep. 519, 97-125 (2012).
- (34) Cui, A. X., Wang, W., Tang, M., Fu, Y., Liang, X. M. & Do, Y. Efficient allocation of heterogeneous response times in information spreading process. Chaos 24, 033113 (2014).
- (35) Gross, T. & Blasius, B. Adaptive coevolutionary networks: a review. J. R. Soc. Interface 5, 259-271 (2008).
- (36) Gross, T. & Sayama, H. Adaptive networks. (Springer, 2009).
- (37) Perra, N., Gonçalves, B., Pastor-Satorras, R. & Vespignani, A. Activity driven modeling of time varying networks. Sci. Rep. 2, 469 (2012).
- (38) Lentz, H. HK., Selhorst, T. & Sokolov, I. M. Unfolding accessibility provides a macroscopic approach to temporal networks. Phys. Rev. Lett. 110, 118701 (2013).
- (39) Zhang, H. F., Xie, J. R., Tang, M. & Lai, Y. C. Suppression of epidemic spreading in complex networks by local information based behavioral responses. Chaos 24, 043106 (2014).
- (40) Gross, T., DâLima, C. J. D. & Blasius, B. Epidemic dynamics on an adaptive network. Phys. Rev. Lett. 96, 208701 (2006).
- (41) Zanette, D. & Risau-Gusmán, S. Infection spreading in a population with evolving contacts. J. Biol. Phys. 34, 135-148 (2008).
- (42) Risau-Gusmán, S. & Zanette, D. Contact switching as a control strategy for epidemic outbreaks. J. Theor. Biol 257, 52-60 (2009).
- (43) Shaw, L. B. & Schwartz, I. B. Fluctuating epidemics on adaptive networks. Phys. Rev. E 77, 066101 (2008).
- (44) Schwartz, I. B. & Shaw, L. B. Rewiring for adaptation. Physics 3, 17 (2010).
- (45) Wang, B., Cao, L., Suzuki, H. & Aihara, K. Epidemic spread in adaptive networks with multitype agents. Physics A 44, 035101 (2011).
- (46) Shaw, L. B. & Schwartz, I. B. Enhanced vaccine control of epidemics in adaptive networks. Phys. Rev. E 81, 046120 (2010).
- (47) Yang, H., Tang, M. & Zhang, H. F. Efficient community-based control strategies in adaptive networks. New J. of Phys. 14, 123017 (2012).
- (48) Anderson, R. M. & May, R. M. Infectious diseases of humans, Vol. 1. (Oxford university press, Oxford, 1991).
- (49) Zschaler, G. & Gross, T. Largenet2: an object-oriented programming library for simulating large adaptive networks. Bioinformatics 29, 277-278 (2013).
- (50) Demirel, G., Vazquez, F., Böhme, G. A. & Gross, T. Moment-closure approximations for discrete adaptive networks. Physica D 267, 68-80 (2014).
- (51) Doedel, E. J. et al. Auto-07p: Continuation and bifurcation software for ordinary differential equations. (Concordia University, Montréal, 2010).
- (52) Moreno, Y., Pastor-Satorras, R. & Vespignani, A. Epidemic outbreaks in complex heterogeneous networks. Eur. Phys. J. B 26, 521-529 (2002).
- (53) Böhme, G. A. & Gross, T. Analytical calculation of fragmentation transitions in adaptive networks. Phys. Rev. E 83, 035101 (2011).

## Acknowledgement

This work was supported by the EPSRC under grant EP/K031686/1, National Natural Science Foundation of China (Grant Nos. 11105025, 91324002) and the Program of Outstanding Ph.D. Candidate in Academic Research by UESTC (Grant No. YBXSZC20131036).

## Author contributions

H. Y., M. T. and T. G. devised the research project. H. Y. performed numerical simulations. H. Y. and T. G. analysed the results. H. Y., M. T. and T. G. wrote the paper.

## Additional information

Competing financial interests: The authors declare no competing financial interests.

## Figure legends

Figure 1: Bifurcation diagram of the adaptive heterogeneous SIS model. Show is the stationary level of disease prevalence as a function of infectivity . When the infectivitiy is decreased the endemic state vanishes in a saddle node bifurcation (‘Persistence threshold’). The disease free state can be invaded by the epidemic, if the infectivity surpasses a point where a transcritical bifurcation occurs (‘Invasion threshold’). Agent-based simulation (circles) and equation-based continuation (lines) provide consistent results on the persistence threshold, but predict different invasion thresholds. This discrepancy appears due to a projection effect, because the state where the disease is extinct is not uniquely defined (see text). In addition to stable solution branches (solid) the continuation also reveals an unstable solution branch (dotted). Parameters: , , , , , .

Figure 2: Comparison of thresholds. The plot shows a very good agreement agreement between equation-based continuation (lines) and agent-based simulations (symbols) for the persistence thresholds (box, dashed). However, a notable difference exists for the invasion thresholds (circle, dotted). Parameters: , , , .

Figure 3: Network heterogeneity in the adapted state, indicated by the degree ratio . The dependence of the degree ratio in agent-based simulations (black circles) closely follows the prediction from integration of the ODE model (solid red line), a very good approximation is also provided by the relationship (blue dashed line), whereas the naive expectation (pink dotted line) overestimates the network heterogeneity significantly. Contrary to expectations the networks following the naive solution (the most heterogeneous case), would be maximally stable against disease invasion. Parameters: , , , . Inset: the magnification of the superimposed part of the main figure.

Figure 4: Comparison of thresholds in self-organized networks. The plot shows a very good agreement between equation-based continuation (lines) and agent-based simulations started in an artificially created adapted state (symbols) for both the invasion thresholds (circle, dotted) and the persistence thresholds (box, dashed). See Fig. 2 for comparison. Parameters: , , , .

Figure 5: Self-organized heterogeneity. Comparison of the degree distributions of the initial unadapted network used in the first set of simulations and the ‘adapted’ self-organized network. The self-organized network is significantly more heterogeneous. However, at the same time it is more resistant to disease invasion. Parameters: , , , , , .