Analysis and Monte Carlo simulations of a model for the spread of infectious diseases in heterogeneous metapopulations

We present a study of the continuous-time equations governing the dynamics of a susceptible-infectedsusceptible model on heterogeneous metapopulations. These equations have been recently proposed as an alternative formulation for the spread of infectious diseases in metapopulations in a continuous-time framework. Individual-based Monte Carlo simulations of epidemic spread in uncorrelated networks are also performed revealing a good agreement with analytical predictions under the assumption of simultaneous transmission or recovery and migration processes.


I. INTRODUCTION
The role of the contact structure in the spread of infectious diseases within populations has long been recognized ͓1-5͔.While the early susceptible-infected-removed ͑SIR͒ and susceptible-infected-susceptible ͑SIS͒ epidemic models were based on the assumption of well-mixed populations, recent models encompass complex contact networks describing the population mixing.Furthermore, several theoretical works have also highlighted the importance of the "contact structure" among populations at a geographical or metapopulation level, i.e., when considering ensembles of local populations with complex spatial arrangements and connected by migration ͑see for instance ͓6-12͔͒.
Recently, a new approach for the study of epidemic dynamics in heterogeneous metapopulations based on the formalism used in the analysis of complex networks has been introduced ͓13-15͔.Under this approach, the structure of the spatial network of patches ͑nodes͒ is encapsulated by means of the connectivity ͑degree͒ distribution p͑k͒ defined as the probability that a randomly chosen patch has connectivity k.Moreover, each patch contains two types of individuals ͑par-ticles͒: S ͑susceptible individuals͒ and I ͑infected individu-als͒.Within each patch, transmission or recovery ͑reaction processes͒ occur between individuals of different type.After reaction, migratory flows take place at once among patches ͑diffusion process͒ at constant rates.
In contrast, in ͓16͔ reaction and diffusion processes are considered to take place simultaneously, which turns out to be the correct assumption for a suitable continuous-time formulation of metapopulation models for the spread of infectious diseases.As in ͓13,15͔, the description of the epidemic spread is made in terms of average number ͑density͒ of susceptible and infected individuals in patches of connectivity k at time t, here denoted by S,k ͑t͒ and I,k ͑t͒, respectively.Reaction takes place under the assumption of a homogenous mixing in each local population.In particular, the spread of the infection within a population follows the dynamics of a SIS model which is described by the reactions

͑1͒
corresponding to the recovery and transmission processes, respectively.Here is the recovery rate and ␤ is the transmission rate across an infective contact.At the same time, migration ͑diffusion͒ of individuals S and I occurs with constant coefficients D S and D I , respectively.Each migrant ͑diffusing particle͒ randomly chooses one of the links departing from the patch.
In the present paper we are concerned with the analysis of the spread of infectious diseases in the modeling framework introduced in ͓16͔.We mainly focus on the case of nonlimited (or density-dependent) transmission of the disease.The latter means that the contact rate c ͑i.e., the per capita rate at which new contacts are made within a population͒ increases with the population size k ͑t͒ = S,k ͑t͒ + I,k ͑t͒ ͓17͔.Such a nonlimited transmission seems to be appropriate for the spread of viruses in computer networks and for the progress of animal diseases in well-mixed populations where the contact rate increases as more individuals are crowded into a given area ͓18,17͔.In other situations, for instance when there is a low number of individuals with which a given individual might interact or when the population density ͑numbers per unit area͒ remains constant as the number of susceptible and infected individual changes, the frequency of contacts per unit of time is independent of the population size.This corresponds to a limited (or frequency-dependent) transmission ͓17͔.
At this point, we must mention that similar equations to those presented in ͓16͔ ͑and analyzed here͒ were also presented in ͓19͔ under a more general modeling framework given by the so-called bosonic reaction-diffusion processes.The study in that paper, however, was focused on the analysis of one species reaction-diffusion processes, in which a single class of particles diffuses and reacts in the system.
Finally, in this paper we have also performed Monte Carlo simulations of the epidemic dynamics on uncorrelated metapopulations described above.The simulations confirm the analytical predictions of the model and, furthermore, show the suitability of the assumption of simultaneity of processes in order to avoid the effects of diffusion appearing under an updating rule in which reaction and diffusion steps alternate sequentially ͓19͔.

