# Fluctuation effects in metapopulation models: percolation and pandemic threshold

## Abstract

Metapopulation models provide the theoretical framework for describing disease spread between different populations connected by a network. In particular, these models are at the basis of most simulations of pandemic spread. They are usually studied at the mean-field level by neglecting fluctuations. Here we include fluctuations in the models by adopting fully stochastic descriptions of the corresponding processes. This level of description allows to address analytically, in the SIS and SIR cases, problems such as the existence and the calculation of an effective threshold for the spread of a disease at a global level. We show that the possibility of the spread at the global level is described in terms of (bond) percolation on the network. This mapping enables us to give an estimate (lower bound) for the pandemic threshold in the SIR case for all values of the model parameters and for all possible networks.

###### keywords:

epidemic spread, metapopulation, percolation, pandemic threshold^{1}

## 1 Introduction

Modeling the spread of a disease between different populations connected by a transportation network is nowadays of crucial importance. Global spread of pandemic influenza is an illustrative example of the importance of this problem for our modern societies. In the context of pandemic spread and airline transportation networks, Rvachev and Longini Rvachev:1985 () proposed a model, now called a metapopulation model, which was used to describe the spread of pandemic influenza Longini:1988 (); Grais:2003 (); Grais:2004 (); Colizza:2006a (), SARS Hufnagel:2004 (); Colizza:2007 (), HIV Flahault:1991 (), and very recently swine flu Vespi:2009a (); Vespi:2009b (). This model describes homogeneously mixed populations connected by a transportation network. Despite its successes, very few theoretical studies are available. Numerical studies Colizza:2006b () showed that weight heterogeneity played a major role in the predictability by creating epidemic pathways which compensate the various scenario possibilities due to hubs Colizza:2007 ().

A major advance in our understanding of metapopulation models was done in the series of papers by Colizza, Pastor-Satorras and Vespignani Colizza:2007a (); Colizza:2007b (); Colizza:2008 (). In particular, these authors studied the condition at which a disease which can spread in an isolated population (and has therefore a basic reproductive rate (see below)) can invade the whole network. In this respect, these authors defined an effective reproductive rate at the global level and showed that its expression for a random network, in the SIR case with , and with a travel rate , is under a mean-field approximation Colizza:2007a ()

(1) |

where the brackets denote the average over the configurations of the network, while represents the fraction of infected individuals in a given city during the epidemic, is the mean time an individual stays infected, and is the mean population of a city. For the disease spreads over an extensive number of cities. For scale-free networks, the second moment being very large, travel restrictions (i.e., decreasing ) have a very mild effect. Equation (1) relies on a number of assumptions: the network is uncorrelated and tree-like, is close to , and more importantly the number of infected individuals going from one city to another is assumed to be given by the constant . This last expression however neglects different temporal effects. Indeed, in the SIR case we can distinguish two different time regimes: the first regime of exponential growth corresponding to the outbreak, and the recovery regime where, after having reached a maximum, the number of infected individuals falls off exponentially. The number of infected individuals going from one city to another thus depends on different factors and can be very different from the expression for given above.

In the present work, we consider the stochastic version of the metapopulation model. We investigate in detail the number of infected individuals going from one site to the other. This allows us to address the problem of the condition for propagation at a global scale without resorting to mean-field approximations or . We will show that the latter can be recast in terms of a bond percolation problem.

So far most of the simulations on the metapopulation model were done in essentially two ways. The first approach consists in using Langevin-type equations, with a noise term whose amplitude is proportional to the square root of the reaction term Gardiner:2004 (); Colizza:2006a (). In this approach, the evolution of the various populations (infected, susceptibles, removed) is described by a set of finite difference equations, with noise and with traveling random terms describing the random jumps between cities Colizza:2006a (). This level of description suffers from several technical difficulties when the finite-difference equations are iterated: truncation of noise, non-integer number of individuals, etc., which render its numerical implementation difficult. Another possibility consists in adopting an agent-based approach, i.e., simulating the motion of each individual Colizza:2006a (). This simulation is free of ambiguities but requires a lot of memory and CPU time.

In contrast, in the present work we will put forward an intermediate level of description, which is free from the technical difficulties of the Langevin-type approach (it does not require to iterate finite-difference equations), and which does not require large amounts of CPU. This approach consists in describing the process at the population level, considering that individuals are indistinguishable. In Statistical Mechanics, this level of description is used for urn models or migration processes (see Godreche:2002 (); Godreche:2007 () for reviews). Let us take for definiteness the example of the stochastic SIS epidemic process. For a single isolated city with population , this process (see for example Allen:2003 ()) is described in terms of two variables: (number of susceptible individuals) and (number of infected individuals), which evolve in continuous time as

