Continuous-time formulation of reaction-diffusion processes on heterogeneous metapopulations

We present the derivation of the continuous-time equations governing the limit dynamics of discrete-time reaction-diffusion processes defined on heterogeneous metapopulations. We show that, when a rigorous time limit is performed, the lack of an epidemic threshold in the spread of infections is not limited to metapopulations with a scale-free architecture, as it has been predicted from dynamical equations in which reaction and diffusion occur sequentially in time.

The analysis of the spread of infectious diseases on complex networks has become a central issue in modern epidemiology ͓1͔ and, indeed, it was one of the main motivations for the development of percolation theory ͓2͔.While the initial approach was focussed on local contact networks ͓3-6͔, i.e., social networks within single populations ͑cities, urban areas͒, a new approach has been recently introduced for dealing with the spread of diseases in ensembles of ͑local͒ populations with a complex spatial arrangement and connected by migration ͓7͔.Such sets of connected populations living in a patchy environment are called metapopulations in ecology, and their study began in 1967 with the theory of island biogeography ͓8͔.
In some recent models of epidemic spreading, the location of the patches in space is treated explicitly thanks to the increasing of computational power ͑see, for instance, ͓9͔͒.In ͓7,10,11͔, however, an alternative approach based on the formalism used in the statistical mechanics of complex networks is presented.Precisely, the topology of the spatial network of local populations ͑nodes͒ is mathematically encoded by means of the connectivity ͑degree͒ distribution p͑k͒, defined as the probability that a randomly chosen node has degree k.Moreover, each node contains two types of particles: A particles corresponding to susceptible individuals, and B particles that correspond to infected ones.Within each node, a transmission process ͑reaction͒ occurs between particles of different type, and migratory flows take place among nodes ͑diffusion͒ at constant rates.Therefore, although the detailed description of the spatial network is lost, the approach offers an elegant description of the epidemic spread in terms of densities of A particles and B particles in patches of degree k at time t, here denoted by A,k ͑t͒ and B,k ͑t͒ respectively.The reaction and diffusion ͑RD͒ processes modeling the spread of an infectious disease are considered as a two-step process in ͓7,11͔.First, inside each network node, the reaction takes place under the assumption of a homogenous mixing and conserving the total number of particles.In particular, the spread of the infection within a population follows the dynamics of a susceptible-infected-susceptible model which is described by the reactions corresponding to the recovering ͑at a rate ͒ and transmission ͑at a rate ␤͒ processes, respectively.Second, after the reaction, fixed fractions 0 ഛ D A ഛ 1 and 0 ഛ D B ഛ 1 of each type of particles ͑A and B, respectively͒ randomly migrate along the links departing from the node.This two-step formulation of the RD processes will play a crucial role in the following and reflects the fact that the progress of the infection at the metapopulation level is conceived as a discrete-time process.The key point is that, since the continuous-time limit of this two-step process is not well defined, different predictions can result under each time formulation.To see this, let us write the discrete-time equations governing the evolution of A,k and B,k from time t to time t + according to the previous two-step formulation.Additionally, following ͓7͔, two different transmission mechanisms will be considered.If the transmissibility of the infection does not depend on the local population size k = A,k + B,k , then the transmission rate is constant ͑␤ k = ␤ 0 ͒ regardless of the node degree ͑type-I transmission͒ whereas, if there is a saturation in the number of contacts per time interval, the transmission rate is ␤ k = ␤ 0 / k ͑type-II transmission͒, which is the usual assumption in epidemiology.Since in a discrete-time setting we must talk about probabilities instead of rates, the model equations read as where and ␤ k are, respectively, the probability of recovering and of transmission during the time interval ͑t , t + ͒, and P͑kЈ ͉ k͒ is the conditional probability that a node of de-* joan.saldana@udg.edugree k is connected to a node of degree kЈ.Finally, note that diffusion from a node of degree kЈ occurs just at the end of the time interval where a fixed fraction D i of i particles move into a randomly chosen neighbor location with probability 1 / kЈ.Hence, diffusion is not modeled as a stochastic process that can happen at any moment within a time interval of length .
A crude continuous-time approximation to the previous equations consists in taking ‫ץ‬ t i,k ͑t͒Ϸ i,k ͑t +1͒ − i,k ͑t͒.This approximation of the time derivative, widely used in the literature, leads to the RD equations presented in ͓7͔ ͑and its supplementary information͒ which have the same equilibria as Eqs.͑2͒ and ͑3͒.However, whereas it works for one-step processes, it fails when dealing with composite Markov processes as the one described by this two-step process for the spread of infectious diseases.A basic assumption when simulating continuous-time processes is that one chooses the length of the time interval sufficiently small to ensure that the probability of more than one event occurring to the same individual is negligible ͓12,13͔.That is, all of the processes are disjoint ͑mutually exclusive͒ events which, in terms of the model, is equivalent to say that a given particle either reacts ͑i.e., becomes infected͒ or diffuses to another node during a time interval .By contrast, terms like ͑1 − D A ͒ B,k occurring in the limit equations obtained using the previous approximation of the time derivative come from the product of probabilities of the events "react" ͑at a rate ͒ of a B particle and "not diffuse" of the resulting A particle.In other words, "react" and "diffuse" are not considered as mutually exclusive events.In particular, this implies that, for both kinds of particles, the limit of ͓ i,k ͑t + ͒ − i,k ͑t͔͒ / as approaches zero is not defined due to the existence of unbounded terms of the form i,k ͑t͒ / in the expression of the incremental quotient for each type of particle, and hence the equations obtained by approximating ‫ץ‬ t i,k ͑t͒Ϸ i,k ͑t +1͒ − i,k ͑t͒ do not correspond to the continuous-time limit of the discrete equations ͑2͒ and ͑3͒.
On the other hand, if reaction and diffusion processes take place simultaneously, one must consider the diffusion probability D instead of the fraction D of diffusing particles.In this case, the limit of the incremental quotients as → 0 becomes well-defined and, after replacing D with D in Eqs.͑2͒ and ͑3͒, these become