A. General equations
The rate at which new infected individuals are produced in a population, sometimes called transmission term ͓͑18͔͒, is equal to the product of the transmission rate ␤ and the total number of infective contacts per unit of time, that is, ␤ ϫ c ϫ S,k ͑t͒ ϫ I,k ͑t͒ / k ͑t͒.The last term in the product gives the proportion of contacts of susceptible individuals which take place with infected ones.In nonlimited transmission, the contact rate c equals k ͑t͒ and the transmission term is ␤ ϫ S,k ͑t͒ ϫ I,k ͑t͒.When limited transmission is assumed, c = 1 and the transmission term reads ␤ S,k ͑t͒ I,k ͑t͒ / k ͑t͒.
According to the derivation in ͓16͔ of the continuous-time formulation for the progress of diseases on metapopulations and assuming nonlimited transmission, the equations governing the dynamics of the disease propagation are

͑3͒
where k denotes the degree of the patches where local populations live ͑k =1, ... ,k max ͒, and P͑kЈ ͉ k͒ is the conditional probability that a patch of degree k has a connection to a patch of degree kЈ.As in classical reaction-diffusion processes, Eqs.͑2͒ and ͑3͒ express the time variation of susceptible and infected individuals as the sum of two independent contributions: reaction and diffusion.In particular, the diffusion term includes the outflow of individuals ͑diffusing par-ticles͒ from patches of degree k and the inflow of migratory individuals from the nearest patches of degree kЈ.
Analogous equations follow for limited transmission replacing ␤ by ␤ / k ͑t͒ in the previous equations ͑see ͓16͔͒.For the sake of brevity, from now on we will consider strictly positive diffusion rates ͑D S , D I Ͼ 0͒.
After multiplying Eqs.͑2͒ and ͑3͒ by p͑k͒, summing over all k, and using the consistency condition between p͑k͒ and P͑kЈ ͉ k͒ given by ͑see ͓20͔͒ kP͑kЈ͉k͒p͑k͒ = kЈP͑k͉kЈ͒p͑kЈ͒ to change the order of the summation indices in the double sums, it is immediate to see that the total density of individuals ͑t͒ = S ͑t͒ + I ͑t͒ remains constant and equal to 0 , the initial average number of individuals per patch.Here j = ͚ k p͑k͒ j,k is the average number of susceptible ͑j = S͒ and infected ͑j = I͒ individuals per patch.
Finally, for networks with a connectivity pattern defined by a set of conditional probabilities P͑kЈ ͉ k͒, we define the elements of the connectivity matrix C as Note that these elements are the average number of individuals that patches of degree k receive from neighboring patches of degree kЈ assuming that one individual leaves each of these patches by choosing at random one of the kЈ connections ͓͑21͔͒.On the other hand, for those degrees k that are not present in the network, P͑kЈ ͉ k͒ =0∀ kЈ.Hereafter in the paper, when talking about degrees, we implicitly mean those degrees that are present in the network.Furthermore, the case with patches having all the same connectivity is excluded from our considerations because, under the present approach, the model equations reduce to those of a singlepatch SIS model.