(2) |

where is conserved, and with initial condition . The first reaction corresponds to the meeting of a susceptible individual with an infected individual, resulting in two infected ones. The second reaction corresponds to the recovery of an infected individual into a susceptible one. The SIR case has an analogous definition (see Section 3).

For two or more cities, we have to take into account the process of traveling between different cities. This diffusion process is described at the microscopic level by the jumping of individuals from one city to another. Let and be two neighboring cities. As described above, at the level of populations (i.e., such that the individuals are indistinguishable), the process inside each city is given by (2), while the elementary traveling processes between these two cities are:

(3) |

where is the probability per unit time and per individual to jump from city to city .

As said above, the numerical implementation of these processes is easy, and it allows us to make simulations with large numbers of large cities, in contrast with agent-based simulations. In the following we will restrict ourselves to uniform symmetric diffusion, defined by the constant rates if cities and are neighbors. These rules are compatible with a stationary state where all cities have the same population () on average. We will consider the initial condition , where is the initial number of infected individuals, which all live in city 0 at time . All the numerical simulations of the SIS and SIR cases will be performed with and , so that (see (6)).

The outline of the paper is as follows. In Section 2 we consider the SIS case in different geometries, starting with the case of a single city, then proceeding to the case of two cities, and finally discussing the propagation on an infinite one-dimensional array of cities. As the number of infected individuals in a given city saturates to a non-zero fraction of the total population, the metapopulation model will always experience a pandemic spread and ballistic propagation. It is however instructive to study this model for the following reasons: (i) it displays the same phase of exponential growth as the SIR model, but with the great simplification that only two compartments ( and ) are present, allowing analytical approaches; (ii) the SIS model has connections with the FKPP equation Fisher:1937 (); Kolmo:1937 (), which is still currently the subject of intense activity. In Section 3 we consider the SIR case along the same lines. We successively consider the simple cases of one and two cities, before addressing the question of the propagation on various extended structures (one-dimensional chain, square lattice, Cayley tree, uncorrelated scale-free network). The key difference with the SIS case is that the time integral of the number of infected individuals in a given city is now finite. We demonstrate that the latter property implies the existence of an effective pandemic threshold for the spread of the disease over the whole network. We also obtain an estimate (lower bound) for this pandemic threshold as a function of all the model parameters and of the bond percolation threshold on the underlying geometrical structure. For scale-free networks our prediction is virtually identical to that of Colizza:2007a ().

## 2 Metapopulation model in the SIS case

In this section we discuss the metapopulation model in its stochastic form, on the example of the SIS process. We start with the case of a single city, then proceed to the case of two cities, which will be the building block of our analysis of the propagation of the epidemic on extended structures. We finally discuss the ballistic propagation on a one-dimensional array of cities.

### 2.1 Stochastic SIS model for one city

We first present a discussion of the stochastic version of the SIS model for one city.

Let be the population of the city and the number of infected individuals at time . As can be seen from the defining reactions (2), the SIS process can be described by the single random variable , which experiences a biased continuous-time random walk on the interval , with variable rates which only depend on the instantaneous position of the walker. This process is known as a birth-and-death process in the probabilistic literature feller (). Let us denote by and the jump rates of the walker, respectively to the right () and to the left (). Here we have

(4) |

Hence the origin () is absorbing while the right end () is reflecting. The walker moves in an effective confining potential, with a restoring force bringing it back to the quasi-equilibrium position corresponding to the equality of the rates and . This position can be interpreted as the number of infected individuals at saturation. It reads

(5) |

The so-called basic reproductive number is defined as

(6) |

and we will hereafter focus on the case of interest where is larger than 1, which is the condition for the occurrence of an outbreak at the level of a single city. The quasi-equilibrium state around is however metastable, since the walker is deemed to be ultimately absorbed at the origin, which corresponds to the extinction of the epidemic in the city. The characteristic absorption time is however exponentially increasing with the population of the city, as will be shown below.

Let us denote by the probability for the number of infected individuals to be equal to the integer at time . This probability contains all the information on the one-time properties of the process (see Godreche:2002 (); Godreche:2007 () for reviews). Its temporal evolution is given by the following master equations, characteristic of birth-and-death processes karlin (); kendall ():

