Large epidemic thresholds emerge in heterogeneous networks of heterogeneous nodes

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.


Introduction
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 [1,2,3,4,5,6].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 [7].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 [8].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.
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 [9,10].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 [11,12,13,14].It can be shown that the ability of a disease to spread is generally related to the leading eigenvalue of the networks adjacency matrix [15,9].Factors that increase this eigenvalue, such as high connectivity, or strong heterogeneity in the number of links per node lead to lower epidemic thresholds.In extreme cases with very heterogeneous connectivity, the epidemic threshold may vanish in the thermodynamic limit such that diseases with arbitrarily low infectivity can still spread [11,16].
With respect to individual heterogeneity in the host population, some questions remain open.Generally intra-individual and link heterogeneity [17,18], such as different levels of resistance, times of contact, differences in infectitivity or recovery rates, reduce the size and risk of epidemics [19,20,21,22,23,24].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 [19,20].Comparable results where found in lattice models and in a biological experiment [21].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 [24].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 [23,25].
The studies referenced above focussed on epidemics on static networks another active line or research concerns the effect of dynamic [26,27,28,29] or temporal networks [30,31].While dynamical and temporal networks are closely related, dynamical networks emphasize the overall statistical effect of changing connectivity [26,27], while temporal networks often focus on the impact of the specific timings of events.For instance, temporal networks can not be approximated by their time-aggregated counterparts [32], but temporal accessibility graphs, which is built by consecutively adding paths of growing length according to the timing of events, can be used as a formalism to analyze temporal networks [33].
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 such measures as face masks, reducing contact with infected individuals [34], 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 induviduals 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.Systems that contain such a feedback loop are called adaptive networks [35,36] and exhibit a number of self-organization phenomena that do not exist in their nonadaptive counterparts.
In the context of epidemics it has been shown in a simple model [37] 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 behavioral feedback loop can also lead to the emergence of epidemic cycles [37,35].
The analysis of wide variety of related models showed that these observations are robust over a wide class of models [38,39,6,40,41,42].Furthermore, studies showed that behavioral response can increase the impact of targeted vaccinations [43], and that the timing of interventions can be more important than in static networks [44].
In this paper, we study the combined effect of heterogenity in intra-individual parameters and the behavioral response.This creates an adaptive network [36] where the heterogeneity in the nodes induces heterogeneous connectivity.We show that, this leads to networks that are more heterogeneous but at the same time have higher epidemic thresholds.Thus in the adaptive network the well known effect of network heterogeneity is reversed.We show that the decisive property governing disease invasion is not network heterogeneity but the heterogeneity of the effective disease risk of agents, which takes both connectivity and intrinsic parameters into account.