B. Uncorrelated networks
In order to obtain analytical results about the epidemic metapopulation dynamics, we need to be precise about the form of P͑kЈ ͉ k͒.The easiest and usual assumption is to restrict ourselves to uncorrelated networks.In this case, we have that P͑kЈ ͉ k͒ = kЈp͑kЈ͒ / ͗k͘ which corresponds to the degree distribution of nodes ͑patches͒ that we arrive at by following a randomly chosen link ͓22͔.
After replacing the expression of P͑kЈ ͉ k͒ into Eqs.͑2͒ and ͑3͒, one obtains the following equations for the epidemic spread in metapopulations described by uncorrelated networks: where ͗k͘ = ͚ k kp͑k͒ is the average network degree.
In these networks, the elements of the connectivity matrix C are simply given by Clearly, C is a rank-one matrix and has the vector with components v k = k as eigenvector of eigenvalue 1.So, if there are n different degrees in the network then the eigenvalues of this matrix are = 0, with algebraic multiplicity n − 1 and = 1 which is a simple eigenvalue.This fact will be used in the stability analysis of the equilibria of the model.

III. DISEASE-FREE EQUILIBRIUM
The equilibria of the model Eqs.͑2͒ and ͑3͒ are the solutions S,k ‫ء‬ , I,k ‫ء‬ to the equations For the analysis of the spread of the infection, it is particularly relevant the disease-free equilibrium.By definition, this is obtained by taking I,k ‫ء‬ = 0 in the previous equations which leads to the fact that the vector of densities S,k ‫ء‬ of susceptible individuals in patches with degree k is, in both types of transmission, an eigenvector of the connectivity matrix C normalized to 0 with associated eigenvalue =1: Since ͚ k Ј P͑kЈ ͉ k͒ =1∀ k, it easily follows that, for any generic network, the disease-free equilibrium is given by Moreover, if C is irreducible and transmission is limited ͑c =1͒, then the disease-free equilibrium is unstable if and only if ␤ Ͼ ͑see ͓21͔͒.

Instability of the disease-free equilibrium in uncorrelated networks with nonlimited transmission
Linearizing Eqs.͑4͒ and ͑5͒ around the disease-free equilibrium ͓Eq.͑9͔͒, it follows that the Jacobian matrix of the system is a block matrix of the form where each block is an n ϫ n matrix with n being the number of degrees in the metapopulation, 0 is the null matrix, I d is the identity matrix, diag͑a k ͒ stands for a diagonal matrix whose kth element is a k , and C is the connectivity matrix given by Eq. ͑6͒.Notice that the upper blocks of the Jacobian matrix are computed differentiating in Eq. ͑4͒ with respect to the susceptible and infected individuals in patches of degree k, respectively, and the lower blocks are computed analogously differentiating in Eq. ͑5͒.The triangular structure of the Jacobian implies that its characteristic polynomial factorizes as p͑͒ = p 1 ͑͒p 2 ͑͒ with p 1 ͑͒ and p 2 ͑͒ being the characteristic polynomials of the main diagonal blocks.So, the spectrum of the Jacobian matrix is the union of the spectra of these blocks.
From the knowledge of the eigenvalues of C and the fact that the spectrum of the matrix D S ͑C − I d ͒ is the one of D S C shifted by −D S , it follows that p 1 ͑͒ = ͑ + D S ͒ n−1 and, hence, that the largest eigenvalue of D S ͑C − I d ͒ is always = 0. On the other hand, from analyzing the changes of sign of the characteristic polynomial of the second main diagonal block .. ,n − 1, while the largest root satisfies n Ͼ k max / ͗k͘␤ 0 − ͑ + D I ͒ ͑see Fig. 1 for a numerical example and ͓23͔ for a general interlacing theorem of eigenvalues for perturbations of a diagonal matrix by rank-one matrices͒.In summary, all the eigenvalues of the Jacobian matrix of Eqs.͑4͒ and ͑5͒ at the disease-free equilibrium are real and the largest one is max = max͕0, n ͖, with Therefore, a sufficient condition for this equilibrium to be unstable is given by 0 Ն ͗k͘ This condition says that, for fixed , D I , and ␤, a high enough density of individuals per patch and/or a large enough maximum connectivity in the metapopulation guarantee the instability of the disease-free equilibrium ͓16͔.In the limit of very large networks with bounded average degree ͗k͘, this sufficient condition implies the lack of an epidemic threshold for any degree distribution with k max → ϱ.
An intuitive interpretation of condition ͑10͒ can be obtained in terms of the basic reproductive number R 0 , i.e., the average number of infections produced by an infected individual in a wholly susceptible population ͓24͔.From this definition, it is clear that the infection can spread in a population only if R 0 Ͼ 1. Rewriting Eq. ͑10͒ as where S,k max ‫ء‬ = k max 0 / ͗k͘, it follows that the left-hand side ͑lhs͒ of this inequality would correspond to the basic reproductive number R 0 in a patch with connectivity k max where infected individuals recover at a rate and migrate at a rate FIG. 1. ͑Color online͒ Eigenvalues of the lower main diagonal block of the Jacobian matrix, labeled by the patch connectivity, for both types of equilibria in uncorrelated networks and nonlimited transmission.The upper half of the picture shows the eigenvalues of this matrix around the ͑unstable͒ disease-free equilibrium.These are all real and increase linearly ͑and becomes positive͒ with the degree k ͑see main text for an explanation͒.The lower half shows the eigenvalues of the same matrix around the ͑stable͒ endemic equilibrium.These are all negative and decrease with k.The dashed line corresponds to Eq. ͑21͒ as a function of k.We have taken a metapopulation with scale-free degree distribution p͑k͒ϳk −3 with ͗k͘ =6 and k min = 3.The parameter values are 0 = 16, ␤ = 0.1, = 20, and D I to patches with different connectivities.Therefore, this expression does not take into account those infected individuals that can reach the patch from other patches with the same connectivity k max where they have also been introduced ͑recall that S,k correspond to the average number of susceptible individuals over all patches with connectivity k͒.So, the actual value of R 0 in patches with k = k max is then greater than the value given by the lhs of Eq. ͑11͒.It is possible, then, to have R 0 Ͼ 1 with the lhs of Eq. ͑11͒ being less than 1.For instance, if all the patches have the same connectivity, at equilibrium it would follow that In conclusion, condition ͑10͒ guarantees the spread of an infectious disease in those patches with the highest connectivity which, in its turn, allows the disease to eventually reach all the patches in the metapopulation.