(7) |

The nonlinear dependence of the rates on the position renders their analytical study difficult. For example the evolution equation of the first moment reads

(8) |

whereas the equation for itself involves and so on. However, eqs. (2.1) are easily implemented numerically, for a given initial condition. In Figure 1 (top) we show a plot of obtained by a numerical integration of eqs. (2.1). We also plot the same quantity as obtained from a numerical simulation of the SIS process. The data were obtained for , , , and , discarding histories leading to absorption. The two sets of data are indistinguishable one from the other. Finally we also show for comparison the deterministic (or mean-field) expression of (see (10) below). Figure 1 (bottom) shows five different realizations of the stochastic process for a single city, together with the deterministic solution (10).

The deterministic (or mean-field) approximation consists in neglecting any correlations. For example, in (8), it amounts to replacing by . This approach is a priori legitimate if the population of the city is large, so that the number of infected individuals is also typically large enough to be considered as a deterministic observable, instead of a random variable. This approximation thus reduces the SIS model for one city to a deterministic dynamical system. Then (8) becomes the following dynamical equation for the temporal evolution of :

(9) |

the solution of which can be cast into the form

(10) |

As demonstrated by Figure 1, the model exhibits an exponential growth regime followed by a saturation regime. The deterministic approximation gives an accurate global description of , even for rather small populations. It however misses the effect of fluctuations in the process, that we now study in both regimes successively.

#### Exponential growth regime

The regime of exponential growth formally corresponds to taking the infinite limit at fixed time . In this limit, (8) and (9) simplify, and their common solution reads

(11) |

Furthermore, the rates become linear in , i.e., and , reducing the model to a solvable birth-and-death process feller (); kendall (). The master equations for the read

(12) |

The generating function

(13) |

satisfies

(14) |

the solution of which is obtained with the method of characteristics and is

(15) |

In particular, the probability for the number of infected individuals to be zero at time , , reads kendall ()

(16) |

For we recover the well-known extinction probability Bailey ()

(17) |

If this extinction probability is equal to 1. Hence the condition ensures the possibility of a local outbreak at the level of a single city. The expressions of the moments of can be extracted from (15). We thus recover (11), while the variance reads

(18) |

In the late stages of the exponential growth regime, the relative variance is proportional to :

(19) |

#### Saturated regime

After the phase of exponential growth is over, the mean number of infected individuals saturates at the value given by (5), as seen in Figure 1. The time to reach saturation therefore scales as .

The saturated state is actually a metastable state, rather than a genuine stationary state. For finite , the system will indeed always end in its absorbing state at . The distribution of the fluctuating number of infected individuals in the metastable state can be evaluated by cutting the link from to which is responsible for absorption. The process thus modified reaches an equilibrium state at long times. In the latter state the distribution satisfies detailed balance, that is

(20) |

where and are defined in (4). This equation yields

(21) |

where the large deviation function reads

(22) |

As expected, the distribution is peaked around for large. Indeed the function takes its maximal value,

(23) |

for . We therefore have the exponential estimate . The lifetime of the metastable state can then be estimated, in the spirit of the Arrhenius law, as . It is therefore predicted to grow exponentially with the population as

(24) |

For we have .

The expression (21) also yields an estimate for the Gaussian fluctuations of around its value in the saturated state. Indeed, expanding to second order around , we obtain the estimate

(25) |

for the variance of the distribution of in the quasi-stationary state, so that the reduced variance

(26) |

is proportional to , whereas it was proportional to in the growth phase (see (19)). These two results provide a quantitative confirmation that relative fluctuations around the deterministic theory become negligible in all regimes, as soon as the number of infected individuals is large. Figure 2 shows a plot of , obtained by integration of the master equations (2.1). The data for small and large times are found to be in very good agreement with the estimates (18) and (25), respectively.

### 2.2 Including travel: two cities

Consider now two cities, numbered 0 and 1. The traveling process is described by (3). We assume symmetric diffusion: , hence the system reaches a stationary state such that . We choose the initial condition . We are primarily interested in the time of occurrence of the outbreak of the epidemic in city 1. This event is governed by the arrival of infected individuals traveling from city 0 to city 1. Returns of infected individuals from city 1 to city 0 can be neglected throughout.

We start by analyzing the distribution of the first arrival time of an infected individual in city 1. Let us denote by the probability that no infected individual exits from city 0 in the time interval , with in particular . We have

(27) | |||||

