Mutually cooperative epidemics on power-law networks.

The spread of an infectious disease can, in some cases, promote the propagation of other pathogens favoring violent outbreaks, which cause a discontinuous transition to an endemic state. The topology of the contact network plays a crucial role in these cooperative dynamics. We consider a susceptible-infected-removed-type model with two mutually cooperative pathogens: An individual already infected with one disease has an increased probability of getting infected by the other. We present a heterogeneous mean-field theoretical approach to the coinfection dynamics on generic uncorrelated power-law degree-distributed networks and validate its results by means of numerical simulations. We show that, when the second moment of the degree distribution is finite, the epidemic transition is continuous for low cooperativity, while it is discontinuous when cooperativity is sufficiently high. For scale-free networks, i.e., topologies with diverging second moment, the transition is instead always continuous. In this way we clarify the effect of heterogeneity and system size on the nature of the transition, and we validate the physical interpretation about the origin of the discontinuity.

INTRODUCTION

Modeling epidemic dynamics plays a key role in predicting disease outbreaks and designing effective strategies to prevent or control them. While traditional theories of disease propagation ignore network effects, there has been a large amount of research in the past decades aimed at understanding how the topological structure of contact networks affects the dynamics upon them [1]. Such studies have focused mainly on the dynamics of a single epidemic disease.

An issue of growing interest in current epidemiological research is how concurrent spreading diseases interact with each other. Such coupled spreading scenarios can occur when either multiple pathogens or multiple strains of the same disease simultaneously propagate in the same population. The complexity of the problem is largely increased already in the two-pathogen case. Two diseases circulating in the same host population can interact in many different ways, with either synergistic or antagonistic effects.

A well-known type of interaction is cross immunity: An individual infected with one disease becomes partially or fully immune to infection by the second one. In this case the two pathogens compete for the same population of hosts. The competition between epidemics that are mutually exclusive or antagonistic was studied in [2-5].

An opposite case that is recently gaining attention is the spreading of two or more cooperating pathogens: In this case, an individual that is already infected with one disease has increased chance of getting infected with another. A notable example is the 1918 Spanish flu pandemic that involved about one-third of the world's population killing tens of millions of people within months [6]. A considerable proportion of the infected were coinfected by pneumonia [7,8], and most deaths where caused not directly by the virus, but by the secondary bacterial infection. Another well-known case is that of HIV, which increases the host susceptibility to other pathogens, in particular to the hepatitis C virus (HCV) [9].

In coinfections, positive feedback between multiple diseases can lead to more rapid outbreaks. One of the most interesting questions is whether cooperation can change the epidemic transition from being continuous to being abrupt when external conditions vary, even slightly, as for a microscopic change in infectivity. This is an extremely relevant problem since the possibility to enact countermeasures critically depends on the type of transition: In the discontinuous case the epidemic can take over a population explosively without any of the warning signs which occur in single epidemics.

Recently, Newman and Ferrario [10] introduced a model of coinfection based on the susceptible-infected-removed (SIR) model [11] with two diseases diffusing over the same contact network. One disease spreads freely while the second can infect only individuals already infected with the first. The authors assume total asymmetry and time-scale separation in the spreading of the two pathogens. The system displays two epidemic thresholds, the first one being the usual SIR threshold, the second occurring when the fraction of the population infected with the first disease becomes large enough to allow the spread of the second. No qualitative changes in the nature of the outbreaks appear. In [12] Chen, Ghanbarnejad, Cai, and Grassberger (CGCG) introduced a generalized SIR model (CGCG model) to include mutual cooperative effects of coinfections. In the CGCG model two different diseases simultaneously spread in a population. Having been infected with one disease gives an increased probability to get infected by the other. The amount of this increase is a proxy of the mutual cooperativity between the two diseases. The authors study the model at the mean-field level and observe that cooperative effects, depending of their strength, can cause a change of the transition from continuous to discontinuous. In [13] Janssen and Stenull showed that the CGCG model is equivalent, in the mean field, to the homogeneous limit of an extended general epidemic process and clarify the spinodal nature of the discontinuous transition observed. In [14,15] CGCG analyzed how in their model the type of transition depends on the contact network topology by simulating the CGCG model on both random networks and lattices. They concluded that a necessary condition for discontinuous transitions to occur, when starting from a doubly infected node, is the relative paucity of short loops with respect to long ones. More in detail they argue that a discontinuous transition occurs if the two epidemics first evolve separately and, only after the independent clusters of singly infected nodes have reached endemic proportions, the two epidemics meet: At that point cooperativity implies that both clusters rapidly become doubly infected. A necessary condition then is that there are few short loops (otherwise, the two pathogens immediately cooperate and the transition is continuous) and there are long loops (otherwise, cooperativity has no effect and one sees only single infections). In agreement with this scenario, CGCG do not observe discontinuous transitions on trees and on two-dimensional (2D) lattices, while they do observe them on Erdös-Rényi (ER) networks, on 4D lattices, and on 2D lattices with sufficiently long-range contacts. On 3D lattices the existence of discontinuous transitions depends on the details of the microscopic realization of the model. Moreover, the observed discontinuous transitions are of hybrid type [16,17], i.e., exhibiting also some features of continuous ones. For other recent work about cooperating infections, see Refs. [18-20].