Heterogeneous adaptive SIS model
We consider a network of N agents connected by K 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, p a , and the proportion of agents of type B, p b = 1 − p a , are constant in time.For conciseness we denote the total number of agents in a given state i ∈ {a, b} by N i = Np i .Moreover, nodes carry an additional internal variable that indicates their epidemic state.We consider a variant of the susceptibleinfected-susceptible (SIS) model of epidemic diseases [45], 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 [S] and [I] = 1 − [S], respectively.
The network is initialized as an Erdős-Renyí random graph, G(N, K).Each node is randomly assigned a type and epidemic state such that desired values of p i and the desired initial values of [S] and [I] 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 βψ i that is dependent on the type i of the susceptible node.
In the following we assume that ψ a > ψ b such that nodes of type A are more susceptible to the disease than nodes of type B. Throughout the paper we balance the parameters ψ a and ψ b such that p a ψ a + p b ψ b = ψ .The variance of susceptibility is In simulation we use an Erdős-Rényi network with N = 10 5 nodes and K = 10 6 links with recovery rate µ = 0.002 and average susceptibility, unless noted otherwise.We vary ψ a as a proxy for heterogeneity.For every choice of ψ a , the parameter ψ b is set such that ψ = 1/2.

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 [46].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, β inv , 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 β per .
The persistence threshold is lower than the invasion threshold, such that a region is created, in which established epidemics can persists, but new diseases cannot invade.This region constitutes a hysteresis loop: If we slowly increase the infectivity the extinct state is stable until β inv , where a jump to the endemic state occurs.If we lower β again, the system persists in the endemic state up to β per , where it collapses back down to the extinct state.
In the following we explore the effect of heterogeneity on the thresholds β inv and β per .To gain analytical insights, we consider a moment expansion of the system [47].We use symbols of the form [X u ] and [X u Y v ] with X, Y ∈ {I, S} and u, v ∈ {a, b} to respectively denote the proportion of agents and per capita density of links between agents of a given type.For instance [I a ] is the proportion of agents that are infected and of type A, and [S a I b ] 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 N. Given the number of infected nodes of a given type we can thus find the 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 [S u I v ] to denote the number of links between a node of state S and type u and a node of state I and type v, 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 where In the equation for [S u S v ], the S-S-links, the first term i.e. the expected number of triplet chains that run through one specific the S-S-link and include an infected node on the S a side, multiplied by the effective infection rate βψ b .To obtain the total rate we multiply by the number of those links, [S a S b ], which cancels the denominator, take the infected on the other side of the S-S-link into account, and multiply κ 1 to avoid double-counting.The final term in the equation accounts for creation of S-S-links by rewiring.Links of the type [S a I b ] 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 [S a S a ] link is created if the newly-chosen partner is of type A. This is the case with probability ). Taking all possible combinations of indices and double-counting into account leads to the term in the equation.
In the equation for [I u I v ]-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-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 where δ is a factor arising from symmetries, such that δ = 4 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 [47].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 model [48].
Using the pair-approximation and writing out the resulting equations in long-form we obtain We solve these equations by numerical continuation of solution branches using AUTO [49].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 annihilates.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 β per and the critical value in the transcritical bifurcation point marks the invasion threshold β inv .
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 are in good agreement.However, the continuation predicts a much higher value of the invasion threshold than the agent-based simulation.
When one encounters a discrepancy between actual micro-level simulations and a theoretical prediction based on an approximation, one is generally inclined to trust the simulation and blame approximation errors for the discrepancy.However, here this is not the case.