where the expression in the last parentheses is the probability that no infected individual exits from city 0 in the infinitesimal time interval . We thus obtain

(28) |

and so

(29) |

with

(30) |

An alternate derivation of this result Eqs. (LABEL:eq:arrivaltime,30) is given in Gautreau:2008 (). We thus conclude that the number of arrivals of infected individuals in city 1 in the time interval is a Poisson process, for which the rate is itself a stochastic process, equal to , as the integrated rate is . Therefore

(31) |

This generalization of the Poisson process is known in the literature as the Cox process Cox (). Hereafter, in order to simplify the presentation, we will implicitly assume that probabilities are conditioned on a single stochastic history .

The mean first arrival time can easily be derived from (LABEL:eq:arrivaltime). It has the simple form

(32) |

This exact expression can not be written in closed form, even if is given its deterministic value (10). The scaling properties of at large can however be obtained using the following simple arguments. In the large regime, we expect , as the number of infected individuals has hardly changed from its initial value in time . In the opposite regime, where is very small, we have (see (5)), and thus . In the intermediate regime, city 0 is in its phase of exponential growth, and the mean arrival time can be estimated by imposing that is of order unity. We thus have

(33) |

The crossovers between these three regimes occur at and . Figure 3 summarizes the above discussion, whereas actual data will be shown in Figure 5.

We now turn to the distribution of the outbreak time , defined as the time of the arrival in city 1 of the first infected individual who induces an outbreak in city 1. If the travel rate is small enough, events where two or more infected individuals coming from city 0 are simultaneously present in city 1 can be neglected. Each infected individual entering city 1 has therefore a chance (see (17)) to disappear before it induces an outbreak. On the other hand, the distribution of the arrival time of the -th infected individual reads

(34) |

The probability that no outbreak occurred up to time is thus given by

(35) | |||||

Comparing this expression with the corresponding one for the first arrival (LABEL:eq:arrivaltime) reveals that taking into account extinction amounts to renormalizing to an effective rate , and accordingly to , multiplying both quantities by the scaling factor :

(36) |

Within this approach, the distribution of the outbreak time in city 1 reads

(37) |

### 2.3 Ballistic propagation on a one-dimensional array

We now consider the stochastic SIS model in the one-dimensional geometry of an infinite array of cities with initial populations . We still choose a symmetric travel rate , hence, in the course of time, these cities remain equally populated on average.

The model is observed to reach very quickly a ballistic propagation regime, where the epidemic invades the whole array by propagating a ballistic front moving at a well-defined finite velocity . Such a front is illustrated by the space-time plot of Figure 4.

The velocity of the front can be estimated as follows. The time at which the epidemic reaches city can be recast as

(38) |

where is the outbreak time in city , i.e., the arrival time of the infected individual which will trigger the outbreak in that city, with the origin of times being set to the outbreak time in the previous city . Modeling the propagation from city to city by the case of two cities studied above, we thus predict that the inverse velocity reads

(39) |

where is the mean outbreak time of the problem of two cities, with the natural initial condition .

In Figure 5 we show a comparison between data for the reciprocal of the ballistic velocity and for , the mean outbreak time for two cities, with . The latter was measured in a Monte Carlo simulation as with renormalized to (red symbols), and computed using the analytic prediction (32), where has the deterministic expression (10), again with renormalized (blue curve). The agreement validates the approximation (39). The prediction (43), (44) of the discrete FKPP theory is also shown.

It is indeed worth confronting the above analysis with yet another approach, based on the analysis of the following deterministic equation for the densities of infected individuals :

(40) |

This equation is the discrete version of the FKPP equation Fisher:1937 (); Kolmo:1937 (). Looking for a traveling-wave solution of the form , propagating with velocity , it follows that the function obeys the differential-difference equation

(41) | |||||

If one assumes that the solution decays exponentially as for , the velocity is found to depend continuously on , according to

(42) |

For a localized initial condition, the actual velocity of the front is known bramson (); derridaspohn () to be obtained by minimizing the above expression with respect to the spatial decay rate . This velocity reads

(43) |

and it is reached for

(44) |

Both above equations give a parametric expression for the ballistic velocity . The usual prediction of the continuum FKPP equation, i.e., , is recovered as , i.e., . The regime of most interest in the present context is the opposite one (, i.e., ), where we have

(45) |

This regime is in correspondence with the estimate (33) which holds in the intermediate growth regime. Figure 5 shows that the discrete FKPP theory provides a good, albeit not quantitative, description of the overall behavior of the ballistic velocity.