IV. ENDEMIC EQUILIBRIUM IN UNCORRELATED NETWORKS
For the case of limited transmission, the expression for an endemic equilibrium of the model is explicit and equals to

͑12͒
with k ‫ء‬ = k 0 / ͗k͘ being the density of individuals ͑of both types͒ at equilibrium in patches of connectivity k.Note that Eq. ͑12͒ makes sense whenever ␤ Ն and at ␤ = this equilibrium bifurcates from the disease-free one.The uniqueness of this endemic equilibrium is proved in ͓21͔ under the assumption of equal migration rates ͑D I = D S ͒.

A. Existence of the endemic equilibrium with nonlimited transmission
When transmission is not limited, the expression of the endemic equilibrium is not explicit and, hence, the conditions for its existence must be derived from the analysis of the equilibrium equations.
An endemic equilibrium of the model is given by a positive solution S,k ‫ء‬ , I,k ‫ء‬ of the system where we recall that j ‫ء‬ = ͚ k p͑k͒ j,k ‫ء‬ is the average number of individuals of type j͑j = S , I͒ per patch at equilibrium.We are going to show that for any high enough 0 Ͼ 0 there exists a unique positive solution of Eqs.In particular, if we assume that migration rates for both types of individuals are equal, i.e., D S = D I , then using the conservation of the total density of individuals ͑at equilibrium͒ ‫ء‬ = S ‫ء‬ + I ‫ء‬ = 0 we get that is, any local population size at equilibrium k ‫ء‬ is proportional to the connectivity of its patch.
On the other hand, using the definition of ¯k ‫ء‬ and Eq.͑15͒ one has that system Eqs.͑13͒ and ͑14͒ is equivalent to which can be explicitly solved in terms of the average number of infected individuals per patch I ‫ء‬ .Indeed, from Eq. ͑16͒ we get with D IS ª D I / D S .Next, putting this expression in Eq. ͑17͒ we arrive at the following quadratic equation for I,k ‫ء‬ , regarding I ‫ء‬ as a parameter: Since the independent term above is positive ͑ I ‫ء‬ Ͼ 0͒, there exists a unique positive solution which is given by Let us point out that, according to Eqs. ͑18͒ and ͑19͒, the prevalence of the infection varies across patches with different connectivity.To end up it suffices to close the loop combining the definition of ‫ء‬ and expression ͑19͒.So, for each solution 0 Ͻ I ‫ء‬ Ͻ 0 of the scalar nonlinear equation there exists an endemic equilibrium given by Eqs.͑18͒ and ͑19͒.Considering the smooth function F͑ I ‫ء‬ ͒ defined by the right-hand side ͑rhs͒ of Eq. ͑20͒ minus I ‫ء‬ , it follows that Eq. ͑20͒ has a positive solution if condition ͑10͒ holds.On the one hand, we have that F͑ 0 ͒ Ͻ 0. On the other hand, whenever condition ͑10͒ holds with inequality, F͑0͒ Ͼ 0. When condition ͑10͒ holds with equality, then F͑0͒ = 0 with FЈ͑0͒ Ͼ 0. So, in both cases there exists a value I ‫ء‬ ͑0, 0 ͒ such that F͑ I ‫ء‬ ͒ = 0, i.e., a positive solution of Eq. ͑20͒.Additionally, the uniqueness of this positive solution is guaranteed if F͑ I ‫ء‬ ͒ does not change its curvature, which occurs when D I Յ D S or 0 is high enough.
So, for any high enough density of individuals per patch 0 fulfilling Eq. ͑10͒ there exists a unique positive solution of Eq. ͑20͒ and, hence, a unique endemic equilibrium of Eqs.͑4͒ and ͑5͒.
Finally, in the limit of very large networks with bounded average degree ͗k͘ and taking I ‫ء‬ as the solution of Eq. ͑20͒, it follows: So the number of infected individuals at equilibrium grows asymptotically as the degree of the patch.From Eqs. ͑18͒ and ͑19͒, on the other hand, we have that the number of susceptible individuals at equilibrium in any patch does not grow proportional to the degree but it remains bounded.More precisely, we have that However, it is not difficult to see that S,k ‫ء‬ is also increasing with degree k as it is I,k ‫ء‬ .Therefore, we can conclude that the prevalence of the infection I,k ‫ء‬ / k ‫ء‬ tends to 1 as the connectivity of the patches tends to infinity.
Moreover, if we assume equal migration rates ͑D I = D S ͒ then these limits have simpler expressions, namely, Let us point out that, in contrast to what happens in the case of limited transmission, the results obtained depend on both the migration rates and the architecture of the network.