Invasion thresholds and heterogeneity
As we show in the following both estimates of the invasion threshold are correct, albeit in different scenarios.Both the agreement of the different methods on the persistence threshold and their disagreement on the invasion threshold are robust features of the system.For instance plotting the thresholds as a function of the parameter ψ a reveals that the agent-based and equation-based persistence thresholds coincide, whereas the agent-based and equation-based invasion thresholds disagree whenever the system is heterogeneous (Fig. 2).
To understand the difference between the two thresholds in the heterogeneous system, consider that they apply in different situations.In the agent-based simulation we start the system from a mixed state, where the node type and the links are placed independently.This is reasonable for the case of a new pathogen arriving in a population, such that the contact structure has not yet adapted to its presence.While the subsequent simulation takes rewiring in the vicinity of infected agents into account, the success of the initial invasion is determined by a network topology whose macroscopic properties have not yet adjusted to the presence of the disease.
Conversely, the continuation results consider a long-term limit in which all properties of the network have adapted.In the following we refer to this case as the adapted state.Considering the invasion threshold in the adapted state is reasonable if the system has been exposed to the present or a closely related pathogen before the current outbreak.(Consider for instance the case of flu: While the next year's flu is likely to be a new strain, people are generally aware of their personal susceptibility to flu-like diseases and may implement measures to reduce avenues of infection, even before the arrival of the pathogen.) Let us discuss the impact of adaptation on the system.Consider first, a network in the initial state, where the nodes' degree and susceptibility are uncorrelated, such that k a = k b .In this network, the nodes of type B spend a lesser proportion of their time in the infected state than nodes of type A, due to their lower susceptibility.This means that nodes of type A have both a higher loss rate of links due to rewiring and a lower gain rate.One could naively expect that this process continues until nodes of both types spend the same proportion of their time in the infected state, i.e. ψ a k a = ψ b k b .However, this is not the case.To see this, let us now suppose we start the network in a state with ψ a k a = ψ b k b then the nodes of type A and B would get infected equally often, as the higher degree of type B nodes balances their lower susceptibility.In this case both types of nodes have the same gain from rewiring, but due to the higher mean degree of nodes of type B they lose links more rapidly than nodes of type A when they get infected.From these arguments we conclude that a loose lower bound for the degree ratio is k b /k a = 1, and an upper bound is Of the two effects that contribute to the higher degree of B-nodes, the lower loss rate is more important than the higher gain rate.For the sake of finding a simple approximation we thus assume that both types of nodes gain links at the same rate.This is exactly true in a variant of the model where rewired links connect to a random new partner independently of the epidemic state, but we do not consider this variant of the model here.Given the assumption, the system will approach a state in which also the loss rates of links of the two node types match up.At this point it is conducive to consider the scaling of the loss rate on a per-node basis.For epidemics close to the threshold, the rate at which a given node of type i is infected is approximately proportional to the product of its degree and susceptibility k i ψ i .Furthermore, once infected the expected loss of links due to rewiring is proportional to k i .All other factors that appear in the balance equation for the degree are independent of node type.Such that in the adapted state we find and therefore Comparing this result with numerical data shows a very good agreement (see Fig. 3).The result 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. See Fig. 2 for comparison.
To verify that the self-organization of the link distribution explains the observed discrepancy between the agent-based and equation-based thresholds, we turn to the agentbased 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 [24] that heterogeneous networks get 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.
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 suggest that in the limit of very heterogeneous networks, the robustness of the network should decline because of the high excess degree [11,50].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 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 where k = p a k a + p b k b is the constant mean degree of the system.Defining q = k b /k a and substituting k a = k/(p a + p b q) and k b = k/(p a q −1 + p b ) we can express the link reproductive number Z 0 as a function of the degree ratio q.Although the resulting expression for Z 0 is relatively complex, we can find it's minimum, i.e. the most robust point, by differentiating which yields 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.

Conclusions
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 degree.The increased resistance thus arises directly from the heterogeneity.A comparable effect would not be possible in homogeneous networks.
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 over many models.
For many common real world epidemics, such as for instance flue, it is plausible that individuals have some intuition about their personal susceptibility to the disease.For such diseases reducing the exposure to potentially infectious agents is feasible.A reduction of contacts can be achieved by measures such as increased hygiene, protective devices (e.g.facemasks), contact avoidance, or even travel to a different geographical region.However, all of these measures involve some cost to the agent.It is thus plausible that agents who know that they are very susceptible to the disease will employ these measures more readily than those who believe themselves to be less susceptible.Thus it is plausible that 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 only be of minor 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.

Figure 1 :
Figure 1: (Color online) Bifurcation diagram of the adaptive heterogeneous SIS model.
accounts for the creation of [S u S v ] links by recovery of an infected node in an S-I-link.The factor κ 1 needs to be included to avoid double counting in the case of identical indices.The second term accounts for the destruction of [S u S v ] 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-Slink is given by the number of triplets [S u S v I w ] and [S v S u I w ].For a single S a -S b -link the infection rate due to infected connected to the S b -node is βψ b ([S a S b I a ]+[S a S b I b ])/[S a S b ],

Figure 2 :
Figure 2: (Color online) Comparison of thresholds.The plot shows a very good agree-

Figure 3 :
Figure 3: (Color online) Network heterogeneity in the adapted state, indicated by the

Figure 4 :
Figure 4: (Color online) Comparison of thresholds in self-organized networks.The plot

Figure 2 :Figure 4 :
Figure 2: (Color online) Comparison of thresholds.The plot shows a very good agree-