͑5͒
These equations recover the classical expression of RD equations as a sum of two independent contributions: Reaction and diffusion, plus the contribution of diffusing particles from the nearest neighbors of nodes of degree k.For the sake of brevity, from now on we will consider strictly positive diffusion rates ͑D A , D B Ͼ 0͒.
In uncorrelated networks, P͑kЈ ͉ k͒ = kЈp͑kЈ͒ / ͗k͘ is the degree distribution of nodes that we arrive at by following a randomly chosen link ͓14͔, and the previous equations become where ͗k͘ = ͚ k kp͑k͒ is the average network degree, and i = ͚ k p͑k͒ i,k is the average number of i particles per network node ͑i = A , B͒.Note that, after multiplying Eqs.͑6͒ and ͑7͒ by p͑k͒ and summing over all k, it follows that that is, the total density of particles ͑t͒ = A ͑t͒ + B ͑t͒ remains constant with t and equal to 0 , the initial average number of particles per node, in agreement with the conservation of the number of particles characterizing reactions ͑1͒.
The equilibria of Eqs.͑6͒ and ͑7͒ satisfy

B,k
On the other hand, there exists an endemic equilibrium whenever the following equality holds: In type-II transmission, condition ͑11͒ leads to the following expressions of the densities of i particles ͑i = A , B͒ in nodes of degree k, with k * = k 0 / ͗k͘ being the density of particles at equilibrium in nodes of degree k.Moreover, according to ͑12͒, the mean density of each type of particle in the network is A * = 0 / ␤ 0 and B * = ͑1− / ␤ 0 ͒ 0 .From these expressions it follows that the condition for the existence of an endemic equilibrium is Moreover, it is not difficult to see that the dominant eigenvalue 1 of the Jacobian matrix of Eqs.͑6͒ and ͑7͒ at the disease-free equilibrium ͑10͒ is 1 = max͑␤ 0 − ,0͒.Hence, the disease-free equilibrium is unstable as long as there exists an endemic equilibrium.
Under type-I transmission, a sufficient condition for the disease-free equilibrium ͑10͒ to be unstable is ͑see the Appendix for details͒ 0 ജ ͗k͘ where k max is the maximum degree in the metapopulation.Hence, for fixed reaction rates and ␤ 0 , a high enough density of particles and/or a large enough maximum degree in the network guarantee the instability of the disease-free equilibrium independently of the topological fluctuations of the network.In the limit of very large networks with bounded average degree ͗k͘, this sufficient condition implies the lack of an epidemic threshold for degree distributions with k max → ϱ as, for instance, exponential and scale-free degree distributions.On the other hand, it is clear that an endemic equilibrium cannot satisfy condition ͑11͒ since, in this type of transmission, ␤ k * = ␤ 0 which implies that the left-hand side ͑lhs͒ of ͑11͒ does not depend on k.Therefore, at an endemic equilibrium, both sides of Eqs.͑8͒ and ͑9͒ must be different from zero.This fact implies that * are linear in k as it is the case in type-II transmission.In the particular case D A = D B it follows, moreover, that k * = k ͗k͘ 0 .To illustrate these predictions under type-I transmission, let us consider a metapopulation with degree distribution p͑1͒ = 0.1, p͑2͒ = 0.3, p͑3͒ = 0.5, p͑4͒ = 0.1.Moreover, let us assume D A = D B =1, =2, ␤ 0 = 1.5, and three different initial densities: ͑i͒ A,k 0 = 0.52, B,k 0 = 0.55; ͑ii͒ A,k 0 = B,k 0 = 0.65; ͑iii͒ A,k 0 = B,k 0 = 3.Note that we are taking Ͼ ␤ 0 which does not allow for the existence of an endemic equilibrium in type-II transmission.In case ͑i͒, 0 = ͚ k=1 4 p͑k͒͑ A,k 0 + B,k 0 ͒ = 1.07 and the only stationary solution corresponds to a stable diseasefree equilibrium.In case ͑ii͒, 0 = 1.3 and it is equal to the right-hand side ͑rhs͒ of condition ͑13͒.For this value of 0 , the existence of a stable endemic equilibrium relatively far from the disease-free equilibrium ͑ B * Ϸ 0.15͒ clearly shows that condition ͑13͒ is sufficient but not necessary.In case ͑iii͒, 0 = 6 which largely satisfies the sufficient condition for the instability of the disease-free equilibrium.From Table I it follows that the densities A,k * are linear in k at the diseasefree equilibrium, whereas they are clearly sublinear in k at the endemic equilibrium.Consequently, at this equilibrium, B,k * are superlinear in k with A,k * + B,k * = k 0 / ͗k͘.In this Brief Report we have analyzed the suitability of the continuous-time approximation of difference equations modeling the time evolution of two processes, reaction and diffusion, in metapopulations which take place sequentially.The conclusion is that having sequential processes during each time interval leads to an undefined continuous approximation of the discrete-time equations.However, if reaction and diffusion processes are assumed to take place simultaneously, one obtains well-defined RD differential equations which, as expected, are consistent with the conservation of the number of particles during the progress of the epidemic disease in a metapopulation.In this sense, it is important to realize that Eqs.͑4͒ and ͑5͒ are not the continuous-time limit of Eqs.͑2͒ and ͑3͒ since the latter are derived under the assumption of dispersal of individuals at the end of each time interval, once reactions ͑1͒ have taken place.In particular, this means that the latter equations cannot be obtained from the former ones by discretizing time.
The difference in the predictions of both formulations of the diffusion process arises for the type-I ͑nonsaturated͒ transmission of infection in uncorrelated networks.On the one hand, in the limit of infinite metapopulation sizes, a lack of an epidemic threshold is predicted for scale-free architectures according to the dynamics of sequential RD processes governed by Eqs.͑2͒ and ͑3͒ with =1 ͑cf.͓7͔͒.On the other hand, topological fluctuations of the network architecture measured in terms of ͗k 2 ͘ turns out to be not essential in simultaneous RD processes described by Eqs.͑4͒ and ͑5͒.In this case, the lack of an epidemic threshold is also predicted in the limit of infinite network sizes for any degree distribution with bounded mean degree but with an unbounded maximum degree.
The explanation for this disagreement between both predictions is that continuous diffusion is more effective for the spread of infectious diseases in metapopulations than that TABLE I. Equilibrium densities in type-I transmission for D A = D B =1, =2, ␤ 0 = 1.5, ͗k͘ = 2.6, and different initial densities 0 of particles.In case ͑i͒, 0 = 1.07 and the disease-free equilibrium is the only stationary solution.In cases ͑ii͒ and ͑iii͒, 0 = 1.3 and 6, respectively, and an endemic equilibrium is the only stable stationary solution.