B. Stability of the endemic equilibrium
Central to our analysis and simulations is the stability of the endemic equilibrium as well as the instability of the disease-free one.In the previous sections we have seen that, for uncorrelated networks and nonlimited transmission, con-dition ͑10͒ assures the instability of the disease-free equilibrium and the existence of an endemic one.
Let us linearize the equations for the epidemic spread around an endemic equilibrium.As in the disease-free case, the Jacobian matrix of systems ͑2͒ and ͑3͒ can be written in four blocks of dimension n ϫ n now given by In uncorrelated networks, assuming equal migration rates ͑D I = D S ͒ and combining block rows and then block columns, the previous Jacobian matrix can be transformed into the following block matrix: Finally, from Eqs. ͑18͒ and ͑19͒ it follows that ␤͑ S,k

͑21͒
So, max is negative if and only if the characteristic polynomial changes sign between max k ͕␤͑ S,k ‫ء‬ − I,k ‫ء‬ ͖͒ − ͑ + D I ͒ Ͻ 0 and zero.Evaluating the polynomial at both values ͑see ͓23͔ for a general expression of the characteristic polynomial for this type of matrix perturbations͒, there is a change in sign whenever Now, rewriting Eq. ͑14͒ as ͗k͘ this expression by p͑k͒ and summing over all k, it follows that which implies that condition ͑22͒ is always fulfilled since I,k ‫ء‬ Ͼ 0 at an endemic equilibrium.See Fig. 1 for an illustra-tion of eigenvalues of the second main diagonal block in both types of equilibria.Therefore, all the eigenvalues of the Jacobian matrix, coming from both main diagonal blocks, are negative except = 0 which is related to the conservation of the total density of individuals 0 .
Summarizing, we have shown that, for each 0 Ͼ 0 and equal migration rates, the endemic equilibrium is asymptotically stable whenever it exists.Moreover, numerical simulations of the model equations show that this equilibrium is globally asymptotically stable and exists for a wide range of parameter values fulfilling Eq. ͑10͒, including the case D I Ͼ D S .

