Effective medium estimate of spreading growth in random networks with Lévy flights
The spreading of epidemic in random spatial metapopulation networks with power-law decaying connectivity is investigated. We perform numerical simulations of reaction-diffusion processes with two compartments on the simplest one-dimensional ring topology of metapopulation graphs and use the effective medium approximation to estimate the spreading velocity of the infection, which we find to grow exponentially in time. We also provide an upper bound for the diameter of the infected populations. Our study outlines a direct and simple approach to estimate possible outbreaks of epidemic spreading in very general mobility patterns.
Understanding the spread of emergent infectious diseases in the geographic space is of fundamental importance in an increasingly connected world. In ancient times, the spreading of epidemics, such as the black death, could be understood in terms of diffusive processes Barabási (2016). In those cases the disease is spread by the agents/hosts that can only travel with low velocities bounded by the local connectivity of the geographical lattice. This gives rise to a wave-front of infected individuals, which travels at a finite speed. The recent great increase of the connectivity among densely populated areas and the correspondent urbanization, has increased the risk that infectious diseases will spread. The complexity of human mobility at all scales, being that urban and inter-urban or world-wide, is reflected in the possibility for the infection to cross arbitrary distances in close to no time. As a consequence, the number of infected sites grows exponentially fast, as opposed to linearly. Similar phenomena are also discussed in the biological context, see Hallatschek and Fisher (2014) and references therein.
The standard approach to model infection spreading on a global scale are metapopulation models. A metapopulation is a geographically localized subpopulation that is assumed to be well mixed. In a fashion similar to reaction-diffusion processes, a local quantity like the density of infection carriers, may diffuse to neighboring metapopulations and in addition interact locally according to certain laws Colizza et al. (2007); Colizza and Vespignani (2008). The metapopulation approach has been successfully used to describe spatially embedded populations, such as cities and urban areas, interacting with each other Van den Broeck et al. (2011). Here, we assume diffusive coupling between the populations as for the case where the individuals travel in the metapopulation and forget about their original population at each time step. Thus the model is Markovian on the metapopulation level. Instead, they perform a random walk over the populations with jumping probabilities that are given according to a travel rate matrix.
This gives rise to a network description of the connections between the populations. Networks have proven to be a fundamental tool to understand complex systems and are very well suited to study the problem of epidemic spreading and human mobility in general Pastor-Satorras et al. (2015); Piontti et al. (2014); Hufnagel et al. (2004); Brockmann et al. (2006); Gonzalez et al. (2008). Paired with numerical simulations, complex network theory provides an accurate representation for the transmission dynamics in the increasingly complex world Halloran et al. (2017); Salathe et al. (2012).
In this work, human motion is modeled by a random network. Random networks, although being complicated objects, can often be characterized by a small amount of parameters due to the universal behavior of a whole ensemble. Furthermore they are versatile enough to characterize the most general transportation system. Often, the analytical treatment of an ensemble of random networks is more viable than of a single representation. We therefore demonstrate that knowledge of the universality class of a transport network is enough to extract crucial information about the infection process. In particular we use the effective medium approximation (EMA) to replace a random network with a typical deterministic representative network, that is analytically tractable.
Effective medium theory is an ancient tool, originally developed to determine material properties in continuous and lattice model of heterogeneous media Bruggeman (1935); Choy (1999); Kirkpatrick (1973). EMA is routinely applied in systems with short-range connections, but only recently was generalized to systems with long-range connections Thiel and Sokolov (2016). The latter theory is employed here to a disordered reaction-diffusion system with long-range connections decaying as power-laws with the geographical distance. This setup serves as a prototype for infection spreading in the human travel network.
Our theoretic reasoning will be laid out in section II. We first give a brief introduction into the SIS metapopulation model, and then show how the spread of an infectious disease can be estimated in general. We proceed with introducing random transition rates and using EMA to estimate the infection spread in the random case. The numerical simulations of section III support our arguments nicely. We close with discussion in section IV. A small calculation can be found in Appendix A.
ii.1 SIS metapopulation model
We will demonstrate our idea in a very simple susceptible-infected-susceptible (SIS) metapopulation model on a line. The SIS model is a two compartment scheme where individuals belongs to only two categories: susceptible individuals that are not affected by the disease and the infected individuals, the one who carry the disease and are able to spread it to the susceptible. The susceptible individuals turn into infected according to the reaction scheme , with infection rate . Furthermore, the infected will spontaneously recover with rate , so that we have the additional reaction . Contrary to the more advanced models with recovery compartments, here the infected individuals will keep participating in the dynamics and eventually turn into susceptible through the cycle . The infection can be sustained as long as the basic reproductive number , defined as the average number of secondary infections caused by a primary case introduced in a fully susceptible population, is larger than unity Anderson et al. (1992). Otherwise the epidemics will eventually die out in the thermodynamic limit of a infinite population for . When a single infective agent will lead to the steady state with a finite average density of infected, which corresponds to the endemic state of the disease, where the fraction of infected population is given by the stationary density of infected individuals Barrat et al. (2008).
In the metapopulation model all populations have the same amount of individuals and are placed at equal distances. Each of them is enumerated by an integer , so that for an infinite system the set of all labels is defined by the set of integers . We allow individuals from to travel not only to adjacent metapopulations , but also to distant metapopulations . The number of people traveling from to per unit time is the product of the rate with the population size . We assume the rates to be symmetric , so that the flux of travelers in one direction is balanced by the flux in the opposite direction. The rates can be interpreted as the weights on links of a weighted graph with node set . We stress that the node index denotes a position in Euclidean space. Therefore a link that stretches only a unit distance in the topological sense, may connect two very remote nodes in Euclidean space. These long-range connections are our model for the fast superdiffusive transport in the international human travel network.
We denote with the transport operator Rvachev and Longini (1985), i.e. the negative Laplacian or Kirchhoff matrix of this graph, defined as
where is the Kronecker symbol.
Let denote the fraction of the population on site at time . The time evolution of is described by the Master equation:
As it will be used later, we associate the initial condition with Eq. (2).
The total number of individuals remains constant. We normalize this constant to unity and denote with the fraction of the infected population at and . The evolution equation for reads
The transport operator accounts for the diffusion of individuals between adjacent populations, while
is the SIS reaction term describing recovery and infection. For the simpler susceptible-infected (SI) model one simply puts . We already used the normalization and identified the susceptible fraction of the population at site with . See Fig. 1 for a sketch of the interplay between reaction and diffusion in the metapopulation model.
Initially we infected a fraction of node ’s population, i.e. .
ii.2 Estimation of infection spreading
Our arguments in this section go parallel to those of Ref. Mancinelli et al. (2003). Using the Feynman-Kac-equation, we can formally write the solution of to this initial value problem as an expectation value Sato (1999)
Here, is a random walk on , that starts at and evolves according to the master equation Eq. (2), and the average is taken over all such random-walk configurations. Hence the transition probabilities of are given by
Although Eq. (5) is a complicated self-consistent equation that can hardly be solved for , it provides a very useful estimate for the latter. Since , we have the obvious estimate . Plugging this inequality into (5) and identifying the transition probability , we obtain:
This remarkable relation shows, that is related to the free diffusion of agents on the graph with rates .
Eq. (7) is useful, if we can actually find an analytical solution for the free motion of agents. Then, the inequality can be inverted for to estimate the spreading of the infection front. To do so, one considers the level set, i.e. the set of nodes , where the fraction of infected is higher than some threshold value , i.e. .
If one assumes that the free diffusion is Gaussian, i.e. with diffusion constant , plugging this in Eq. (7) one ends up with:
which yields after partially solving for :
One can then relate the diffusion constant to the rates appearing in the Laplacian to obtain a velocity that increases monotonically with the rates Belik et al. (2011). This shows that if are the solutions for the standard diffusion problem then there is an upper bound for the front propagating at constant speed.
The crucial assumption for Eq. (10) is the Gaussianity of (see also the discussion in Mancinelli et al. (2003)). This is justifiable in the normal diffusive case, when the infected individuals’ displacements have a finite scale and are in that sense “small”. Their trajectory, being a sum of many of those small displacements, is Gaussian by virtue of the central limit theorem. However, one of the most prominent features of the human travel network is a violation of this principle. Nowadays, disease carriers can travel long distances in short times with ease. One should no longer assume that the individuals’ displacements have a finite scale. The Gaussian form of needs to be replaced with a stable (also called Lévy-) distribution, which plays the same role in the generalized central limit theorem Feller (1971). An important characteristic of these distributions is a power-law decay for large :
where is a sort of generalized diffusion constant, and characterizes the class of stable distribution. Plugging this assumption into Eq. (7), one has for large :
Every node obeying this inequality has a higher infected fraction than . Assuming the infection started on the origin, the diameter of this set is given by twice of above bound, due to symmetry. Hence the diameter of the infected region can be estimated as
We see that the infection region grows exponentially fast, which is a consequence of our power-law assumption on .
ii.3 Random networks
Consider now the case, that the rates are random variables. We assume that the rates of different links are statistically independent of each other. Hence , but both are independent from . We will assume that the rates decay over long distances like a power-law. In particular, we assume that
where are i.i.d. random variables chosen for each link. In practice we will sample the from a uniform distribution, although in principle any distribution with a finite mean is viable. For example a distribution for that mimics a more realistic human mobility patterns would take into account the scarceness of the long-range connections (flights), and the abundance of the short-range connections (cars, bus, metro and bicycle).
ii.4 Effective medium approximation
Using the effective medium approximation (EMA), we can replace the ensemble of the random graphs with a deterministic graph, that is analytically tractable, that exhibits the same qualitative transport behavior, and that enables us to make quantitative predictions on the propagation speed. In a certain sense, EMA enables us to find the “average” Kirchhoff matrix to the random in such a way that the average transport fluxes remain unchanged. The graph corresponding to is called the “effective medium”; it has the deterministic transition rates . One is free to choose the effective medium as long as the two necessary conditions are fulfilled: (i) one never erases a link that could be present in the random network, and (ii) one respects the spatial scaling of the random transition rates Thiel and Sokolov (2016). In our case this means, that the effective medium is fully connected and that decays with the same power law as in Eq.(13). The effective medium is determined by the following condition: If one link in the effective medium were to be replaced with its random original link, then the transport flux along this link should not change on the average. This prescription leads to a self-consistent equation for that has an easy solution in our special case Thiel and Sokolov (2016):
denotes the disorder average over the ensemble of the random networks, i.e. over the distribution of the in our case. Eq. (2) has the following form in the effective medium (here we use the special form of ’s diagonal entries):
We show in Appendix A that defined like this indeed decays like a power-law:
This defines the constant in Eq. (12) and provides an effective medium estimate for the infection diameter. The later is defined as the number of infected populations at time and is given by
The last equation is one of our main results. Let us restate our reasoning. We are given an infection process in a random environment represented by a symmetric weighted graph with random transition rates. We want to estimate the prefactor in the time dependence of , which can be interpreted as a generalized velocity of the infection spreading. The latter is obtained from the diameter of the infected region as
The upper bound for the diameter is determined from the effective medium approximation of the random network by means of the Feynman-Kac equation and it also gives the upper bound for the condition
Thus the upper bound for the generalized infection velocity in the fully connected topology with long-range Lévy jumps is independent on the infection and recovery rates and as well as on the threshold concentration. Consequently, it is independent from the infection strength on the local scale, which is quantified by the reproductive number . It only depends on the initial concentration of the infected seed population and on the EMA parameters.
This implies that the topology plays the key role in identifying the threshold for the global outbreak. The model employed here mimics the long-range connections found in the human travel network. In the following section we will demonstrate our prediction Eq. (17) in numerical simulations.
Iii Numerical results
We start by solving numerically Eq. (3) without a recovery term, i.e. with on a ring with sites, which correspond to the simplest SI model with a saturation to all individuals infected.
The transition rates have been determined for each link by Eq. (13) with and drawn from a uniform distribution over . Therefore in our simulations. It is important to note that the metric of a ring in one dimension is used
This determines the upper triangle of ; the lower triangle is given by the symmetry condition of . Its diagonal elements are the negative sum of all other elements in the respective column.
Using and a Runge-Kutta method of order 5, we solve the equation with a localized initial condition , where , for each realization of and obtain a collection of . Then we compute as the average over the single environment realizations, i.e. the average over the disorder. From this quantity, given the infection threshold value , we determine the infected cluster by finding the set of populations , and then we compute the diameter (in circular distance) of this set.
A comparison of this quantity with Eq. (17) is given in Fig. 2(a), where we also show the generalized velocity estimated from Eq. 18. However the time series for is truncated very early, because it saturates soon after. (Then the infection has basically reached the other side of the ring.) is always below the theory upper bound given by Eq.( 19), confirming our prediction. Also the long-range connections make finite size effects in our simulations very severe. Since the diameter obtained from simulation follows the analytic prediction Eq. (17) pretty well, we can extract some of the parameters from the exponential fit
Comparison with Eq. (17) gives:
In principle the EMA parameters and and the net infection rate could be extracted this way. This may however be a hard task, because the non-linear term is not easy to detect in the exponential fit. The combination on the other hand is easy to obtain and it can also be measured from the slope of the tail of . In our simulations of the SI metapopulation model we obtain , which gives . This is a reasonably close value to the Lévy exponent used for generating the graph realizations in the first place.
As a next step, we repeated the whole procedure in a SIS model with , and , which gives a basic reproductive number of new infection cases per time step in a fully susceptible population. The infection prevalence curves obtained with the SIS metapopulation model are shown in Fig. 3, with the reference concentration threshold . The time gap between the outbreaks of the first and last population infected is about time steps, and the absolute global infection time is time steps. As expected the results are similar to the previous case, see Fig. 2(b). However in this case the estimation of the Lévy exponent is worse than in SI case. Here it is given by the fit parameter , which results in , a value more than twice the real value . This overestimation is in reality a purely finite size effect of our numerical simulations due to the extreme long-range connections, as we show in Fig. 4. Increasing the number of populations (nodes), the difference between growth rate evaluated from the numerical fit and the theoretical value decays faster than a power-law. In particular the total error decrease as a power-law with the system size . For the maximum probed value of , for , we find , which is a good estimation of the theoretical value.
To investigate sample to sample fluctuations, we measured in each environment. That means we obtained the infection diameter from each realization and fitted a line through the tail. The slopes of these lines and the corresponding value of obtained from the infection rate are presented respectively in the box-plots in Fig. 5.
The comparison of the exponential growth coefficient between theory and fit value as a function of is given in Fig. 6. For small value of the numerical fit estimation is very poor, while after the value the theoretical and numerical value eventually merge.
The overall agreement with the theoretical prediction of Eq. 17 for the probed range in is shown in Fig. 7. The inequality becomes more and more close to a proper equality by increasing . For the higher dimensional case we can expect that this behavior is qualitatively the same, since the growth is still governed by the exponential parameter given by .
The goal of this paper was to present a new analytical tool for reaction-diffusion problems in random media. We wanted to advocate the use of effective medium theory that provides a deterministic representative for an originally random network. Together with the standard Feynman-Kac argument we provide an upper bound for the infection spread in a simple SIS model. This way we demonstrated that EMA is still relevant even beyond the short-range connection paradigm.
With the human travel network in mind, we presented a meta-population model with random long-range connections. Due to the presence of long-range connections we find the exponential growth of the infection diameter. The growth rate depends on both the infection and recovery rates and as well as on the topology encoded in the Lévy exponent of the statistical decay of the link strength, see Eq. (13). Other characteristics of the transition rates (like their mean) only determine the prefactor. Notice, that long-range links with a “weak” power law – i.e. in Eq. (13) – would eventually lead to a ballistic growth of the infection front. These results are also discussed in Ref. Hallatschek and Fisher (2014).
EMA is known to nicely reproduce the transport behavior of a random system provided it is far away from the percolation threshold. The networks we treated here are very well connected due to the presence of the long-range links. Therefore they are generically far from percolation threshold which partly explains the success of our approach.
It might seem that the model considered here lacks realism as we are considering a one dimensional array of populations. However the long-range connections dominating the dynamics make our results quite general and applicable to other underlying topologies. In particular a measurement of , available from a simple linear fit of the logarithm of the infection diameter, can give a rough estimate of the exponent of the transition rates. The numerical fit estimation of the EMA parameters becomes reliable for the SIS model only for sufficiently large networks as finite size effects are quite severe in this case.
We believe that modifications of long-range EMA will remove some restrictions that we imposed in this article. For example, modifications of EMA are necessary to deal with random networks where the transition rates are not symmetric so that variable population sizes can be taken into account. Furthermore, when it is possible to deal with correlated links, more realistic models than a simple grid of the metapopulation’s locations can be included.
Future work and applications of EMA to reaction-diffusion systems include the generalization to arbitrary heterogeneous connectivity networks, such as real-world networks of human mobility.
Acknowledgements.The authors are indebted to A. Vulpiani for insightful discussions. This work has partially been founded by the DFG / FAPESP, within the scope of the IRTG 1740 / TRP 2015/50122-0.
Appendix A Asymptotic behavior of
The equation is solved using Fourier transform, i.e. we multiply on both sides and sum overall . Defining , we obtain:
Here is the polylogarithm function. The expression on the right hand side before defines the Fourier symbol of the transport operator. Using the polylogarithm’s expansion around (obtained by Mathematica) one finds the following small wave-vector expression for :
Note that the sign of is negative for all . As we now have , the solution is given by
where we used here the initial condition , which gives and identified
The expansion also shows that is asymptotically equal to a symmetric stable distribution whose Fourier transform is exactly given by our stretched exponential. The PDF of such a random variable decays like a power law for large , Feller (1971):
Here denotes the inverse Fourier transform. Using this equation and the definition of , with , and , we recover Eq. (16) from the main text.
- Barabási (2016) A.-L. Barabási, Network science (Cambridge university press, 2016).
- Hallatschek and Fisher (2014) O. Hallatschek and D. S. Fisher, Proceedings of the National Academy of Sciences 111, E4911 (2014).
- Colizza et al. (2007) V. Colizza, R. Pastor-Satorras, and A. Vespignani, Nature Physics 3, 276 (2007).
- Colizza and Vespignani (2008) V. Colizza and A. Vespignani, Journal of theoretical biology 251, 450 (2008).
- Van den Broeck et al. (2011) W. Van den Broeck, C. Gioannini, B. Gonçalves, M. Quaggiotto, V. Colizza, and A. Vespignani, BMC infectious diseases 11, 37 (2011).
- Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Reviews of modern physics 87, 925 (2015).
- Piontti et al. (2014) A. P. Y. Piontti, M. F. D. C. Gomes, N. Samay, N. Perra, and A. Vespignani, Network Science 2, 132 (2014).
- Hufnagel et al. (2004) L. Hufnagel, D. Brockmann, and T. Geisel, Proceedings of the National Academy of Sciences of the United States of America 101, 15124 (2004).
- Brockmann et al. (2006) D. Brockmann, L. Hufnagel, and T. Geisel, Nature 439, 462 (2006).
- Gonzalez et al. (2008) M. C. Gonzalez, C. A. Hidalgo, and A.-L. Barabasi, Nature 453, 779 (2008).
- Halloran et al. (2017) M. E. Halloran, K. Auranen, S. Baird, N. E. Basta, S. Bellan, R. Brookmeyer, B. Cooper, V. DeGruttola, J. Hughes, J. Lessler, E. T. Lofgren, I. Longini, J.-P. Onnela, B. Ozler, G. Seage, T. A. Smith, A. Vespignani, E. Vynnycky, and M. Lipsitch, bioRxiv (2017), 10.1101/198051.
- Salathe et al. (2012) M. Salathe, L. Bengtsson, T. J. Bodnar, D. D. Brewer, J. S. Brownstein, C. Buckee, E. M. Campbell, C. Cattuto, S. Khandelwal, P. L. Mabry, et al., PLoS computational biology 8, e1002616 (2012).
- Bruggeman (1935) D. A. G. Bruggeman, Annalen der Physik 24, 636 (1935).
- Choy (1999) T. C. Choy, Effective Medium Theory – Principles and Applications, International Series of Monographs on Physics (Oxford University Press, New York, 1999).
- Kirkpatrick (1973) S. Kirkpatrick, Reviews of Modern Physics 45, 574 (1973).
- Thiel and Sokolov (2016) F. Thiel and I. M. Sokolov, Physical Review E 94, 012135 (2016).
- Anderson et al. (1992) R. M. Anderson, R. M. May, and B. Anderson, Infectious diseases of humans: dynamics and control, Vol. 28 (Wiley Online Library, 1992).
- Barrat et al. (2008) A. Barrat, M. Barthelemy, and A. Vespignani, Dynamical processes on complex networks (Cambridge university press, 2008).
- Rvachev and Longini (1985) L. A. Rvachev and I. M. Longini, Mathematical biosciences 75, 3 (1985).
- Mancinelli et al. (2003) R. Mancinelli, D. Vergni, and A. Vulpiani, Physica D Nonlinear Phenomena 185, 175 (2003).
- Sato (1999) K.-I. Sato, Lévy Processes and Infinitely Divisible Distributions (Cambridge University Press, Cambridge, 1999).
- Fisher (1937) R. A. Fisher, Annals of Human Genetics 7, 355 (1937).
- Tikhomirov (1991) V. Tikhomirov, in Selected works of AN Kolmogorov (Springer, 1991) pp. 242–270.
- Belik et al. (2011) V. Belik, T. Geisel, and D. Brockmann, Physical Review X 1, 011001 (2011).
- Feller (1971) W. Feller, An Introduction to Probability Theory and Its Applications. Volume II, second edition ed. (John Wiley & Sons, Inc., 1971).