Equilibrium densities
Case occurring at discrete times.In terms of the asymptotic behavior of the model, although this higher effectiveness is not noticeable when there is saturation in the transmission of the infection, it becomes crucial when such a saturation is not present.
Since the simultaneity of the reaction and diffusion processes seems to be the case when migratory flows occur continuously in time among local populations, from a modeling point of view the present analysis and discussion about the predictions obtained under two different approaches is not only a technical issue but one that we think must be considered in the analysis of more refined infection models defined on heterogeneous metapopulations.
The author thanks J. L. Garcia-Domingo for pointing out the model in Ref. ͓7͔ and for discussions on complex networks.This work has been partially supported by Grant No. MTM2005-07660-C02-02 of the Spanish government and by the project FP6-2003-NEST-PATH-1 "Unifying Networks for Science and Society" of the Sixth European Framework Programme.

APPENDIX
The Jacobian matrix of the system ͑6͒ and ͑7͒ at the disease-free equilibrium given by ͑10͒ is of the form where A, B, C, and 0 are k max ϫ k max matrices with 0 being the zero matrix.Therefore, the characteristic polynomial of J DF factorizes as p J ͑͒ = p A ͑͒p B ͑͒.In type-I transmission, p A ͑͒ = ͑ + D A ͒ k max −1 and the roots k p B ͑͒ are simple and satisfy for k =1, ... ,k max − 1, while the larger one satisfies Finally, since the roots of p A ͑͒ are nonpositive, the sufficient condition ͑13͒ follows from imposing the nonnegativity of the rhs of ͑A1͒, which is clearly sufficient but not necessary for the positivity of k max .