An open question left by the analysis in [14,15] has to do with what happens for coinfections on generic broadly degree-distributed networks. Simulations performed on a Barabasi-Albert (BA) topology seem to indicate a continuous transition even for strong cooperativity. However, it is not clear whether these results are affected by finite-size effects and what happens for coinfections on other broadly distributed networks. Theoretical approaches [12,13] deal only with homogeneous networks.

In this paper we elucidate these issues. We first numerically show that for power-law-distributed networks of large size the transition is asymptotically discontinuous, if cooperativity is sufficiently high; however, strong size effects may conceal the real nature of the transition in finite systems, so that the transition appears continuous even for large networks. We then present an analytical heterogeneous mean-field approach to the problem, which allows us to derive a number of predictions, including the position of the threshold, the nature of the transition and the associated critical exponent. Numerical simulations validate the analytical predictions. The theoretical approach indicates that for scale-free networks, with diverging second moment of the degree distribution, the transition can only be continuous.

The paper is organized as follows. In Sec. II we describe the CGCG model in detail. Section III is devoted to a numerical investigation of the nature of the transition for power-law degree distributions, showing how finite-size effects can hinder the emergence of the discontinuity in scale-rich networks. Section IV is devoted instead to the theoretical analysis of the coinfection dynamics by means of a heterogeneous mean-field approach. Theoretical results are compared with numerical simulations, revealing a satisfactory agreement. We present conclusions and an outlook in Sec. V.

MODEL

In the classical SIR model in discrete time, each individual can be in one of three different states: susceptible (S, individuals that are healthy, neither infected nor immune), infected (I, individuals who can transmit the disease), or removed (R, dead or recovered and immunized individuals). At each time step each infected individual spontaneously decays with probability into the removed state, while she transmits the infection to each of her susceptible neighbors with probability .

The CGCG model is a modification of the SIR model with two circulating diseases, A and B. The infection probability for one disease is increased if the individual currently has, or has had in the past, the other disease: Individuals uninfected by either disease get infected (with either A or B) by any infective neighbor with probability , while a node that is or has been infected by one disease has a higher probability to get infected by the other pathogen. When recovering from one disease an individual becomes immune to that disease, but she can still be infected by the other. The model is totally symmetric with respect to A and B. Since each individual can be in one of three possible states (S, I, R) with respect to each of the two diseases (A, B), there are nine possible states for each individual, denoted by S, A, B, AB, a, b, aB, Ab, and ab, where, for each disease, capital letters refer to the infected state, while lowercase letters refer to the removed state. States denoted by single letters (a, b, A, B) refer to states where the individual is still susceptible with respect to the other disease.

We simulate the CGCG model on power-law distributed networks , generated according to the uncorrelated configuration model [21], for various values of the exponent. We fix the recovery probability to throughout the paper. In the various simulations, a number (ranging from 200 to 2000) of independent realizations are performed. We consider an initial condition with all individuals initially in the susceptible (S) state except for a single, randomly chosen individual who is in the doubly infected AB state.

FINITE-SIZE EFFECTS AND THE NATURE OF THE TRANSITION FOR HIGH COOPERATIVITY

We start by analyzing the behavior of the CGCG model on power-law degree-distributed networks. In Ref. [14] only continuous transitions were observed for the BA scale-free network. This is in contrast with the behavior of the same model on Erdös-Rényi networks, which exhibits a discontinuous transition for sufficiently high cooperativity. One may wonder whether the continuity of the transition for BA networks persists in the infinite size limit and whether the same occurs for other broadly distributed networks, which are, however, not scale-free.

