key: cord-0043384-zykknzel
authors: Iannelli, Flavio; Sokolov, Igor M.; Thiel, Felix
title: Reaction-diffusion on random spatial networks with scale-free jumping rates via effective medium theory
date: 2018-09-28
journal: nan
DOI: 10.1103/physreve.98.032313
sha: 3963a5bd0c5217113cd6533eb96bf7e3da5047d0
doc_id: 43384
cord_uid: zykknzel

We study epidemic processes using a metapopulation approach on the line featuring random transport rates between arbitrarily distant sites. An average transport network is found using a recently developed variant of the effective medium approximation (EMA) that is capable of dealing with these long-range connections. Using a Feynman-Kac argument in the effective medium, we derive an estimate on the size of the infected domain, and reproduce the known result of its exponential growth in time. We hereby demonstrate the applicability of long-range EMA to dynamical processes on networks more intricate than simple diffusion.

Network science has emerged in recent years as a fundamental theoretical framework for the modeling and understanding of large complex systems with many interdependent subunits [1] . Virtually any relation or interaction among any set of agents can be represented as a graph. Previous works have considered either the structural properties of the network itself [2, 3] or certain dynamics placed on the network, e.g., coupled oscillators [4, 5] or diffusive transport [6] [7] [8] [9] [10] [11] . The latter is the focus of this paper, in particular we will consider the transport of infectious pathogens.

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 [1] . In those cases the disease is spread by the agents/hosts that can only travel with bounded velocities between neighboring locations. 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 a different biological context, see Ref. [12] and references therein. Network theory's success stems not only from its versatility, but even more from the fact that many dynamical processes can be characterized and understood from the * iannelli.flavio@gmail.com † thiel@posteo.de underlying connectivity properties of the graph [13] [14] [15] [16] [17] [18] . Often, exact knowledge of the network connectivity eludes the theorist, because they are either quickly changing-as is the case in temporal networks [19] -or are hard to assess. In this case a practical approach is to model our ignorance with a random network.

Many techniques have been developed to deal with dynamical processes on random networks, among them the heterogeneous mean field [20] and the annealed adjacency matrix approximation [21] . The main rationale of statistical physics applies: many dynamical details of random networks are determined by a few parameters of the whole ensemble [2] .

To properly understand transport processes it is important to embed the network into the geographic space, i.e., one has to consider spatial networks [22, 23] . The simplest spatial networks are of course lattices. In the context of percolation theory, the so-called effective medium approximation has been developed to describe the diffusion on disordered-i.e., random-lattices [24, 25] . The idea is to replace the initially random transport rates between the nodes with fixed deterministic ones. This deterministic average network is called the effective medium, it is characterized by an effective diffusion coefficient. EMA is not a blind average of the transport rates, rather it is determined in a self-consistent manner. Would a link in the effective medium be replaced by its random original, the transport flux along this particular link would not change on average. Hence EMA is particularly suited for systems with independent links. In this paper, we employ EMA for infection spreading in the global human traffic network.

The spatial embedding of the network is especially important in global human traffic, as two topologically adjacent nodes (e.g., airports) may be geographically very far apart. Crucially, empirical observations show that human mobility lacks a definite scale [10, 26] and features long-range connections which have been a major limitation for EMA. Recently [27] , we have developed an EMA variant that overcomes this restriction and provides an analytical technique to deal with random spatial networks-a model of tremendous complexity. The goal of this paper is to demonstrate that and how EMA can be used in reaction-diffusion systems on random spatial networks with long-range connections. Contemporary fields where our proposed theory may become relevant are epidemic spreading in the global mobility network [28] [29] [30] [31] [32] or dispersal phenomena in biological contexts [12] .

For the remaining of the paper we are concerned with epidemic processes in a metapopulation, where the subpopulations are placed on an (ideally) infinite line with random long-range transport rates between them. The metapopulation approach has been successfully used to describe spatially embedded subpopulations, such as cities and urban areas, interacting with each other [33] . Here, we assume diffusive coupling between the subpopulations. The individuals travel in the metapopulation and forget about their original subpopulation at each time step. Thus the model is Markovian on the metapopulation level. Each individual performs a random walk over the subpopulations with jumping probabilities that are given according to a travel rate matrix.