V. MONTE CARLO SIMULATIONS OF THE MODEL
In order to check the predictions of the mean-field reaction-diffusion Eqs.͑4͒ and ͑5͒, we have carried out Monte Carlo simulations of this system of interacting individuals distributed in randomly generated uncorrelated networks of a fixed size N.
In each discrete time step of a Monte Carlo simulation, the rates D S , D I , ␤, and have to be, respectively, replaced by probabilities D S ⌬t, D I ⌬t, ␤⌬t, and ⌬t, where ⌬t is a small enough time interval ͑which will be recalculated in each iteration, see below͒.At each iteration, and for each infected individual, the algorithm randomly chooses one of the following three mutually disjoint events: ͑i͒ either the individual recovers with probability ⌬t ͑ii͒ or migrates with probability D I ⌬t ͑iii͒ or nothing happens with probability 1 − ͑ + D I ͒⌬t.
Analogously, for each susceptible individual: ͑i͒ either the individual becomes infected with probability ͑t , ⌬t͒ or ͑ii͒ migrates with probability D S ⌬t or ͑iii͒ nothing happens with probability 1 − ͑t , ⌬t͒ − D S ⌬t.
The form of the probability depends on the transmission type: limited or nonlimited.In the limited transmission, one assumes that each susceptible individual interacts with one randomly selected individual within the local population.Therefore, ͑t,⌬t͒ = I,k ͑t͒ S,k ͑t͒ + I,k ͑t͒ ␤⌬t.͑23͒ On the other hand, in the nonlimited transmission it is assumed that each susceptible individual interacts with all the infected local population.Thus, the probability that a susceptible individual becomes infected is equal to the probability of acquiring the infection at least once from the infectious contacts made during ⌬t, that is, ͑t,⌬t͒ = 1 − ͑1 − ␤⌬t͒ I,k ͑t͒ .͑24͒ This scheme makes sense if and only if ␤⌬t Ͻ 1, ͑ + D I ͒⌬t Ͻ 1, and ͑t , ⌬t͒ + D S ⌬t Ͻ 1 hold.The inequality ͑t , ⌬t͒ + D S ⌬t Ͻ 1 leads to ⌬t Ͻ ͑␤ + D S ͒ −1 when Eq. ͑23͒ holds.When Eq. ͑24͒ holds instead of Eq. ͑23͒, the previous inequality yields to ⌬t Ͻ x, where x is the smallest positive solution of the equation D S x = ͑1−␤x͒ n , with n being the maximum number of infected individuals per patch in the nth iteration.Summarizing, the algorithm is well defined, independently of the transmission type if, at each iteration, we take ⌬t smaller than min ͭ 1