The interpretation for the existence of a discontinuous transition put forward in Ref. [14] provides some hint on the issue. The origin of the discontinuity is traced back to the relative abundance of long loops compared to the scarcity of short loops in the topology. It is well known [22] that the clustering coefficient of uncorrelated random networks is Hence, the clustering coefficient vanishes for large size for any . However, this occurs differently depending on the value of . The factor in parentheses in Eq. (1) is finite for , but its value becomes larger as gets smaller. Hence, we do expect that for systems of fixed size the transition will be discontinuous for large , but the residual clustering will make the transition appear continuous as is reduced. This expectation is tested by simulating the CGCG dynamics on networks of fixed size (), large cooperativity (), and growing values of starting from (Fig. 1). As expected, upon increasing , the transition becomes discontinuous: Short loops become less abundant as grows so that the two infections evolve separately and give rise to cooperative effects only after they meet following one long loop in the network.

FIG. 1.

The six panels represent the total density of population recovered from both A and B in the final state plotted versus on power-law networks with , built by applying the uncorrelated configuration model. The system size is , and . The results are averaged over a number of independent realizations, that is, 200 for to 5, and is increased to 2000 for to 8.

Further evidence about the same effect is obtained by simulating the CGCG dynamics on a network with fixed and increasing system size from to (see Fig. 2): When the size is big enough, the discontinuity appears.

FIG. 2.

Total density of population recovered from both A and B in the final state, plotted versus on power-law networks of different sizes (, respectively), with . The parameter is 0.99. The results are averaged over independent realizations in panels (a), (b), (c), and (d), respectively.

The transition presented in Fig. 2 is discontinuous, but Fig. 3 shows that it is, in fact, of hybrid type, as discussed in Ref. [14]. This means that already at the transition point there may be giant infected clusters occupying a finite fraction of the system [see Figs. 3(a) and 3(c)]; this is at odds with what happens in the usual, continuous, percolation transition. However, the probability that one of these clusters appears undergoes a continuous transition at the threshold [Fig. 3(b)], as in normal percolation.

FIG. 3.

(a) Mass distribution of doubly infected clusters at the threshold . The number of realizations is at least . The solid straight line represents the decay . (b) Probability of forming a giant cluster (with size larger than ) as a function of and system sizes ranging (top to bottom) from to . The number of realizations is . (c) Plot of the size of the giant cluster as a function of for , showing a linear growth. All data are for and .

For values , Eq. (1) points out that, since diverges as , the decay of with the system size becomes , slower than . This qualitative change does not allow one to draw firm conclusions on the nature of the transition for . Asymptotically, the clustering coefficient vanishes also in this case, but it does so slowly so that it might be impossible for the two epidemics to grow initially separated clusters. Therefore, we cannot say whether the continuous transition observed in Ref. [14] is a finite-size effect or not. More on this issue will be provided by the analytical results in the next section.

HETEROGENEOUS MEAN-FIELD THEORY FOR UNCORRELATED NETWORKS

In this section we present an analytical approach to the cooperative dynamics on uncorrelated networks, based on the heterogeneous mean-field (HMF) theory [23,24]. It is a refinement of the homogeneous mean field used in Ref. [12], which allows for quantities to depend on the degree of the node considered and is thus suitable also for topologies with a broad connectivity distribution. For example, is the fraction of individuals in state A, restricted to nodes of degree .

Following the same formalism of [12], and exploiting the full symmetry among A and B pathogens, we define the quantities is the fraction of nodes of degree , which have contracted one of the infections and are susceptible with respect to the other; is the fraction of nodes of degree , which are currently infected by one of the pathogens, regardless of the state with respect to the other infection; is instead the fraction of nodes which have recovered from one infection, regardless of the state with respect to the other. Equations (4) of Ref. [12] can thus be generalized to

which conserve the total probability . The rates for the transmission of the first infection, of the second infection and for recovery are , and , respectively. is the probability that any given edge points to an infected node and is capable of transmitting the disease. For SIR dynamics one must take into account the fact that an infected individual cannot transmit the infection through the edge from which she was infected; hence [25], In the following we set the recovery rate , with no loss of generality. Notice that Eqs. (5) are for continuous dynamics, while the dynamics considered in simulations are discrete, with the lifetime of the infected state deterministically equal to 1. In order to compare results between theory and simulation, we must therefore consider the mapping between the rates and the probabilities and (see Appendix A). We also note that, since HMF considers probabilities averaged over all possible realizations of the dynamics, it cannot predict hybrid transitions.

Considering the initial condition with only one doubly infected node , and , Eqs. (5) are readily integrated, yielding

where the auxiliary function is defined as By deriving Eq. (8) with respect to time, we obtain

At the end of the spreading process and . Hence, , the asymptotic value of , obeys

Besides the zero solution , a nonzero solution representing endemic outbreaks can be obtained only if the derivative of the right-hand side of Eq. (12) with respect to , evaluated for , is larger than 1; i.e., Consequently, the epidemic threshold is Notice that the threshold does not depend on and is equal to the threshold for the SIR spreading of a single disease.