This gives rise to a network description of the connections between the subpopulations. We leverage the Feynman-Kac formula to derive a bound for the diameter of the infected region in a deterministic model. Then we proceed to show that this estimate is as well realized in random models, where it can be computed from the effective medium approximation.

The rest of the paper is organized as follows: We start with Sec. II where we review how to obtain the infection diameter in a deterministic short-range metapopulation model. We proceed to review the necessary amendments for deterministic long-range models in Sec. III and finally consider random long-range models in Sec. IV. This is also where we explain the effective medium approximation. Numerical confirmation of our theory is presented in Sec. V. Discussion and concluding remarks are found in Sec. VI.

As an introductory example, we consider the susceptibleinfected-susceptible (SIS) metapopulation model on a line. Each lattice site in the metapopulation is occupied by a subpopulation, each with the same population size. The density of infected individuals at site x and time t is denoted by ρ x (t ), where x ∈ Z and t 0. The column vector of all ρ x is denoted without the subscript, i.e., ρ = (. . . , ρ x−1 , ρ x , ρ x+1 , . . .) T . Individuals travel to an adjacent subpopulation with a constant transport rate W , i.e., Wρ x infected will travel from x to x + 1 and to x − 1, per unit time. We assume that these travel rates are symmetric, so that the total number of individuals in each subpopulation does not change in time. The change in density of infected individuals can be expressed in terms of the following continuous-time reaction-diffusion equations:

The matrix is the transport operator that describes jumps between adjacent lattice sites [28] . In the current example it is equal to the graph Laplacian of the line:

x,y := W δ x,y−1 + W δ x,y+1 − 2W δ x,y .

(2)

Initially we infect a fraction c 0 of the site on the origin, i.e., ρ x (t = 0) = c 0 δ x,0 . Locally in each subpopulation, the infection dynamics take place. It is described by the reaction term

The second term describes the recovery of infected individuals with rate μ, I μ → S, when the infected become susceptible again. As the total number of individuals per subpopulation is conserved, the density of susceptible individuals is given by 1 − ρ x (t ). This reveals the first term in Eq. (3) as the infection of a susceptible individual with rate β, S + I β → 2I . Note that the local reaction rate is bounded from above by

The ratio R 0 = β/μ is the basic reproductive number and denotes the average number of secondary infections caused by a primary case in a fully susceptible population. The infection can be sustained locally in the long-time limit only when R 0 > 1, or equivalently when β > μ, which we will assume throughout the text. In this case a single infected agent in an otherwise susceptible population will lead to a steady state with a non-zero infection density given by (β − μ)/β [34] . A sketch is given in Fig. 1 where we show the interplay between reaction and diffusion in the metapopulation model. One way to investigate the spreading process is to apply the Feynman-Kac representation to solve Eq. (1). This approach was used, e.g., in Ref. [35] . Here, one considers a random walk X(t ), with X(t = 0) = 0 whose pdf P x (t ) is determined by Eq. (1) without the reaction term:

Then the solution of Eq. (1) is given by the following implicit equation:

where the average is taken with respect to the random-walk realizations {X(t )}. The Feynman-Kac equation relates the characteristic function of certain integrals of random processes with a partial differential equation and vice versa. The rationale is the same like in the construction of path integrals in quantum mechanics. Although it is hard to solve Eq. (6) explicitly, it is easy to find an upper bound. This is done by replacing the exponent via Eq. (4), by plugging in the initial condition ρ x (0) = c 0 δ x,0 and by recognizing δ X(t ),x = P x (t ):

For the second part of the expression we used the fact that P (t ) approaches a Gaussian for large times. An alternative way to derive the inequality is to linearize Eq. (1). The inequality Eq. (7) is useful, when one considers thē c-level set of the infected fraction, i.e., all sites x, such that 

which can be solved for x and yields a (time-dependent) radius of the infected region:

The diameter of the infected region is twice of the above radius and is asymptotically bounded:

This shows that for large times the infection spreads no faster than ballistically [36, 37] , with a velocity that grows monotonically with the transport rate W , which is indeed the case [38] . For a diffusion-limited infection this means that there is an upper bound for the front propagation speed. This is a consequence of the Gaussianity of P x (t ) (see the discussion in [35] ) which in turn is related to the lack of long-range connections.

With fixed reaction dynamics, the above reasoning can be extended in two directions: (i) the introduction of transport beyond the nearest-neighbour population, and/or (ii) make the transport rates a random quantity. This will be done in the next two sections.

To model the fast multi-scale human mobility, one might consider the introduction of more than nearest neighbor connections in the transport operator. Instead of Eq. (2) one might consider

where the transition rates W x,y are symmetric, i.e., W x,y = W y,x , and decay with the distance, so that the sum in the diagonal terms of converges. Consider first the example, when the transport rates decay like a power law with distance:

with α ∈ (0, 2) and where K plays the role of an anomalous diffusion constant. All reasoning from Sec. II can be repeated up to Eq. (7). However, the random walk generated by the transport operator of Eq. (11) with rates given by Eq. (12), is very different from before. Due to the possibility of long-range jumps that lack a finite variance, it will not converge to a Brownian motion, but instead to an α-stable distribution which is characterized by power-law tails instead of a Gaussian decay [39] :

These scale-free random walks are known as Lévy flights, and α ∈ (0, 2) is the Lévy exponent. A derivation of the previous equation is reproduced in the Appendix. Using this power law in Eq. (7) and solving for |x| allows us to estimate the diameter of the infected region:

which, contrary to the ballistic growth Eq. (10) found for bounded jumps, grows exponentially fast. One might argue that the assumed power law decay in the transition rates is rather specific and far off the measured travel rates. To overcome the this problem, we model our ignorance about the actual travel rates with chance.

We now consider the reaction diffusion Eq. (1) with a transport operator Eq. (11) that features random rates W x,y . It is assumed that the rates are random variables independently placed on each link (x, y). They are symmetric W x,y = W y,x and decay with the distance between the nodes |x − y| such that the diagonal terms in Eq. (11) are well defined. Hence, there is a family {p (x,y) (w)} of probability density functions, that describe the distribution of W x,y and that in total describes the ensemble of random networks.

The Feynman-Kac equation could still be used, but it would involve the random walk in a random network generated by the random operator . Since this is a rather hopeless venture, we will first employ EMA to compute an average transport operator . For details we refer to Ref. [27] . We call the average network described by the effective medium. The deterministic rates W x,y of the EMA operator have to be chosen such, that (i) any link that is present in some network of the ensemble will be present in the effective medium, albeit with possibly different strength; and such that (ii) the distance scaling in p (x,y) (w) is preserved. These conditions are necessary for the effective medium to be well defined. The transport rates are determined by the following set of self-consistency equations:

Here, the average is taken over the distribution of one fixed transport rate W x,y . R x,y is the so-called resistance distance [40] computed from the (pseudo-)inverse of :

The expression in Eq. (15) describes the average change in the stationary transport flux upon replacement of the effective medium link along (x, y) with its random original W x,y . EMA requires this change to vanish on average. For this reason it is very successful in reproducing the diffusive properties of the random network ensemble. It is important to note that the actual choice of the effective medium graph is mostly arbitrary, as long as the two conditions given above are respected. Equation (15) constitutes a set of equations for each class of links that share the same distribution. It simplifies considerably, if one assumes scaling behavior between distance and rates. We will focus here on the simplest case, when the transport rates are given by some independent and identically distributed (i.i.d.) random number divided by a power of the distance

where α ∈ (0, 2) and Z x,y is a family of i.i.d. random variables. In Ref. [27] it was shown that the actual distribution of the rescaled transition rates Z x,y does not influence the qualitative behavior of the effective medium as the effective medium transition rates are given by