͑25͒
The simulations have been performed on two uncorrelated networks of N = 5000 patches: a scale-free network generated via a Barabasi-Albert growing mechanism ͑or homogeneous linear preferential attachment, see ͓25͔͒ and an exponential network generated via a random nonpreferential growing mechanism.In both cases, the average degree is ͗k͘ = 6 and the minimum degree is 3.The lack of correlations has been tested by checking that the assortativity coefficient ͑introduced by Newman in ͓26͔͒ is of the order of 10 −3 for both networks.The initial number of individuals in patches of degree k is either equal to k +10 ͑which amounts to 0 =16͒ or to k + 100 ͑which amounts to 0 = 106͒.A 50% of infected individuals in all patches has been taken as the initial condition for all simulations.Finally, a time step ⌬t =2ϫ 10 −4 has been used in all cases.This value of ⌬t fulfills condition ͑25͒ for all the generated networks and for the initial conditions used in the simulations.To see why, note that the probability will be maximized when all the individuals living in patches with the highest connectivity are infected and, moreover, that the total number of individuals per patch increases linearly with the degree at equilibrium.For instance, considering the maximum degree occurring in the scale-free generated networks ͑k max = 214͒, taking D I = D S =1, ␤ = 0.6, and 0 = 106, and replacing the rhs of the equation for x by its linear approximation around x =0 ͑which yields a quite conservative estimate of x͒, one obtains a value of ⌬t of about 4 ϫ 10 −4 for the non-limited transmission, which is small enough even for variations in the number of infected individuals much larger than those occurring in the simulations.
The mean prevalence in nodes of degree k is the statistics we use to compare the model predictions with the Monte Carlo simulations: where k i , n i ͑t͒, and n I,i ͑t͒ are, respectively, the degree, the number of individuals, and the number of infected individuals in node i at time t, and N k is the number of nodes of degree k in the network.For large t, the output of the simulations is consistent with the predictions of the continuous-time equations.Precisely, for each degree k, the time average ͑from t =50 to t = 100 with a time series of one hundred values͒ of P k ͑t͒ is very close to that of the stable endemic equilibrium predicted by the model and given by Eqs.͑18͒-͑20͒.In addition, there is a low dispersion of the data with respect to the time-averaged value of P k ͑t͒ for each degree ͑see Figs.2-4͒.