To conclude, in the one-dimensional SIS model, the epidemic always spreads ballistically, with a non-zero velocity . We expect this result to hold on any extended structure, since it is merely a consequence of the existence of a very long-lived saturated regime where a finite fraction of individuals are infected. In a metapopulation model, this will for sure trigger an infection which will spread over the whole network as long as .

## 3 Metapopulation model in the SIR case

We now turn to the case of the SIR model for the growth of an epidemic in a given city. In contrast with the SIS model, the number of infected individuals in a given city now falls off to zero at large times. As a consequence, if the traveling rate is too small compared to the typical inverse duration of the epidemic phase, the epidemic may die out in a city before propagating to the neighboring ones. This phenomenon will lead to the existence of a non-trivial threshold for the travel rate, below which the disease will not spread in the system. We will show that the problem can be recast in terms of a bond percolation problem. This picture greatly simplifies the analysis, and leads to an estimate (lower bound) of the pandemic threshold.

### 3.1 Stochastic SIR model for one city

We first investigate the stochastic SIR process in the case of one isolated city of population . As in the SIS case, the process is described at the population level (i.e., the individuals are indistinguishable), but we now have three variables , , and (number of recovered individuals), with . The stochastic process is described by the reactions

(46) |

with initial condition and . The second reaction corresponds to the recovery of an infected individual, that stays immune later on. The master equations for the process can be written in terms of two independent integer random variables, say and , and integrated numerically.

The deterministic equations describing the SIR process read

(47) |

At variance with the case of SIS, the above equations cannot be solved in closed form. Figure 6 shows a comparison between the numerical solution of these equations and a numerical simulation of the stochastic process. The upper panel shows five different trajectories in the - plane. We observe relatively important fluctuations around the deterministic solution. In particular the epidemic stops after a finite random time. The lower panel demonstrates that mean quantities are however well described by the deterministic approach.

A quantity of central interest for the sequel is the final size of the epidemic, defined as the number of recovered individuals after the epidemic has stopped:

(48) |

We refer the reader to Bailey:1953 (); Ball:1994 (); Loef:1998 () for studies on the distribution of this quantity. For a city with a large population (), can be approximately determined in the framework of the deterministic approach. It follows from (47) that

(49) |

from which we obtain

(50) |

Using the initial condition given above, and the final condition , i.e., , we obtain a transcendental equation relating the densities and :

(51) |

which reduces for , i.e., , to

(52) |

The density starts rising linearly as

(53) |

and reaches unity exponentially fast for , as . For , the value used in numerical simulations, we have .

### 3.2 Including travel: two cities

As in the SIS case, it is natural to start the study of propagation by the case of two neighboring cities. The SIR process inside each city is described in terms of three random variables , , and (), with , which evolve stochastically as in (46). The traveling process between these two cities is given by the reactions

(54) |

where stands for , or . We again assume symmetric diffusion. A typical history of the system is shown in Figure 7.

The first arrival time of an infected individual in city 1 is still distributed according to (LABEL:eq:arrivaltime) and (30). There is however a crucial difference between SIS and SIR. As already discussed, in the present case of SIR, the integrated rate converges to a finite limit as . For a large population , the deterministic approach yields

(55) |

As a consequence, there is now a non-zero probability that no infected individual travels from city 0 to city 1 during the whole epidemic in city 0. In other words, the event that at least one infected individual reaches city 1, i.e., that the time is finite, only occurs with probability

(56) |

The distribution of is said to be defective feller ().

Taking into account the fact that single infected individuals only trigger an outbreak with probability , the distribution of the outbreak time in city 1 is still given by (LABEL:tout), as a result of the renormalization procedure which led us to replace by the effective rate (see (36)). In particular, the probability of occurrence of an outbreak in city 1 reads

(57) |

This probability is non-trivial, i.e., less than one, in contrast with the SIS case.

### 3.3 Propagation on extended structures

We can now consider the general case where the cities are connected so as to form a network or any other kind of extended structure (one-dimensional array, regular finite-dimensional lattice, regular tree). Individuals can travel by performing diffusion along the links of the network, allowing thus the disease to spread over different cities.

The probability given in (57) can be interpreted as the probability that the disease propagates through one given link from a city to one of its neighbors. In the SIS case, since the integrated rate diverges with time, the probability is trivially equal to unity. In contrast, for the SIR process, this probability is less than unity, so that the disease can stop invading the network.