As long as the mean transition rate is finite, the effective medium is exactly the deterministic long-range system of Sec. III with K = E[Z x,y ]. Recently, it was proven under certain regularity conditions that this is the correct self-averaging limit of the random walk in the random network [41] . We can draw the same conclusions for the random model as we did for the deterministic one, namely that the diameter of the infected regions grows exponentially, just like in Eq. (14) . Although, this behavior is known in the literature [35, 42, 43] , EMA opens a new way to analytically compute the speed of the infection spreading or even other quantities of desire. Importantly, the method presented here is not limited to the simple topology and the simple choice of transport rates that we used in Eq. (17) . In our example, the effective medium transport rates are simple averages of the original rates and the EMA result becomes equal to the annealed adjacency matrix approximation of Ref. [21] . This is however a consequence of the high connectivity and the power-law in Eq. (17) and doeas not have to hold in general. For more general topologies or other scaling relations, a different effective medium has to be chosen. This is already seen in the traditional EMA examples, e.g., a random short-range model (the so-called random barrier model, see e.g., [44] ), where only next-neighbor transport is allowed. The equation system Eq. (15) reduces to a single equation for the effective medium diffusivity K:

For a barrier model in one dimension d = 1, one finds K = E[1/W x,x+1 ] −1 . The effective medium diffusivity is given by the reciprocal of the harmonic mean, instead of by the arithmetic mean of the transport rates.

Since EMA reproduces the diffusive behavior of random systems pretty well [27] , it is a good candidate to produce a disorder-averaged random walk that can be used in Eq. (6) . Using EMA, one can make predictions about the reactiondiffusion system with a random transport operator, as we demonstrate numerically in the next section.

To validate our theory, we consider a ring of subpopulations with transport rates defined by Eq. (17). As mentioned above, the actual distribution of Z x,y does not matter, hence we sampled them uniformly from the interval [0,1]. Therefore, K = E[Z] = 0.5 in our simulations. It is important to note that the metric of a ring in one dimension is used: 14) is given in Fig. 2 . The numerical data respects the bound nicely.

Since the numerical diameter D(t ) shows a nice exponential growth pattern like in Eq. (14), we can extract some of the parameters from the exponential fit

Comparison with Eq. (14) would give measured values for α, β − μ and the diffusivity K = E[Z x,y ]. This may however be a hard task, because the non-linear term t B is not easy to detect in the exponential fit. The diameter's growth rate,

however, is easy to obtain, as it can also be measured from the slope of the tail of ln D(t ). In our simulations of the SI metapopulation model we obtain C = 0.076, which gives α fit = 1.622. This is a reasonably close value to the Lévy exponent α = 1.5 used for generating the graph realizations in the first place. D(t ) is only presented before the saturation sets in, and before the whole ring is infected. We now consider the SIS model in N = 8000 subpopulations with β = 0.2 and μ = 0.1, which gives a basic reproductive number of R 0 = 2. The correct time frame to assess D(t ) is visible in a prevalence plot; see Fig. 3 . In this figure the curves ρ x (t ) for each x are plotted against time; the stationary value ρ x (∞) = 1 − μ/β as well as the time when ρ x (t ) >c can be read from such a plot. For N = 8000 subpopulations the time gap between the outbreaks of the first and last subpopulation infected is 124 time steps, and the absolute global infection time is 193 time steps. As expected, the results are similar to the SI case. For the estimation of the Lévy exponent at this reproductive number we find C = 0.041, which results in α fit = 1.454, i.e., in only 3% error of the theoretical value.

Varying the reproductive number and measuring the growth rate C or the Lévy exponent α, respectively, leads to When the theoretical Lévy exponent α is varied and the growth rate is measured, the agreement appears much worse; see Fig. 4 

This mismatch is easily explained as pure finite-size effects as we show in Fig. 5 . There, we plotted the difference |C the − C fit | between the measured growth rate C fit and its theoretical prediction from Eq. (23) in a double logarithmic fashion against the system size N . The figure shows that the error decays at least like a power law and will vanish in the thermodynamic limit N → ∞. Due to the extreme long-range connections, ρ x (t ) saturates very quickly. This leads to very short time frame in which D(t ) grows exponentially that makes a correct estimation of C difficult. The effect becomes worse as α decreases, which also explains the slightly worse agreement for small α in Fig. 4(a) . For this reason we concentrated our numerical studies to the range α ∈ (1, 2).

We found in Fig. 4 (b) that C fit > C the in that data range, which should not be possible as Eq. (14) represents an upper