VI. CONCLUSIONS
In this paper we have characterized the endemic equilibrium of the continuous-time SIS model for the spread of diseases in complex metapopulations introduced in ͓16͔.In this framework and assuming uncorrelated networks and nonlimited transmission, we have shown that the prevalence of the infection in local populations is not constant across the metapopulation but increases with the patch connectivity.Only when transmission is limited, the mean prevalence P k ͑t͒ is constant and equals 1 − / ␤ for all k, as predicted from Eq. ͑12͒ ͑see Fig. 3͒.
It is also remarkable the role of the migration rate of infected individuals D I when transmission is nonlimited.Diminishing the value of D I causes a reduction in the prevalence in patches of low connectivity, while those patches with higher connectivity have an even greater level of preva-lence ͑see Figs. 2 and 4͒.This is the only case we have observed in which the infection prevalence changes nonuniformly across the metapopulation when varying the value of a parameter.
The influence of the metapopulation architecture on the progress of the disease with nonlimited transmission is reflected in the fact that patches with higher connectivities are inhabited by more individuals.Since there is a long-term persistence of infected individuals in the SIS model, having large populations with R 0 Ͼ 1 inhabiting highly connected patches guarantees an outflow of infected individuals from those patches to the rest of the patches and, hence, the instability of the disease-free equilibrium.In particular, a lack of an epidemic threshold is predicted in very large metapopulations with the maximum connectivity tending to infinity ͓16͔.
No other especial topological features of the degree distribution ͑except the finiteness of the mean degree͒ are needed ͑cf.͓13͔͒.
This feature of the epidemic spread on metapopulations is likely to be reinforced in assortative networks, i.e., networks with high degree correlations, as it is the case of the geographical graphs in ͓27͔ inferred from data of a real metapopulation.In these networks and with nonlimited transmission, the role of the patches with the highest connectivity as sources of infected individuals is strengthened because these tend to be coupled to other highly connected patches and, then, the persistence of the infection increases in all of them.
In addition, we have also seen that, when reaction and diffusion are considered as mutually exclusive events over a short enough time interval ⌬t, individual-based Monte Carlo simulations are in agreement with the theoretical predictions.In particular, this means that the problem of having different behaviors of j,k ͑j = S , I͒ as a function of k when observed  after reaction or after diffusion is no longer present under the present simulation approach ͑cf.͓19͔͒.In other words, when reaction and diffusion take place simultaneously, individualbased algorithms for Monte Carlo simulations as the one used in the present paper lead to the same results as alternative algorithms based on "reaction-based" schemes where reactions instead of individuals are chosen using Monte Carlo techniques.Examples of the latter are, for instance, the algorithm introduced in ͓28͔ for the stochastic time evolution of coupled chemical reactions or the more recent one introduced in ͓19,29͔ to carry out similar simulations taking also into account diffusion of particles on complex networks and regular lattices, respectively.
* S,k − ρ * I,k ) − (µ + DI ) ͑13͒ and ͑14͒ with I ‫ء‬ Ͼ 0. First of all notice the following property of any equilibrium S,k ‫ء‬ , I,k ‫ء‬ of the system.Let us define the quantity ¯k ‫ء‬ ª D S S,k ‫ء‬ + D I I,k ‫ء‬ , which is interpreted as the ͑average͒ emigration flow in patches of connectivity k.Now the sum of Eqs.͑13͒ and ͑14͒ yields that ¯k ‫ء‬ is proportional to the degree k: ¯k ‫ء‬ = k D S S ‫ء‬ + D I I ‫ء‬ ͗k͘ .͑15͒ which is easier to analyze and has the same spectrum.As in Sec.III for the case of the disease-free equilibrium, the eigenvalues are just computed from both main diagonal blocks and the computations are analogous.The first block D I ͑C − I d ͒, depending only on the connectivity matrix, has zero as the largest eigenvalue and all the rest are equal to −D I Ͻ 0. On the other hand, from analyzing the changes in sign of the characteristic polynomial of the second main diagonal blockD I ͑C − I d ͒ + diag͓␤͑ S,k ‫ء‬ − I,k ‫ء‬ ͒ − ͔, itfollows that this polynomial has n simple real roots satisfying Ͻmax k ͕␤͑ S,k ‫ء‬ − I,k ‫ء‬ ͖͒ − ͑ + D I ͒ for all roots except for the largest one max which satisfies max Ͼ max k ͕␤͑ S,k ‫ء‬ − I,k ‫ء‬ ͖͒ − ͑ + D I ͒.
FIG. 2. ͑Color online͒ Prevalence of the infection in nodes of degree k of an uncorrelated scale-free network with nonlimited transmission.The continuous line corresponds to the prevalence I,k ‫ء‬ / k ‫ء‬ predicted by the model whereas dots are the time average of the Monte Carlo simulations P k ͑t͒ with its standard error.The mean prevalence per patch mpp = ͚ p͑k͒ I,k ‫ء‬ / k ‫ء‬ is also shown.Parameter values: the average number of individuals per patch 0 = 16, 106 ͑top/bottom͒, the migration rate D I = 0.1, 10 ͑left/right͒, and ␤ = 0.1, = 5, and D S =1.
FIG. 3. ͑Color online͒ Prevalence of the infection in nodes of degree k of an uncorrelated scale-free network with limited transmission.See caption of Fig. 2 for details.Here the predicted prevalence is constant and equals 1 − / ␤ across the network.Parameter values: 0 = 106, D I = 0.1, and ␤ =5, =2, D S =1.