For a given network, denoting the bond percolation threshold by , the condition for the disease to spread is therefore

(58) |

This equation defines the pandemic threshold: for large, (58) yields , where

(59) |

This static prediction is not claimed to give an exact value for the pandemic threshold. It can however be argued to provide a lower bound for the threshold, which is also meant as a reasonable and useful estimate. Indeed the picture of static percolation, where links are occupied with constant probability , corresponds to an ideal infinitely slow propagation, where the epidemic can take an infinitely long time to cross some of the links. In real situations, propagation is rather observed to take place with finite velocity. This velocity can be thought of as providing a time cutoff for propagation across every single link of the network. The effect of this cutoff time is to decrease the integrated rate from to the smaller value , and hence to increase the threshold value of from the theoretical prediction (59) to a higher value.

Let us now illustrate this result on some examples of networks and other extended structures.

#### Cayley tree

We consider first the geometry of a regular Cayley tree (or Bethe lattice) with coordination number , so that .

For the simplest example of , and an epidemic initially located at the center of a tree with 15 generations, Figure 8 shows a plot of the fraction of infected cities at a very long time. We indeed observe a pandemic threshold near . For , , , hence , and , the predicted static threshold is . The observed threshold is thus some two times larger than the static one.

#### Square lattice

We now consider the two-dimensional situation of propagation on the square lattice. This lattice also has bond percolation threshold , so that the above static threshold of still holds for the same parameter values.

On a array, the data of Figure 9 clearly demonstrates the existence of a threshold near , in good quantitative agreement (within ten percent, say) with the predicted static value. In figure 10 we provide an example of the observed shape of the infected region for a large but finite time.

#### One-dimensional array

We finally address the case of propagation along an infinite array of cities. The percolation threshold in the one-dimensional case is . Our analysis therefore suggests that the disease will always stop, after having invaded only a finite range of typical size , but not the whole array. The static percolation approach suggests that diverges as in the same way as the static correlation length in percolation theory, namely . We thus find the exponential growth

(60) |

We indeed observe a symmetric ballistic front on both sides of the seed. Each branch suddenly dies at a random time, and hence in a randomly located city, say number . Such a front is shown in Figure 11. For given parameter values, the right stopping point is observed to be exponentially distributed, with a cumulative distribution falling off as (see Figure 12). This exponential distribution confirms our intuition that a local mechanism is responsible for the arrest of the propagation. The and dependence of the characteristic length thus measured is shown in Figure 13. The data demonstrate an increase of with both and . The agreement with the static prediction however remains very qualitative.

#### Scale-Free networks

We close up our analysis of the propagation on a SIR epidemic by considering the case of a scale-free network. For an uncorrelated network, the percolation threshold is known to be given by (see (69)). As most scale-free networks are such that , we have , so that our static prediction (59) for the pandemic threshold reads

(61) |

In the same regime, it is worth recasting the result of Colizza:2007a (), where the expression for is recalled in (1), in the form of a threshold value of the travel rate. Using the fact that for SIR, we obtain

(62) |

Our result thus only differs from that of Colizza:2007a () by an inessential multiplicative factor involving the mean degree and the reproductive rate . The coincidence is all the more striking that both approaches are largely different. Both expressions agree to predict that the threshold diverges quadratically as as the basic reproductive rate . Indeed vanishes linearly (see (53)).

## 4 Conclusion

In this paper we have considered the metapopulation models for the SIS and SIR processes in their stochastic versions. We have put forward a description considering individuals as indistinguishable, but keeping the full stochastic character of the populations in the various compartments. This intermediate level of description, along the lines of the theory of urn models and migration processes, allows very efficient numerical simulations. We have considered with great care temporal effects in the spread of a disease from one city to another and investigated its propagation at a global scale on scale-free networks and other extended structures such as regular trees and lattices.

For the SIS case we always observe a ballistic propagation. For the SIR case, even after an infinitely long time, the disease only gets transmitted with probability through every link of the network. This picture allows us to rephrase the possibility of global invasion as a static bond percolation problem. The resulting prediction is an estimate (lower bound) for the pandemic threshold, expressed as a threshold value for the travel rate . In the case of scale-free networks, our prediction is virtually identical to that of Colizza:2007a (). Our approach is however not limited to the case , and takes into account both temporal and topological fluctuations. For other geometries (Cayley tree, square lattice), our static prediction yields a reasonable estimate for the pandemic thresholds observed in numerical simulations.