This conclusion is supported by the simulation results illustrated in Fig. 4. As the cooperativity grows, the nature of the transition passes from continuous to discontinuous; however, the position of the threshold remains unchanged. Notice that the independence of the threshold value from is another piece of evidence in favor of the interpretation of the transition put forward in [14]. If the discontinuous transition is originated by the encounter of macroscopic singly infected clusters, such an event is possible only when is larger than the single-infection threshold. This is what is observed in Fig. 4.

FIG. 4.

Large red symbols represent the final fraction of doubly infected nodes on power-law networks of size , with . They are plotted vs for various fixed values of the cooperativity . The small blue symbols are the results of the numerical integration of Eq. (5). The plotted values are computed as . The vertical black line is .

In the same figure we also plot the results of the numerical integration of the HMF Eqs. (5). The agreement between the theoretical and numerical results is only qualitative: We observe the same behavior (continuous or discontinuous) in both cases, but the analytical approach does not predict precisely the position of the threshold. This is not a surprise, as it is known that the HMF theory does not give accurate estimates of the SIR threshold for homogeneous networks (see Appendix A).

For single-infection SIR the critical value predicted by the HMF for continuous-time dynamics has exactly the same value of the exact threshold of discrete time dynamics. The same seems to occur for SIR coinfections: The threshold found in numerical simulations of the discrete dynamics is, again by coincidence, very close to the value predicted by HMF theory and given by Eq. (14).

The analytical approach allows also to derive the conditions under which the transition is continuous and the associated critical exponent , defined as

The detailed calculations, reported in Appendix B, give the following results, which provide a complete picture of the model behavior.

For the transition is continuous if the cooperativity is sufficiently small, i.e., and the associated exponent is . When instead the transition is discontinuous. This is in agreement with what is found in Fig. 4. In the marginal case , the transition is still continuous but with critical exponent for and for . This last result is consistent with the mean-field treatment in Ref. [12].

For the limit value for the transition to be continuous is always larger than 2 and assumes a nontrivial dependence on (see Fig. 5) that continuously interpolates between 2 for and infinity for . Also the critical exponent changes continuously . Also in this range of values things are different for marginal cooperativity : In this case regardless of the value of .

FIG. 5.

The value of the minimum cooperativity needed for a discontinuous transition, as predicted by HMF, as a function of .

For the threshold vanishes, the transition is continuous for any value of the cooperativity , and the growth of the order parameter for small is , where .

CONCLUSIONS

In this paper we have performed an analysis of the implications of cooperativity for SIR-like epidemic spreading on power-law degree-distributed networks. First, we have tested the interpretation put forward in Ref. [14] about the physical origin of the discontinuous epidemic transition occurring in this kind of system for large cooperativity. Via numerical simulations, we have shown that strong finite-system-size effects may hide the discontinuity, but that if the system size is large enough the true, hybrid, nature of the transition becomes evident. This provides a strong validation of the physical picture proposed by Cai et al. [14].

In the second part of the paper we have applied the heterogeneous mean-field theoretical approach to the coinfection dynamics on uncorrelated networks. In this way we have derived an expression for the epidemic threshold, which does not depend on cooperativity, in agreement with simulations. Moreover, the approach allows us to derive, for power-law degree-distributed topologies, the range of cooperativity values for which the transition is continuous and the associated critical exponent. We find a nice confirmation that the discontinuous transition occurs when cooperativity is sufficiently high and local clustering is instead small. An interesting point in this respect regards what happens for scale-free networks, i.e., for . In such a case HMF theory predicts that the transition is continuous for any value of the cooperativity . This is in agreement with the simulation results for BA networks presented in Ref. [14,15] for networks of moderate size. However, one cannot, in principle, exclude that also for the continuous transition observed in simulations is just a finite-size effect. While HMF theory seems to indicate the contrary, the global clustering coefficient Eq. (1) vanishes as the system size diverges also for , so that asymptotically there are no short loops in the system. Numerical simulations cannot be conclusive about this question. The issue could be clarified only by a careful analysis of the balance between the relative abundance of short loops (encoded by ) and long loops in the network. It remains as a challenging goal.

It is also worth noticing that our theoretical approach is fully general with respect to the value of the infection rates ratio . Therefore, our results provide predictions also for competing pathogens, for which . In that case our theory predicts always a continuous transition, for any perfectly analog to the case of weak cooperativity .

Mutual cooperation is a very important ingredient also for nonbiological spreading processes, such as the diffusion of ideas or the adoption of innovations. The investigation of cooperative effects in general complex contagion dynamics remains a very interesting avenue for future research activity.