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

TL;DR

This paper analyzes how heterogeneity in recovery rates influences epidemic spreading on complex networks, revealing shifts in critical thresholds and the interplay between network structure and dynamical parameters.

Contribution

It introduces an analytical model incorporating arbitrary recovery rate distributions, highlighting their impact on epidemic thresholds and network behavior.

Findings

Heterogeneous recovery rates significantly alter epidemic critical points.

Power-law networks can mimic homogeneous behavior through recovery rate tuning.

Recovery rate heterogeneity affects spreading localization properties.

Abstract

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. Additionally, we show 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. We then illustrate the role of dynamical--structural correlations, which allow for a complete change in the critical behavior. We show that it is possible for a power-law network topology to behave similarly to a homogeneous structure by an appropriate tuning of its recovery rates, and vice versa. Finally, we show how heterogeneity in recovery rates affects the network localization properties of the spreading process.

Equations

$$\rho =\frac{1}{N}\sum _i^N⟨Y_i⟩,\chi =\frac{⟨{n_I}^2⟩-{⟨n_I⟩}^2}{⟨n_I⟩},$$

$$\frac{d{y}_i}{dt}=-{\delta }_i{y}_i+(1-y_i)\sum _j{\lambda }_{ij}{\mathbf{A}}_{ij}{y}_j.$$

$$\frac{d\mathbf{y}}{dt}=\left(\lambda \mathbf{A}\circ \mathbf{W}-\Delta \right)\mathbf{y},$$

$$\mathbf{y}=\lambda {\Delta }^{-1}\left(\mathbf{A}\circ \mathbf{W}\right)\mathbf{y},$$

$${\lambda }_c={\left({\Lambda }_{max}(\mathbf{Q})\right)}^{-1}.$$

$$\frac{\Vert \mathbf{A}{\Vert }_2}{{\Lambda }_{\max }(\Delta )}\le \Vert {\Delta }^{-1}\mathbf{A}{\Vert }_2\le \frac{\Vert \mathbf{A}{\Vert }_2}{{\Lambda }_{\min }(\Delta )},$$

$$\Vert \mathbf{Q}{\Vert }_{\infty }=\underset{1\le i\le N}{\max }\left(\sum _{j=1}^N\frac{{\mathbf{A}}_{ij}}{{\delta }_i}\right)=\underset{1\le i\le N}{\max }\left(\frac{k_i}{{\delta }_i}\right).$$

$$\underset{N\to \infty }{\lim }\Vert \mathbf{Q}{\Vert }_{\infty }=\underset{N\to \infty }{\lim }\left[\underset{1\le i\le N}{\max }\left(\frac{k_i}{{\delta }_i}\right)\right]=c,$$

$$f(\delta ;\alpha ,\beta )=\frac{{\delta }^{\alpha -1}{e}^{-\frac{\delta }{\beta }}}{{\beta }^{\alpha }\Gamma (\alpha )}$$

Full text

Heterogeneity, whether in the nature of the components or in 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’ critical behavior  Pastor-Satorras and Vespignani (2001); Mieghem et al. (2009); Goltsev et al. (2012); de Arruda et al. (2017); Pastor-Satorras et al. (2015); de Arruda et al. (2018). As an illustration, consider two classical contagion models, the Susceptible-Infected-Susceptible (SIS) and the Susceptible-Infected-Recovered (SIR) models, which, when evolving on an homogeneous network, present a non-vanishing critical point Pastor-Satorras et al. (2015); de Arruda et al. (2018). However, the introduction of structural heterogeneity –in the form of broad degree distributions of the nodes– can result in a vanishing critical point Pastor-Satorras and Vespignani (2001); Chatterjee and Durrett (2009); Pastor-Satorras et al. (2015); de Arruda et al. (2018). More specifically, in the thermodynamic limit, a divergence of the second moment of the degree distribution Pastor-Satorras and Vespignani (2001); Boguñá et al. (2013); Pastor-Satorras et al. (2015); de Arruda et al. (2018) or a divergence in the maximum degree Pastor-Satorras and Vespignani (2001); Boguñá et al. (2013); Pastor-Satorras et al. (2015); de Arruda et al. (2018) imply a vanishing critical infectivity, which in turn has important practical implications for real-world networks as many display very broad  Chatterjee and Durrett (2009); Boguñá et al. (2013); Clauset et al. (2009) –or even, scale-free– degree distributions Voitalov et al. (2018). While structural heterogeneity is widely accounted for, heterogeneity in the dynamical parameters has received considerably less attention until recently. Indeed, it was mainly studied for the SIR model. A message passing formalism was proposed in Karrer and Newman (2010); Sherborne et al. (2018) and an heterogeneous mean-field approach in Gou and Jin (2017). In the latter, the authors also performed numerical experiments showing that the population can be more vulnerable in this scenario. More recently, this problem was examined on temporal networks in Darbon et al. (2018). Here, we focus on a different type of dynamical heterogeneity by providing the first characterization of the SIS critical point when recovery rates are distributed heterogeneously across the population. This case is empirically relevant because heterogeneous recovery rate distributions have been associated with biological differences between individuals Segal and Hill (2003); Fryer et al. (2010), demographic characteristics Dorjee et al. (2013) and social differences that result in non-homogeneous access to the health system Barber et al. (2017). In addition, we consider also the case in which correlations arise between structure and dynamics and we show, analytically and numerically, that such correlations can result in the extreme cases of power-law networks with non-vanishing critical points, and of homogeneous networks with a vanishing critical point. Our results complement previous evidence on the SIR model Gou and Jin (2017) and imply that a 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, and resource management.