As in most previous theoretical studies, we assumed that the populations of all cities were equal and that the traveling rate per link and per individual was constant. It would be interesting to extend our result and to test for the relevance of travel and/or population heterogeneities on the existence of a pandemic threshold.

Finally, it is worth putting the present work in a broader perspective. The relationship between epidemic spreading and percolation has already been discussed in the context of contact networks. Table 1 summarizes the parallel between the contact network and the metapopulation approaches. It has been emphasized first by Grassberger that the spreading of an epidemic over a contact network can be mapped onto directed percolation Grass:1983 (). It was then realized that, as a general rule, epidemic models without immunization belong to the universality class of directed percolation, whereas those with immunization, such as SIR, belong to the universality class of dynamical percolation (see hinrichsen () for a review). In the latter case, clusters of immune individuals have the same critical properties as usual percolation clusters. The effective description of epidemic spreading put forward in this work demonstrates that epidemic spreading is also intimately related to static percolation at the metapopulation level.

Network | Nodes | Links |
---|---|---|

contact | individuals | contacts |

metapopulation | cities | travel |

## Appendix A Pandemic threshold on a network with arbitrary travel rates

In this appendix we determine the pandemic threshold on an uncorrelated network where travel is described by arbitrary hopping rates.

We consider an uncorrelated network, denote by be the degree distribution of the nodes, and use the framework of degree classes, assuming that all nodes with given degree are equivalent. Travel is defined by the rate for an individual to hop from a node with degree to a neighboring node of degree . The stationary populations obey the balance equation

(63) |

where

(64) |

is the degree distribution of a neighboring node of a given node. For instance, uniform rates yield uniform populations , whereas the rates of ordinary random walk yield populations . More generally, separable (i.e., factorized) travel rates of the form yield stationary populations .

The basic quantity in the metapopulation model is the probability that the epidemic will propagate (in an infinitely long time) from a node of degree to a node of degree . Dropping the star for conciseness, this probability reads

(65) |

In order to determine the static pandemic threshold of SIR on the network, we are thus led to consider the problem of directed bond percolation on the network defined by the above probabilities that an oriented link is open from a node of degree to a node of degree . Let us introduce the probabilities that an oriented link from a node of degree to a node of degree leads to a node belonging to the giant component. These quantities obey

(66) |

where runs over the neighbors of the node of degree which are not the initial one and are their degrees. Near the threshold, the probabilities are expected to be small. Linearizing the above relation, we obtain

(67) |

Setting , the above condition reads , with

(68) |

The percolation threshold is therefore given by the condition that the largest (Perron-Frobenius) eigenvalue of the positive matrix equals unity.

The usual percolation problem corresponds to the case where every link is occupied with probability , i.e., for all and . We thus recover the known result

(69) |

In particular, for a regular Cayley tree with coordination number , we have

(70) |

The reasoning leading to (66) can be extended to other cases. Let us mention the example of a bipartite tree, where nodes with degree are neighbors of nodes of degree . We then obtain

(71) |

For separable probabilities, of the form , we obtain the following threshold condition

(72) |

In a more general case, we have to resort to a numerical computation of the largest eigenvalue of the matrix .

### Footnotes

- journal: Journal of Theoretical Biology

### References