bound. An overview of the agreement with the theoretical bound for the probed range in α is shown in Fig. 6 . We find that for some of these large values of α the numerical data overestimate the EMA bound. This, however, only happens in an intermediate time regime and not in the long time limit, in which we derived Eq. (14) . In fact, we find that the upper bound is respected in the long-time limit for all values of α. The predictions given by EMA are rigorously valid in the thermodynamic limit N → ∞, when the infection propagates indefinitely and saturation is never reached.

The goal of this paper was to present a new analytical tool for reaction-diffusion problems in random long-range networks. 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 that is well respected in the long time limit of our numerical simulations. We also demonstrated that certain parameters, like the Lévy exponent α, can be extracted from data, thus verifying that a made assumption on the random network's ensemble is correct. This way we demonstrated that EMA is still relevant even beyond the shortrange connection paradigm.

With the human travel network in mind, we presented a simple metapopulation model with random long-range connections. We reproduced the exponential growth of the infection diameter, that is known in the literature [35, 42, 43] . Our EMA prediction of 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. (17) . Other characteristics of the transition rates (like their mean) only play a minor role in the dynamics. Notice, that long-range links with a "weak" power law-i.e., α > 2 in Eq. (17)-would eventually lead to a ballistic growth of the infection front. These results are also discussed in Ref. [12] .

The main restrictions of EMA are currently the necessity of independent and symmetric transport rates W x,y = W y,x . Future modifications of EMA are necessary to deal with asymmetric rates, and can thus take variable subpopulation sizes into account. Furthermore, when it is possible to deal with correlated links, more realistic models than a simple grid of the subpopulation's locations can be included.

In its current form, EMA could already be used to tackle more involved models than the one considered here. For example, an extension of our argument to d dimensions is possible without major change and would only lead to a different growth rate of C = (β − μ)/(d + α). Internal dynamics on the nodes (like commuting agents) could be considered by replacing the subpopulations by small networks themselves. The EMA method is not restricted to the simple model considered here. In particular, one can overcome the strong finite-size effects, that we encountered in our work by considering a finite-size effective medium instead of an infinite one, as we did here for simplicity.

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. Note that the case discussed in Ref. [45] violates both assumptions of absence of correlations and deviation from percolation transition and leads to a very different behavior termed paradoxical diffusion.

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. We believe that EMA will develop to a great practical tool for the analysis of dynamics on networks.

Network Science

Synchronization: A Universal Concept in Nonlinear Sciences

Anomalous Transport: Foundations and Applications

Effective Medium Theory-Principles and Applications, International Series of Monographs on Physics

Dynamical Processes on Complex Networks

Selected Works of A. N. Kolmogorov

An Introduction to Probability Theory and Its Applications

Graphs and Matrices

The authors are indebted to A. Vulpiani for insightful discussions. This work has partially been funded by the DFG/FAPESP, within the scope of the IRTG 1740/TRP 2015/50122-0.

In this Appendix we show how to obtain Eq. (13) from Eq. (12) . The computation follows closely Ref. [27] . We start by plugging the transport rates into Eq. (5) and using their symmetry to reorder the summation:The equation is solved using Fourier transform, i.e., we multiply e ikx on both sides and sum overall x. Defining P (k; t ) = x∈Z e ikx P x (t ), we obtaiṅ (1)].Here Li ν (z) = ∞ n=1 z n /n ν is the polylogarithm function and S(k) is the Fourier symbol of the transport operator defined asUsing the polylogarithm's expansion around k = 0 (obtained by Mathematica) one finds the following small wave-vector expression for S(k):Note that the sign of S(k) is negative for all α ∈ (0, 2). As we now have˙ P (k; t ) = S(k) P (k; t ), the solution is given by P (k; t ) = exp[S(k)t] ∼ exp(−a|k| α ), where we used the initial condition P x (t = 0) = δ x,0 , which gives P (k; t = 0) = 1 and identifiedThe expansion also shows that P x (t ) 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 |x|, [39] :Here F −1 denotes the inverse Fourier transform.Using this equation and the definition of a Eq. (A4), with (α + 1) = α (α), 2 sin(y) cos(y) = sin(2y) and (y) (−y) = π/ sin(πy), we recover Eq. (13) from the main text.