SIS model with heterogeneous recoveries. We start by considering a population composed of $N$ individuals with an arbitrary pattern of connections, which can be represented as a network and is described by its adjacency matrix $\mathbf{A}$, which is usually assumed to be symmetric. Each individual can be in one of two states: (i) infected ($Y_{i}=1$) or (ii) susceptible ($X_{i}=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, for each directed edge, $i\sim j$, emanating from the infected individual $i$, we associate a Poisson process with rate $\lambda_{ij}$, $N^{\lambda_{ij}}(t)$ ($Y_{i}+X_{j}\rightarrow Y_{i}+Y_{j}$). Additionally, for each infected individual, we associate a Poisson process with rate $\delta_{i}$, $N^{\delta_{i}}(t)$, modeling the recovery ($Y_{i}\rightarrow X_{i}$). This system is statistically described using the order parameter, $\rho$, and the susceptibility, $\chi$, defined as

[TABLE]

where $n_{I}$ is the number of infected individuals. Both quantities can be directly estimated using Monte Carlo methods, in particular, the quasi-stationary method (QS) de Arruda et al. (2018) and the Gillespie algorithm de Arruda et al. (2018), where each of the aforementioned processes are simulated and the state of the nodes is evaluated de Arruda et al. (2018).

In the quenched mean field approach (QMF) one implicitly assumes that $\left\langle X_{i}Y_{j}\right\rangle\approx\left\langle X_{i}\right\rangle\left\langle Y_{j}\right\rangle$. Physically, this corresponds to neglecting dynamical correlations. Thus, denoting $y_{i}=\left\langle Y_{i}\right\rangle$, we have

[TABLE]

The standard SIS model considers that there is no variance in the dynamical parameters, i.e., $\lambda_{ij}=\lambda$ and $\delta_{i}=\delta$. As a consequence, it is possible to re-scale time and reduce the parameter space by defining $\tau=\frac{\lambda}{\delta}$. Eq. 2 is thus an upper bound Mieghem et al. (2009) of $y_{i}$ and, consequently, a lower bound for the critical point, which is calculated as $\tau^{\mathrm{QMF}}_{c}=\left(\Lambda_{max}(\mathbf{A})\right)^{-1}$. Here $\Lambda_{max}(\mathbf{A})$ is the leading eigenvalue of the adjacency matrix. Note that, for power-law networks (PL), in the thermodynamic limit, the critical point goes to zero if the maximum degree is a growing function of the network size. On the other hand, in the case of a contact process (CP), the spreading rate is defined as $\lambda_{ij}=\frac{\lambda}{k_{i}}$, and is thus described by the probability transition matrix $\mathbf{P}_{ij}=\frac{\mathbf{A}_{ij}}{k_{i}}$. In this case, the critical point is finite and $\tau^{\mathrm{QMF,CP}}_{c}=1$, regardless of the underlying structure.

Next, we focus on the case of heterogeneous recovery rates. To this end, consider the most general set-up with heterogeneous parameters allowing an arbitrary distribution of $\delta_{i}$ and $\lambda_{ij}$. Denoting $\mathbf{y}=\left[y_{i}\right]$, near the critical point, we can perform a linear stability analysis of Eq. 2, hence

[TABLE]

where $\Delta$ is a diagonal matrix, whose diagonal elements are $\Delta_{ii}=\delta_{i}$, $\lambda_{ij}=\lambda\mathbf{W}_{ij}$ is the rate matrix, and $\mathbf{A}\circ\mathbf{W}$ is the Hadamard product between the rate and the adjacency matrices. At the steady state, i.e., when $\dfrac{d\mathbf{y}}{dt}\rightarrow 0$, we have

[TABLE]

that will be positive iff the spreading rate is larger than the leading eigenvalue of $\mathbf{Q}=\Delta^{-1}\left(\mathbf{A}\circ\mathbf{W}\right)$, which in turn yields a critical point given as

[TABLE]

Observe that the elements of $\mathbf{Q}$ are the expected number of contacts before recovery. Obviously, the critical point simplifies to $\tau^{\mathrm{QMF}}_{c}=\left(\Lambda_{max}(\mathbf{A})\right)^{-1}$ in the homogeneous case, i.e., when $\delta_{i}=\delta$ and $\lambda_{ij}=\lambda$. The same applies to the CP. Note that, similarly to the homogeneous case, this prediction is an upper bound for the heterogeneous recovery rate scenario, because it relies on the independence of the random variables. In other words, if $i\sim j$, then $\mathbb{P}\left(Y_{i}=1|Y_{j}=1\right)\geq\mathbb{P}\left(Y_{i}=1\right)=y_{i}$, then the nodal probability is always overestimated (see  Mieghem et al. (2009) for a similar argument). From here onward, we fix $\lambda_{ij}=\lambda$ and focus on the effect of the recovery rate distribution on the critical point. Note that Eq. 5 can be bounded using a matrix norm. Thus, using the 2-norm, we obtain our first result,

[TABLE]

where $\|\mathbf{A}\|_{2}=\sqrt{\Lambda_{\max}\left(\mathbf{A}^{T}\mathbf{A}\right)}$ and, more specifically, $\|\mathbf{A}\|_{2}=\Lambda_{\max}\left(\mathbf{A}\right)$ for undirected networks. Eq. 6 therefore provides bounds on the leading eigenvalue of $\mathbf{Q}$ given by the structure and the variance of $\delta_{i}$.

Synthetic networks. To further characterize the critical behavior of our model, we first consider an Erdös – Rényi network (ER) with $N=10^{5}$ and $\langle k\rangle\approx 10$ (therefore $\tau_{c}^{\mathrm{QMF}}\approx 0.1$), which has a homogeneous structure and allows us to analyze the structural and dynamical effects independently. In Krylova and Earn (2013); Clancy (2014) the authors showed evidence in real data that the inter-infection time follows a gamma distribution. Consequently, the rate distribution must follow an inverse-gamma distribution. Therefore, we impose the recovery rates to have an inverse-gamma distribution, $\delta\sim\Gamma^{-1}(\alpha,\beta)$,where $\alpha$ and $\beta$ are the shape and scale parameters, respectively. Its mean is $\left\langle\delta_{i}\right\rangle=\frac{\beta}{\alpha-1}$ and its variance is $\text{Var}(\delta_{i})=\frac{\beta^{2}}{(\alpha-1)^{2}(\alpha-2)}$, for $\alpha>2$. In order to allow the comparison between different distributions, we restrict the distributions to unitary mean. In Fig. 1 we present the critical behavior of an ER network for different shapes, $\alpha$. As $\alpha$ decreases, the variance of $\delta$ and, consequently, its maximum, also increases. Consistently with Eq. 6, the critical point also moves toward zero. The insets in the top panel emphasize the behavior of the predicted critical point as a function of $\alpha$ and its comparison with the estimations from the Monte Carlo simulations. As expected, for sufficiently large values of $\alpha$ the dynamics behaves similarly to the standard SIS model with uniform $\delta$, where the predicted threshold coincides (see top inset of Fig. 1). The agreement between analytical and simulated critical points is very good, as can be seen in the top inset of Fig. 1. We remark that a similar experiment was carried out considering a gamma distribution, whose results are presented in Appendix A.

Real-world networks. To obtain further insights into real-world epidemics we also consider an inverse-gamma distribution for the recovery rates using a real network. Figure  2 presents the results of simulations in two real networks: (i) the UC Irvine messages social network Opsahl and Panzarasa (2009); Kunegis (2013) and (ii) the open flights network kon (2016); Kunegis (2013). These two networks represent different scales of a similar spreading process: the social network corresponds to a spatially localized network, while the open flights one captures a wider spatial scale. In the top inset of the bottom panel, we observe that the critical point predictions are remarkable for inverse-gamma recovery rates, even for these real networks.

From figures 1 and 2 we observe that the critical point decreases as we increase the variance of the recovery rate distribution. The critical point predictions for the standard process using QMF are a lower bound. However, when we consider the heterogeneous case, assuming an average recovery rate in the QMF is not enough to provide an adequate characterization of the process. In fact, it is not a lower bound anymore (see Fig. 1). The proper correction for the QMF predictions are given by Eq. 5, which is a lower bound for the underlying process.

Effects of dynamics-structure correlations. The bounds in Eq. 6 implicitly assume that there is no correlations between structure and dynamics. From the Gershgorin circle theorem we know that every eigenvalue of $\mathbf{Q}$ lies at least in one of the disks $D(Q_{ii},R_{i})$, centered in $Q_{ii}$ with a radius given as $R_{i}=\sum_{i\neq j}|\mathbf{Q}_{ij}|$. Therefore, considering a symmetric matrix, $|\Lambda_{k}|\leq\mathbf{Q}_{ii}+R_{k}$, hence $\Lambda_{\max}\leq\|\mathbf{Q}\|_{\infty}$, where the infinity norm is defined as

[TABLE]

If the structure and the dynamics are uncorrelated, Eq. 6 is a better bound. However, if they are correlated, Eq. 7 might give us further insights. For instance, for the PL case, the leading eigenvalue of $\mathbf{A}$ diverges in the thermodynamic limit leading to a vanishing critical point. Conversely, using Eq. 7 and a proper choice of $\delta_{i}$’s, we can change this behavior. In fact, assuming that $\delta_{i}(k_{i})\propto k_{i}$ in the thermodynamic, we have

[TABLE]

where $c<\infty$ is a finite real constant. This radically changes the critical behavior of the dynamics. Note that both the CP and the $\delta_{i}=k_{i}$ cases are described, at first order, by the probability transition matrix, $\mathbf{P}$, yielding to $\tau_{c}^{CP}=\lambda_{c}=1$.

In Fig. 3 we analyze the structure-dynamics correlation effects. In Fig. 3 (a) and (b) we perform a finite size analysis, comparing both processes on top of PL networks. In (a) we present the CP, where we can already observe that the critical point predictions are not as accurate as for the previous case, in alignment with the predictions reported in Ferreira et al. (2011); Mata et al. (2014). Thus, the mismatch between prediction and estimated critical point, in this case, seem to be related to dynamical correlations, which is neglected in the QMF. In Fig. 3 (b) we consider the recovery rate distribution as $\delta_{i}=k_{i}$, whose results suggest a finite critical point. However, the convergence seems to be slower if compared with the CP case.

The previous results show that it is possible, for the same structure, to have a vanishing critical point in the standard model and a non-null critical point when recovery rates are distributed. The opposite scenario is also possible. To show this, we consider an ER network with $\langle k\rangle\approx 10$ with $\delta_{i}=\frac{k_{i}}{k_{PL}}$, where $\mathbb{P}\left(k_{PL}\right)\sim k_{PL}^{-\gamma}$. That is, we now have a homogeneous structure and a heterogeneous recovery rate distribution. In Fig. 3 (c) we show a finite size analysis for this configuration varying $\gamma=2.1,2.7,3.5$. We observe that for $\gamma=2.7$ and $\gamma=3.5$ the underlying structure plays an important role, maintaining a non-vanishing critical point (see inset of Fig. 3 (c), where both curves have a slope close to zero). However, for $\gamma=2.1$ our results indicate the existence of a vanishing critical point (see Fig. 3 (c) inset). It seems reasonable to hypothesize that the scenario observed when $\gamma=2.1$ is due to the fact that, in the steady-state, the infection probabilities are inversely proportional to the nodal recovery rate and thus, that the evaluation of the recovery time at both ends of every edge enables an infection-reinfection mechanism. What are the necessary and sufficient conditions to observe this phenomenology needs, however, further exploration.

Furthermore, note that our model plays an important role in the localization of the leading eigenvector. As shown in Pastor-Satorras Romualdo and Castellano Claudio (2016) and recently formalized in Liu and Mieghem (2018), the eigenvector is localized in a sub-extensive portion, i.e., $IPR\sim O(N^{-\nu})$, where $0<\nu<1$ and $IPR=\sum_{i}^{N}v_{i}^{4}$ and $v$ is the normalized leading eigenvector. Thus, in the fully delocalized case, $IPR\sim O(N^{-1})$. As it can be seen in the top inset of Fig. 3 (c), when heterogeneous recovery rates are considered, the localization of the disease might also change. In one limiting case, the leading eigenvector of $\mathbf{P}$ is homogeneously distributed, therefore $IPR\sim O(N^{-1})$ and fully delocalized. However, as shown in the same inset, one can consider a structure that is delocalized for the standard case, but that becomes localized when a recovery rate of the form $\delta_{i}=\frac{k_{i}}{k_{PL}}$ (as introduced above) is considered. We remark that the control of the localization of diseases is still an open problem.

In summary, here we have analyzed the impact of heterogeneity in the recovery rates, allowing it to be arbitrarily distributed. We showed that dynamical heterogeneity is as important as the structural one, and that it can induce drastic changes in the SIS critical properties. Furthermore, an important consequence of our results is that the QMF approach provides a lower bound for the standard SIS, and hence gives a conservative prediction of the critical threshold. However, the standard formulation is not a lower bound for the heterogeneous case, i.e., when assuming $\delta=\left\langle\delta_{i}\right\rangle$ in the classical formulation. To solve this inconsistency, we proposed a solution that relates the structural and dynamical features by the spectral properties of the new matrix $\mathbf{Q}$. Thus, the new formulation opens the path for future research regarding resource allocation, as $\delta_{i}$ can be associated to the availability of resources to recover individuals. Aside from 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, our results might also relate to the predictability of complex systems, and in particular, of diseases Scarpino and Petri (2017). At the same time, our findings raise intriguing questions about the consequences of potential heterogeneities in spreading rates, $\lambda_{ij}$ (as suggested by Eq. 5), and also in other dynamical processes. For example, in the case of Kuramoto oscillators, correlations between the natural frequencies and node degrees change the nature of the transition. Uncorrelated natural frequencies present a second-order phase transition, while correlations might introduce a first-order phase transition on PL networks Gómez-Gardeñes et al. (2011). This phenomenology contrasts with what is observed for the SIS model, where, while the phase transition is still second-order, the usual vanishing critical point changes to a well-defined transition.

Acknowledgement

Research carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (grant 2013/07375-0). Y. M. acknowledges partial support from the Government of Aragón, Spain through grant E36-17R, and by MINECO and FEDER funds (grant FIS2017-87519-P).

Appendix A Gamma distributed recovery rates

To further characterize the critical behavior of our model, we first consider an Erdös – Rényi network (ER) with $N=10^{5}$ and $\langle k\rangle\approx 10$ (therefore $\tau_{c}^{\mathrm{QMF}}\approx 0.1$), which has an homogeneous structure and allows us to analyze the structural and dynamical effects independently. We impose the recovery rates to have a Gamma distribution, $\delta\sim\Gamma(\alpha,\beta)$, whose p.d.f is expressed as

[TABLE]

where $\alpha$ and $\beta$ are the shape and scale parameters respectively and $\Gamma(\alpha)$ is the gamma function evaluated at $\alpha$. Moreover, its mean is $\left\langle\delta_{i}\right\rangle=\alpha\beta$ and its variance is $\text{Var}(\delta_{i})=\alpha\beta^{2}$. In order to allow the comparison between different distributions, we restrict the distributions to unitary mean by setting $\beta=\alpha^{-1}$ (hence $\text{Var}(\delta_{i})=\alpha^{-1}$). In Fig. 4 we present the critical behavior of an ER network for different shapes, $\alpha$. As $\alpha$ decreases, the variance of $\delta$ and, consequently, its maximum, also increases. Consistently with Eq. 6, the critical point also moves toward zero. The insets in the top panel emphasize the behavior of the predicted critical point as a function of $\alpha$ and its comparison with the estimations from the Monte Carlo simulations. As expected, for sufficiently large values of $\alpha$ the dynamics behave similarly to the standard SIS model with uniform $\delta$, where the predicted threshold coincides (see top inset of 4). Although the agreement between analytical and simulated critical points decreases for very heterogeneous rate distributions, the analytical values for the critical points are always below the simulated ones and thus provide a –safe– lower bound on the critical threshold.