- L.A. Rvachev and I.M. Longini, A mathematical model for the global spread of Influenza, Math. Biosci. 75, 3-23 (1985).
- I.M. Longini, A mathematical model for predicting the geographic spread of new infectious agents, Math. Biosci. 90, 367-383 (1988).
- R.F. Grais, J. Hugh Ellis, and G.E. Glass, Assessing the impact of airline travel on the geographic spread of pandemic influenza, Eur. J. Epidemiol. 18, 1065-1072 (2003).
- R.F. Grais, J. Hugh Ellis, A. Kress, and G.E. Glass, Modeling the spread of annual influenza epidemics in the US: the potential role of air travel Health Care Manage. Sci. 7, 127-134 (2004).
- V. Colizza, A. Barrat, M. Barthélemy, and A. Vespignani, The Modeling of Global Epidemics: Stochastic Dynamics and Predictability, Bull. Math. Biol. 68, 1893-1921 (2006).
- L. Hufnagel, D. Brockmann, and T. Geisel, Forecast and control of epidemics in a globalized world, Proc. Natl Acad. Sci. (USA) 101, 15124-15129 (2004).
- V. Colizza, A. Barrat, M. Barthélemy, and A. Vespignani, Predictability and epidemic pathways in global outbreaks of infectious diseases: the SARS case study, BMC Medicine 5, 34 (2007).
- A. Flahault and A.-J. Valleron, A method for assessing the global spread of HIV-1 infection based on air-travel, Math. Pop. Studies 3, 1-11 (1991).
- P. Bajardi, C. Poletto, D. Balcan, H. Hu, B. Goncalves, J.J. Ramasco, D. Paolotti, N. Perra, M. Tizzoni, W. Van den Broeck, V. Colizza, and A. Vespignani, Modeling vaccination campaigns and the Fall/Winter activity of the new A(H1N1) influenza in the Northern Hemisphere, Emerging Health Threats Journal 2, 11 (2009).
- D. Balcan, H. Hu, B. Goncalves, P. Bajardi, C. Poletto, J.J. Ramasco, D. Paolotti, N. Perra, M. Tizzoni, W. Van den Broeck, V. Colizza, and A. Vespignani, Seasonal transmission potential and activity peaks of the new influenza A(H1N1): a Monte Carlo likelihood analysis based on human mobility, BMC Medicine 7, 45 (2009).
- V. Colizza, A. Barrat, M. Barthélemy, and A. Vespignani, The role of the airline transportation network in the prediction and predictability of global epidemics, Proc. Natl Acad. Sci. (USA) 103, 2015-2020 (2006).
- V. Colizza and A. Vespignani, Invasion threshold in heterogeneous metapopulation networks, Phys. Rev. Lett. 99, 148701 (2007).
- V. Colizza and R. Pastor-Satorras, Reaction-diffusion processes and metapopulation models in heterogeneous networks, Nature Phys. 3, 276-282 (2007).
- V. Colizza and A. Vespignani, Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations, J. Theor. Biol. 251, 450-467 (2008).
- W.C. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry, and Natural Sciences, 3rd edition (Springer, New York, 2004).
- C. Godrèche and J.M. Luck, Nonequilibrium dynamics of urn models, J. Phys.: Condens. Matter 14, 1601-1615 (2002).
- C. Godrèche, From Urn Models to Zero-Range Processes: Statics and Dynamics, in Ageing and the Glass Transition (Lecture Notes in Physics 716) (Springer, Berlin, 2007), p. 261.
- L.J.S. Allen, An Introduction to Stochastic Processes with Applications to Biology (Pearson, Prentice Hall, New Jersey, 2003).
- R.A. Fisher, The wave of advance of advantageous genes, Ann. Eugenics 7, 353-369 (1937).
- A. Kolmogorov, I. Petrovsky, and N. Piscounov, Etude de l’équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique, Moscow Univ. Bull. Math. 1, 1 (1937).
- W. Feller, An introduction to probability theory and its applications (Wiley and Sons, New York, 1968).
- S. Karlin and H.M. Taylor, A First Course in Stochastic Processes (Academic Press, 1975).
- D.G. Kendall, On the generalized birth-and-death process, Annals of Math. Stat. 19, 1-15 (1948).
- N.T.J. Bailey, The Mathematical Theory of Infectious Diseases and its Applications, 2nd ed. (Griffin, London, 1975).
- A. Gautreau, A. Barrat, and M. Barthélemy, Global disease spread: statistics and estimation of arrival times, J. Theor. Biol. 251, 509-522 (2008).
- D.R. Cox, Some statistical methods connected with series of events, J. R. Statist. Soc. Ser. B 17, 129-164 (1955).
- M.D. Bramson, Convergence of solutions of the Kolmogorov equation to travelling waves, Mem. Amer. Math. Soc. 285, 1-190 (1983).
- B. Derrida and H. Spohn, Polymers on Disordered Trees, Spin Glasses, and Traveling Waves, J. Stat. Phys. 51, 817-840 (1988).
- N.T.J. Bailey, The total size of a general stochastic epidemic, Biometrika 40, 177-185 (1953).
- F. Ball and I. Nasell, The shape of the size distribution of an epidemic in a finite population, Math. Biosci. 123, 167-181 (1994).
- A. Martin-Loef, The final size of a nearly critical epidemic, and the first passage time of a Wiener process to a parabolic barrier, J. Appl. Prob. 35, 671-682 (1998).
- P. Grassberger, On the critical behavior of the general epidemic process and dynamical percolation, Math. Biosc. 63, 157-172 (1983).
- H. Hinrichsen, Non-equilibrium critical phenomena and phase transitions into absorbing states, Adv. Phys. 49, 815-958 (2000).