Impact of the distribution of recovery rates on disease spreading in complex networks

PHYSICAL REVIEW RESEARCH 2, 013046 (2020) Impact of the distribution of recovery rates on disease spreading in complex networks Guilherme Ferraz de Arruda,1Giovanni Petri,1Francisco A. Rodrigues,2and Yamir Moreno3,1 1ISI Foundation, Via Chisola 5, 10126 Torino, Italy 2Departamento de Matemática Aplicada e Estatística, Instituto de Ciências Matemáticas e de Computação, Universidade de São Paulo, Campus de São Carlos, Caixa Postal 668, 13560-970 São Carlos, São Paulo, Brazil 3Institute for Biocomputation and Physics of Complex Systems and Department of Theoretical Physics, University of Zaragoza, 50018 Zaragoza, Spain (Received 14 February 2019; revised manuscript received 12 November 2019; published 14 January 2020) We study a general epidemic model with arbitrary recovery rate distributions. This simple deviation from the standard setup is sufficient to prove that heterogeneity in the dynamical parameters can be as important as the more studied structural heterogeneity. Our analytical solution is able to predict the shift in the critical properties induced by heterogeneous recovery rates. We find that the critical value of infectivity tends to be smaller than the one predicted by quenched mean-field approaches in the homogeneous case and that it can be linked to the variance of the recovery rates. Our findings also illustrate the role of dynamical-structural correlations, where we allow a power-law network to dynamically behave as a homogeneous structure by an appropriate tuning of its recovery rates. Overall, our results demonstrate that heterogeneity in the recovery rates, eventually in all dynamical parameters, is as important as the structural heterogeneity. DOI: 10.1103/PhysRevResearch.2.013046 I. INTRODUCTION Heterogeneity, whether in the nature of the components or the pattern of connections, is a key characteristic of complex systems. This is particularly evident in the case of the spreading of a disease in a networked population, where the inclusion of structural heterogeneity has long been known to radically change the process’s critical behavior [1–6]. As an illustration, consider two classical contagion models, the susceptible-infected-susceptible (SIS) and the susceptibleinfected-recovered (SIR) models. On a homogeneous network, they both present a nonvanishing critical point [5,6]. However, the introduction of structural heterogeneity, in the form of broad degree distributions of the nodes, can result in a vanishing critical point [1,5–7]. More specifically, in the thermodynamic limit, a divergence of the second moment of the degree distribution [1,5,6,8] or a divergence in the maximum degree [1,5,6,8] implies a vanishing critical infectivity. This in turn has important practical implications for real-world networks, because many of them display very broad [7–9] (or even scale-free) degree distributions [10]. References [11–18] considered, in contrast, homogeneous structures but accounted for arbitrary times in the state transitions. Interestingly, previous works studied models characterized by a series of infected states (without a biological interpretation) that described the global behavior accurately [12–15,17]. Notably, the model proposed in [12] is general, but each transition has Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. an assigned function, which complicates the model, making its practical use limited. These works encounter, however, a fundamental limitation, because they describe the population at the mean-field level [11–18]. While structural and dynamical heterogeneities are independently accounted for, modeling both types of heterogeneities together has received considerably less attention until recently. Indeed, it was mainly studied for the SIR model: A message-passing formalism was proposed in [19,20] and a heterogeneous mean-field approach in [21]. In the latter, the authors also performed numerical experiments showing that the population can be more vulnerable in the scenario with dynamical heterogeneity. More recently, this problem was investigated on temporal networks [22] using an SIS process, which was described within the quenched mean-field (QMF) formalism and mainly focusing on spreading rates [23,24]. A similar approach was considered in [25], where however the authors highlighted a different aspect, i.e., the allocation of resources during an outbreak. Here we investigate a different type of dynamical heterogeneity by characterizing the critical properties of an SIS model when recovery rates are distributed heterogeneously across the population. Heterogeneous recovery rate distributions can be associated with biological differences between individuals [26,27], demographic characteristics [28], and social differences that result in nonhomogeneous access to the health system [29]. Additionally, we also consider the case in which correlations arise between the structure and dynamics. We show, both analytically and numerically, that such correlations can induce opposing and unexpected dynamical outcomes, for example, power-law (PL) networks displaying nonvanishing critical points and conversely homogeneous networks displaying vanishing critical points. Indeed, we show that the standard QMF predictions, which are a lower bound in the standard scenario, no longer provide such a bound in the 2643-1564/2020/2(1)/013046(8) 013046-1 Published by the American Physical Society DE ARRUDA, PETRI, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013046 (2020) heterogeneous rate case. We then propose a simple formulation to overcome this and provide a different lower bound for the process. Our results complement previous evidence on the SIR model [21] and imply that proper characterization of the dynamical parameters is of utmost importance not only for a better understanding of spreading processes, but also for many practical applications, such as surveillance, forecasting, resource management, and network reconstruction, among many others. II. THE SIS MODEL WITH HETEROGENEOUS RECOVERIES We start by considering a population composed of N individuals with an arbitrary pattern of connections in a single connected component. These connections can be represented as a network and are described by the network’s (usually symmetric) adjacency matrix A. Each individual can be in one of two states: (i) infected (Yi=1) or (ii) susceptible (Xi=1). Using a Markovian approach, the epidemic process is modeled as a collection of independent Poisson processes. In order to model the spreading of the disease through the network of contacts, with each directed edge i∼j, emanating from the infected individual i, we associate a Poisson process with rate λij,Nλij(t) (whose transitions are Yi+Xj→Yi+ Yj). Additionally, with each infected individual, we associate a Poisson process with rate δi,Nδi(t), modeling the recovery (Yi→Xi). This system is statistically described using the order parameter ρand the susceptibility χ, defined as ρ=1 N N  i Yi,χ=nI2−nI2 nI,(1) where nIis the number of infected individuals. Both quantities can be directly estimated using Monte Carlo methods, in particular, the quasistationary method and the Gillespie algorithm, where each of the processes mentioned above is simulated, and the state of the nodes is evaluated [6]. In the QMF approach, one implicitly assumes that XiYj≈ XiYj. Physically, this corresponds to neglecting dynamical correlations. Defining yi=Yi, the process is described as dyi dt =−δiyi+(1 −yi) j λijAijyj.(2) Thus, defining ii =δi,W=λij, and Q=−1(A◦W), the critical point is given as λQMF c=[max(Q)]−1,(3) where max(Q) is the leading eigenvalue of Q. Note that the elements of Qare the expected number of contacts before recovery. Obviously, the critical point simplifies to τQMF,std c= (λ δ)c=[max(A)]−1in the homogeneous case, i.e., when δi= δand λij =λ. Thus, in the thermodynamic limit, the critical point for PL networks goes to zero if the maximum degree is a growing function of the network size [5,6,8]. On the other hand, we can consider a scenario that is equivalent to the contact process (CP) by setting the spreading rates as λij =λ ki, which is thus described by the probability transition matrix Pij =Aij ki. In this case, the critical point is finite and 0.00 0.02 0.04 0.06 0.08 0.10 λ 10− 1 100 101 102 χ 2 3 4 5 10 100 1000 α 100101102 α 10−3 10−2 10− 1 λQMF c 0.0 0.1 λc 0.0 0.1 λQMF c FIG. 1. Susceptibility curves extracted from the Monte Carlo simulations for an Erd˝ os-Rényi network with N=105and k≈10 considering that the recovery rate distribution follows an inversegamma distribution, whose shape parameter αis color coded. The insets show the dependence of λQMF con α(top) and a comparison with numerical results (bottom), also for different α. τQMF,CP c=1, regardless of the underlying structure. Note that, similarly to the homogeneous case, this prediction [Eq. (3)] is an upper bound for the heterogeneous recovery rate scenario, because it relies on the independence of the random variables: If i∼j, then P(Yi=1|Yj=1) ⩾P(Yi=1) =yi, which implies that the nodal probability is always overestimated (see [2] for a similar argument). From here onward, we set λij =λ and focus on the effect of the recovery rate distribution on the critical point. Considering an undirected network, from the matrix norm we can bound Eq. (3) using the standard QMF predictions as min(δi)τQMF,std c⩽λQMF c⩽max(δi)τQMF,std c.(4) We can see that Eq. (4) suggests that the standard QMF predictions might not be a lower bound to the alternative process. Note that the uncertainty assuming the standard QMF prediction increases as the variance also increases. III. SYNTHETIC NETWORKS To further characterize the critical behavior of our model, we first consider an Erd˝ os-Rényi (ER) network with N=105 and k≈10 (therefore τQMF,std c≈0.1). This graph has a homogeneous structure and allows us to analyze the structural and dynamical effects independently. Mounting evidence shows that infectious times in real epidemics follow a gamma distribution [17,30,31]. Consequently, the rate distribution must follow an inverse-gamma distribution. We impose therefore the recovery rates to have an inverse-gamma distribution δ∼−1(α,β), where αand βare the shape and scale parameters, respectively. Its mean is δi= β α−1and its variance is Var(δi)=β2 (α−1)2(α−2) for α>2. To facilitate the comparison between different distributions, we restrict the distributions to unitary mean. In Fig. 1we present the critical behavior of an ER network for different shapes α. The top inset emphasizes the behavior of the predicted critical point as a function of αand its comparison with estimations from Monte Carlo 013046-2 IMPACT OF THE DISTRIBUTION OF RECOVERY RATES … PHYSICAL REVIEW RESEARCH 2, 013046 (2020) simulations. As expected, for sufficiently large values of α, the dynamics behave similarly to the standard SIS model with uniform δ, where the predicted threshold coincides (dashed line in the top inset). The agreement between analytical and simulated critical points is excellent, as it can be seen in the bottom inset of Fig. 1. Furthermore, as αdecreases, the variance of δand consequently the uncertainty of the bounds are also enlarged. However, we observed that the critical point systematically moved towards zero in our simulations, which is also consistent with Eq. (4). Note that the inverse-gamma distribution is asymmetric with respect to its average and centered at δi⩽δi. This implies that the bounds in Eq. (4)are also asymmetric. Since the recovery rates are sampled from this distribution, we expect that the number of individuals that take longer to recover is greater than the number of individuals that recover fast. This suggests that, via infection/reinfection mechanisms, the disease can survive for lower values of λ as compared to the standard QMF predictions. In addition to the critical properties, we also show a different supercritical behavior; for example, the low-αregime of a network can be similar to, and even in some regions it can be mistaken for, the high-αregime on top of a network with a different structure. This multiplicity of supercritical behaviors also raises questions in network reconstruction models based on disease dynamics. IV. REAL-WORLD NETWORKS We confirmed that our results hold in real-world networks. As before, we consider an inverse-gamma distribution. Figure 2shows results of simulations in two real networks: the UC Irvine messages social network [Fig. 2(a)][32,33] and the OpenFlights network [Fig. 2(b)][33,34]. These networks represent different spatial scales of a similar spreading process: The social network corresponds to smaller scales and spatially localized systems, while the OpenFlights network captures a wider spatial scale. In the top inset of Fig. 2(b) we show that the critical point predictions are remarkably good for inverse-gamma recovery rates, even for these real networks. Figures 1and 2reveal that the critical point decreases as we increase the variance of the recovery rate distribution. Thus, if we consider heterogeneous networks, assuming an average recovery rate in the QMF is not enough to provide an adequate characterization of the process. It is not a lower bound anymore, as can be seen in Fig. 1, where the numerically estimated critical point is always slightly lower than the value predicted by the standard QMF predictions (τQMF,std c≈0.1). The proper correction for the QMF predictions is given by Eq. (3), which is a lower bound for the underlying process (see the bottom inset in Fig. 1). V. EFFECTS OF DYNAMICS-STRUCTURE CORRELATIONS The bounds in Eq. (4) implicitly assume that there are no correlations between the structure and dynamics. From the Gershgorin circle theorem we know that every eigenvalue of Qlies at least in one of the disks D(Qii,Ri) centered in Qii with a radius given as Ri=i= j|Qij|. Therefore, considering a symmetric matrix, |k|⩽Qii +Rk, hence 0.000 0.005 0.010 0.015 0.020 0.025 0.030 λ 10−1 100 101 χ (a) 0.000 0.005 0.010 0.015 0.020 λ 10−1 100 101 χ (b) 0.00 0.02 λ 0.0 0.1 0.2 ρ 0.00 0.02 λ 0.00 0.05 0.10 ρ 2 3 4 5 10 100 α 0.01 0.02 λc 0.01 0.02 λQMF c FIG. 2. Results for the susceptibility when the dynamics are on top of real networks. We use (a) the UC Irvine messages social network [32,33](•) and (b) the OpenFlights network [33,34](). In both cases, we considered the undirected version of the giant component. In the main figure of each panel, we present the susceptibility for different values of αand λ. In the bottom right insets, we present the order parameter. In the top inset of (b), we present a comparison between the QMF estimated and predicted critical points. max ⩽Q∞, where the infinity norm is defined as Q∞=max 1⩽i⩽N⎛ ⎝ N  j=1 Aij δi⎞ ⎠=max 1⩽i⩽Nki δi.(5) If the structure and the dynamics are correlated, Eq. (5) might give us further insight. For instance, for the PL case, the leading eigenvalue of Adiverges in the thermodynamic limit, leading to a vanishing critical point. Conversely, using Eq. (5) and a proper choice of δi, we can change this behavior. In fact, assuming that δi(ki)∝kiin the thermodynamic limit, we have lim N→∞ Q∞=lim N→∞ max 1⩽i⩽Nki δi=c,(6) where c<∞is a finite real constant. This radically changes the critical behavior of the dynamics. Note that both the CP (λij =λ ki) and the δi=kicases are described, at first order, by the probability transition matrix P, yielding τQMF,CP c= λQMF c=1. 013046-3 DE ARRUDA, PETRI, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013046 (2020) 1.0 1.2 1.4 1.6 1.8 λ 100 101 102 χ (a) γ=3.5γ=2.7γ=2.1 N=10 3 N=10 4 N=10 5 N=10 6 0.1 0.2 0.3 0.4 0.5 λ 10−1 100 101 102 χ (b) γ=3.5γ=2.7γ=2.1 103104105106 N 0.2 0.3 0.4 λc FIG. 3. Finite-size analysis considering structure-dynamics correlations. We present the susceptibility (colors represent the sizes) in (a) the heterogeneous recovery rates, where δi=kion top of the power-law networks, whose exponents are γ=2.1,2.7,3.5, and (b) the SIS model with heterogeneous recovery rates, considering an Erd˝ os-Rényi network and the recovery rates as δi=ki kPL ,wherekPL is a discrete power law, more specifically, the degree sequence of the networks used in (a). In the inset we present the critical point as a function of the system size in log-log scale. In Fig. 3(a) we perform a finite-size analysis when the dynamics occurs on top of PL networks and recovery rates are δi=ki, finding evidence of a finite critical point, i.e., the value of λcorresponding to the peak of χdoes not vanish. For comparison, the results for the CP are reported in Appendix E. We remark that the convergence towards the critical value for growing Nseems to be slower in the CP case. In summary, our results show that a network with a power-law degree distribution may show a finite critical point for the SIS dynamics if degrees and recovery rates are appropriately correlated. Next we show that the inverse scenario is also possible. Consider an ER network with k≈10 and δi=ki kPL , where kPL is a discrete power law, i.e., P(kPL)∼k−γ PL . That is, we now have a homogeneous structure and a heterogeneous recovery rate distribution. In Fig. 3(b) we show a finite-size analysis for this configuration with varying γ=2.1,2.7,3.5. We observe that for γ=2.7 and 3.5 the underlying structure plays an important role, suggesting a nonvanishing critical point or considerable reduction in the scaling exponent [see Fig. 3(b) inset, where both curves have a slope close to zero]. However, for γ=2.1 our results indicate the existence of a vanishing critical point [see Fig. 3(b) inset]. It seems reasonable to hypothesize that the scenario observed when γ=2.1 is due to the fact that, in the steady state, the infection probabilities are inversely proportional to the nodal recovery rates and thus that the evaluation of the recovery time at both ends of every edge enables an infection-reinfection mechanism. Understanding what the necessary and sufficient conditions to observe such a mechanism are and the interplay between structure and dynamics requires, however, further exploration. In [8,35,36] it was shown that, depending on the network structure, the transition can be triggered by different activation mechanisms, namely, (i) collective (e.g., ER networks or the CP), (ii) kcore (uncorrelated power law with 2 <γ <2.5), or (iii) hub (uncorrelated power law with γ>2.5). Note that for uncorrelated power laws, this coincides with the different regimes of leading eigenvector localization [37][forthe analysis of the inverse participation ratio (IPR) of the leading eigenvector of Qwe refer the reader to Appendix F]. Thus, we conjecture that it may be possible to alter these structural mechanisms with a proper recovery rate distribution, possibly considering dynamics-structure correlations. Note that we implicitly showed [Fig. 3(a)] how to transform k-core (γ=2.1) and hub-activated (γ=2.7,3.5) mechanisms into a collective behavior phenomenology. Importantly, our findings may lead to alternative prevention/intervention techniques that take advantage of the phenomenology reported here. We also highlight that the consequences of our findings are not limited to the cases explored here. For instance, it is natural to conjecture the existence of Griffith’s phase in our setup. In this type of transition, we have an extended critical region instead of a single critical point, which was studied in complex networks in [36,38–42]. More specifically, in [39] the authors showed that slow dynamics on a weighed treelike structure can occur in a contact process. Its similarities to our scenario thus suggest that similar phenomenology is also possible in our case. VI. CONCLUSION We have analyzed the impact of introducing heterogeneity in the recovery rates of an SIS disease dynamics. We showed that dynamical heterogeneity is as significant as structural heterogeneity and that it can induce drastic changes in the SIS critical properties. Furthermore, our results show that the standard QMF approach does not provide a lower bound for the heterogeneous case anymore. To solve this inconsistency, we proposed a solution that relates the structural and dynamical features via the spectral properties of a different matrix Q. This formulation presents opportunities for future research. For example, our findings raise questions about the consequences of heterogeneities in spreading rates and the interplay between spreading and recovery rates. Furthermore, with regard to control/containment strategies, heterogeneity in the recovery rates can be considered for intervention strategies (or even immunization, δi→0). In contrast, heterogeneity of spreading rates would be related to prevention. In this context, 013046-4 IMPACT OF THE DISTRIBUTION OF RECOVERY RATES … PHYSICAL REVIEW RESEARCH 2, 013046 (2020) our findings might also be of potential interest. In addition to the specific conclusions drawn here, there are others that concern more general aspects of disease spreading processes as well as the characterization of complex systems in general. For instance, differences in the localization properties induced by dynamical heterogeneities might influence the predictability of complex systems, and in particular of diseases [43], or in the reconstruction of networks from the dynamics as proposed in [44]. Finally, we also stress that our results might also have an impact on information spreading processes since individuals have different activity timescales, which could ultimately be related to the recovery time distribution. ACKNOWLEDGMENTS G.F.d.A., G.P., and Y.M. acknowledge partial support from Intesa Sanpaolo Innovation Center. Y.M. acknowledges partial support from the Government of Aragón, Spain through Grant No. E36-17R and from MINECO and FEDER funds (Grant No. FIS2017-87519-P). Research was carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry funded by FAPESP (Grant No. 2013/07375-0). APPENDIX A: EQUATION (4) AS A BOUND FOR EQ. (3) The critical point given by Eq. (3) can be bounded using a matrix norm. Thus, using the 2-norm, we obtain A2 max()⩽−1A2⩽A2 min(),(A1) where A2=max(ATA). Equation (A1) therefore provides bounds on the leading eigenvalue of Qgiven by the structure and the variance of δi. Next, for undirected networks we have that A2=max(A)=(τQMF c)−1. Moreover, the Perron-Frobenius theorem for non-negative matrices states that every non-negative matrix can be written as a limit of positive matrices, and therefore one has the existence of an eigenvector with non-negative components and the corresponding eigenvalue will be non-negative and greater than or equal (in absolute value) to all other eigenvalues [45]. Therefore, from the norm definition, Q2=σmax(Q)=max i[i(Q)] (A2) since Qis non-negative. Finally, inverting the equations and writing it in terms of the critical properties, we obtain Eq. (4). APPENDIX B: ANALYZING EQ. (4) FOR THE INVERSE-GAMMA DISTRIBUTION The bounds obtained in Eq. (4)or(A1) cannot be further improved without knowledge about the relationship between the matrix Aor Wand −1[for more details about the effects of dynamics-structure correlations, see Eq. (5)]. Alternatively, one can use the asymmetries in the recovery rate distribution to understand how the lower and upper bounds behave. More specifically, the cumulative distribution for the inverse-gamma 101102103104 0.3 0.4 0.5 0.6 0.7 α δ≤ P ( 1) δ > P ( 1) P FIG. 4. Evaluation of P(δ⩽1) and P(δ>1) as a function of α, the shape parameter of an inverse-gamma distribution with β= α−1. distribution is given as P(δ⩽x)=α, β x (α)=Qα, β x,(B1) where (a,b) is the upper incomplete Gamma function, (a) is the Gamma function, and Q(a,b) is the regularized gamma function. Next, keeping the average fixed, i.e., setting β= α−1, we can analyze how these bounds behave as a function of the shape α. We can calculate the probability that δis larger or smaller than a give value x. Formally, P(δ⩽x)=Qα, α−1 x,(B2) P(δ>x)=1−Qα, α−1 x.(B3) Next, evaluating (B2) and (B3) for the average value δi= x=1, we obtain the fraction of individuals that have a recovery rate that is lower and higher, respectively, than average. Since Q(α, α−1 x)⩾1 2the inverse-gamma distribution is peaked below the average value, therefore suggesting that the lower bound decreases faster than the upper bound. Interestingly, α→∞implies that P(δ⩽1) =P(δ>1) =1 2, and thus the uncertainty is symmetrical. This is also depicted in Fig. 4. We remark that the discussion of this Appendix concerns the uncertainty of the bounds and their asymmetry, which only suggests (but does not prove) that the critical point might decrease. This hypothesis is supported by our simulations. APPENDIX C: UNIFORM SYMMETRIC RECOVERY RATE DISTRIBUTION In the main text we focused on the inverse-gamma distribution. This distribution is asymmetric. In this Appendix we focus on a symmetric scenario considering the uniform distribution U(a,2−a) whose average is δi=1 and variance is Var(δi)=(1−a)2 3. We remark that a normal distribution was not considered since it is defined for negative numbers, while rates are positive. Figure 5shows the critical point estimations 013046-5 DE ARRUDA, PETRI, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013046 (2020) 0.00 0.02 0.04 0.06 0.08 0.10 λ 10− 1 100 101 102 χ 0.01 0.1 0.2 0.4 0.6 0.8 1.0 0.0 0.1 λc 0.0 0.1 λQMF c a FIG. 5. Susceptibility curves extracted from the Monte Carlo simulations for an Erd˝ os-Rényi network with N=105and k≈10 considering that the rate distribution follows a uniform distribution (symmetric) U(a,2−a), where the parameter acontrols the variance and is color coded. in this scenario. The same conclusions we obtained for Fig. 1 also apply here. The same trend is observed, in which the critical point moves to the left. The main observed difference is that, in this case, the critical predictions are poorer when the variance is larger. Although the error is larger in these cases, our predictions are still better than the standard QMF formulation. APPENDIX D: SUPERCRITICAL BEHAVIOR Aside from the critical behavior, discussed in the main text, it is also instructive to evaluate the supercritical behavior of our dynamics. Figure 6shows the phase diagram for Erd˝ osRényi and PL networks with γ=2.1,2.7,3.5. It is clear that not only the critical but also the supercritical behavior changes for different values of α. Note that, in some cases, the lower-αregime is similar to the higher-αregime of a different structure. This is particularly clear for the PL case of γ=3.5 and α=2.0, which is very similar to γ=2.7 and α→∞. A similar effect is also visible for other parameters. This is a particularly important issue if one uses the dynamical process to infer the network structure [44]. Furthermore, for sufficiently large values of λall the curves tend to converge; however, this convergence is arguably slower for PL networks, as emphasized in the inset of Fig. 6. APPENDIX E: CONTACT PROCESS In the continuous formulation, the CP is modeled as a set of Poisson processes, where each edge has the spreading rate defined as λij =λ ki. The recovery rate is kept fixed for all nodes. In this dynamical process, the average number of contacts of each node is the same for all individuals in the network regardless of their degree. This process is described by the probability transition matrix P, hence τQMF,CP c=1. In Fig. 7we can observe that the critical point predictions are not as accurate as for the previous case we studied in the main text [Fig. 3(a)]. This was anticipated in [46,47]. Interestingly, a better prediction for the critical point in the CP is obtained using the heterogeneous mean field, which (a) (b) FIG. 6. Phase diagram for (a) Erd˝ os-Rényi (•)networksand (b) PL networks. In (b), from top to bottom, we have γ=2.1(), γ=2.7(), and γ=3.5(). In the inset in (a) we present the behavior for larger λ’s in all the networks. They comprise N=104 nodes. The recovery rate distribution is given as an inverse-gamma distribution, whose shape parameter αis color coded. predicts τHMF, CP c=k k−1[47]. Thus, the mismatch between prediction and estimated critical points, in both cases, seems to be related to dynamical correlations. Due to its similarities, it is instructive to compare this model with our dynamic-structure correlated model, i.e., δi= ki.InFig.7we show the finite-size analysis for the CP case on top of power-law networks with exponents γ=2.1,2.7,3.5. We remark that this process has a finite critical point. Thus, comparing with Fig. 3(a), we observe that both processes present a similar behavior as a function on N. However, note that both models have slightly different estimated critical points. Aside from that, our results also suggest that the convergence to the asymptotic value seems to be faster in the CP case. We can thus conclude that although both processes have the same predictions for the critical point, the secondorder effects are different in both of them. APPENDIX F: LOCALIZATION PROPERTIES On top of some structures and just above the critical point, the mean-field theory predicts that the disease can be restricted (localized) in a subgraph [3,48]. This phenomenon is known as metastable localization. In the standard scenario, this property depends only on the structure, more specifically, 013046-6 IMPACT OF THE DISTRIBUTION OF RECOVERY RATES … PHYSICAL REVIEW RESEARCH 2, 013046 (2020) 100 101 102 χ γ=3.5γ=2.7γ=2.1 (a) 1.0 1.2 1.4 1.6 1.8 λ 10−5 10−4 10−3 10−2 10−1 ρ (b) N=10 3 N=10 4 N=10 5 N=10 6 FIG. 7. Finite-size analysis for the CP for (a) susceptibility and (b) the order parameter. The colors represent the sizes and from bottom to top, the curves are grouped by the power-law exponents γ=2.1,2.7,3.5, respectively. on the leading eigenvector of A. This property is very robust with respect to the spreading mechanism. In [48] the authors showed that bursts of infection are not sufficient to turn a localized spreading into a delocalized spreading. Interestingly, in our model, we can show that the recovery rates are enough to alter the localization. We remark that the control of the localization of diseases is still an open and interesting problem, as also pointed out in [48]. As anticipated, our model plays an important role in the localization of the leading eigenvector. More specifically, the matrix Qmight present a different localization pattern if compared to the adjacency matrix. To quantify localization we use the so-called inverse participation ratio, initially used in epidemics in [3]. Formally, it is given as IPR = N  i v4 i,(F1) FIG. 8. The IPR versus Nfor the Erd˝ os-Rényi, PL networks and the correlated case −1A, where i=ki kPL for (a) γ=2.1and (b) γ=3.5. Fifty networks were considered in each case. where vis the normalized leading eigenvector, i.e., v= 1. As shown in [37] and recently formalized in [48], the eigenvector is localized in a subextensive portion, i.e., IPR ∼ O(N−ν) with 0 <ν<1. Thus, in the fully delocalized case, IPR ∼O(N−1), and in the fully localized case, IPR ∼O(1). As expected, Fig. 8shows that in an ER network the IPR scales with ν≈1, ν<1 for a PL network with γ=2.1 [Fig. 8(a)], and a PL network with γ=3.5 is localized [Fig. 8(b)]. Complementarily, when δi=ki, the leading eigenvector of P=Qis homogeneously distributed; therefore IPR ∼O(N−1) and the disease is fully delocalized, which suggests that this is an effective strategy for full delocalization. Furthermore, Fig. 8also shows that when heterogeneous recovery rates are considered, the IPR might also change. It is possible that a network that exhibits full delocalization in the standard case (light pink circles for an ER network in Fig. 8) will instead exhibit localization when the recovery rates and structure are correlated (blue triangles in Fig. 8). Thus, by introducing localization and/or tuning, its scaling exponent νis not trivial. More specifically, one might naively use the recovery rates as δi=ki kPL as a strategy. In this case, kiwould delocalize the matrix, while kPL would control the new IPR towards the distribution kPL. To evaluate this strategy, we fixed the structure to an ER network. Next we defined kPL as a PL degree distribution and opted for a discrete approach for kPL, which allowed for a comparison with the localization properties of a network with the same degree distribution, i.e., kPL. However, as we observed in Fig. 8, this was insufficient to control localization. In the evaluated cases, the resulting matrix Qis localized. Note that kPL ∼k−2.1is partially delocalized [dark pink squares in Fig. 8(a)], but its composition with the delocalized structure results in a localized IPR. The case of kPL ∼k−3.5is less surprising, but also emphasizes the role of heterogeneity in the recovery times in disease localization. Therefore, with this simple strategy, we could not control the localization for an arbitrary ν. [1] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86,3200 (2001). [2] P. V. Mieghem, J. Omic, and R. Kooij, IEEE/ACM Trans. Netw. 17,1(2009). [3] A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F. Mendes, Phys. Rev. Lett. 109,128702 (2012). [4] G. F. de Arruda, E. Cozzo, T. P. Peixoto, F. A. Rodrigues, and Y. Moreno, Phys. Rev. X 7,011014 (2017). 013046-7 DE ARRUDA, PETRI, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013046 (2020) [5] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys. 87,925 (2015). [6] G. F. de Arruda, F. A. Rodrigues, and Y. Moreno, Phys. Rep. 756,1(2018). [7] S. Chatterjee and R. Durrett, Ann. Probab. 37,2332 (2009). [8] M. Boguñá, C. Castellano, and R. Pastor-Satorras, Phys. Rev. Lett. 111,068701 (2013). [9] A. Clauset, C. Shalizi, and M. Newman, SIAM Rev. 51,661 (2009). [10] I. Voitalov, P. van der Hoorn, R. van der Hofstad, and D. Krioukov, Phys. Rev. Res. 1,033034 (2019). [11] H. W. Hethcote, H. W. Stech, and P. van den Driessche, J. Math. Biol. 13,185 (1981). [12] Z. Feng and H. R. Thieme, SIAM J. Appl. Math. 61,803 (2000). [13] Z. Feng and H. R. Thieme, SIAM J. Appl. Math. 61,983 (2000). [14] A. L. Lloyd, Proc.R.Soc.B268,985 (2001). [15] A. L. Lloyd, Theor. Popul. Biol. 60,59 (2001). [16] Z. Feng, D. Xu, and H. Zhao, Bull. Math. Biol. 69,1511 (2007). [17] O. Krylova and D. J. D. Earn, J. R. Soc. Interface 10,20130098 (2013). [18] D. Clancy, Stat. Probab. Lett. 85,1(2014). [19] B. Karrer and M. E. J. Newman, Phys. Rev. E 82,016101 (2010). [20] N. Sherborne, J. C. Miller, K. B. Blyuss, and I. Z. Kiss, J. Math. Biol. 76,755 (2018). [21] W. Gou and Z. Jin, Infect. Dis. Model. 2,353 (2017). [22] A. Darbon, D. Colombi, E. Valdano, L. Savini, A. Giovannini, and V. Colizza, R. Soc. Open Sci. 6,181404 (2019). [23] B. Qu and H. Wang, IEEE Trans. Netw. Sci. Eng. 4,177 (2017). [24] B. Qu and H. Wang, Physica A 472,13 (2017). [25] S. Ottaviano, F. De Pellegrini, S. Bonaccorsi, and P. Van Mieghem, J. Complex Netw. 6,800 (2017). [26] S. Segal and A. V. Hill, Trends Microbiol. 11,445 (2003). [27] H. R. Fryer, J. Frater, A. Duda, M. G. Roberts, The SPARTAC Trial Investigators, R. E. Phillips, and A. R. McLean, PLOS Pathogens 6,e1001196 (2010). [28] S. Dorjee, Z. Poljak, C. W. Revie, J. Bridgland, B. McNab, E. Leger, and J. Sanchez, Zoonoses Public Health 60,383 (2013). [29] R. M. Barber, N. Fullman, R. J. D. Sorensen, T. Bollyky, M. McKee, E. Nolte, A. A. Abajobir, K. H. Abate, C. Abbafati, K. M. Abbas et al.,Lancet 390,231 (2017). [30] M. Eichner and K. Dietz, Am. J. Epidemiol. 158,110 (2003). [31] H. J. Wearing, P. Rohani, and M. J. Keeling, PLOS Med. 2,e174 (2005). [32] T. Opsahl and P. Panzarasa, Social Netw. 31,155 (2009). [33] J. Kunegis, Proceedings of the 22nd International Conference on World Wide Web, Rio de Janeiro, 2013 (ACM, New York, 2013), pp. 1343–1350. [34] OpenFlights network data set, KONECT, 2016. [35] C. Castellano and R. Pastor-Satorras, Sci. Rep. 2,371 (2012). [36] W. Cota, A. S. Mata, and S. C. Ferreira, Phys.Rev.E98,012310 (2018). [37] R. Pastor-Satorras and C. Castellano, Sci. Rep. 6,18847 (2016). [38] M. A. Muñoz, R. Juhász, C. Castellano, and G. Ódor, Phys. Rev. Lett. 105,128701 (2010). [39] G. Ódor and R. Pastor-Satorras, Phys.Rev.E86,026117 (2012). [40] R. Juhász, G. Ódor, C. Castellano, and M. A. Muñoz, Phys. Rev. E85,066125 (2012). [41] H. K. Lee, P.-S. Shim, and J. D. Noh, Phys. Rev. E 87,062812 (2013). [42] G. Ódor, R. Dickman, and G. Ódor, Sci. Rep. 5,14451 (2015). [43] S. V. Scarpino and G. Petri, Nat. Commun. 10,898 (2019). [44] T. P. Peixoto, Phys. Rev. Lett. 123,128301 (2019). [45] Matrix Analysis and Applied Linear Algebra,editedby C. D. Meyer (Society for Industrial and Applied Mathematics, Philadelphia, 2000). [46] S. C. Ferreira, R. S. Ferreira, and R. Pastor-Satorras, Phys. Rev. E83,066113 (2011). [47] A. S. Mata, R. S. Ferreira, and S. C. Ferreira, New J. Phys. 16, 053006 (2014). [48] Q. Liu and P. Van Mieghem, IEEE Trans. Netw. Sci. Eng. 6,983 (2019). 013046-8