Impact of density-dependent migration flows on epidemic outbreaks in heterogeneous metapopulations
aa r X i v : . [ q - b i o . P E ] J u l Impact of density-dependent migration flowson epidemic outbreaks in heterogeneous metapopulations
J. Ripoll, A. Aviny´o, M. Pellicer, and J. Salda˜na
Departament d’Inform`atica, Matem`atica Aplicada i Estad´ıstica,Universitat de Girona, 17071 Girona, Catalunya (Spain) (Dated: May 10, 2019)We investigate the role of migration patterns on the spread of epidemics in complex networks. Weenhance the SIS-diffusion model on metapopulations to a nonlinear diffusion. Specifically, individ-uals move randomly over the network but at a rate depending on the population of the departurepatch. In the absence of epidemics, the migration-driven equilibrium is described by quantifyingthe total number of individuals living in heavily/lightly populated areas. Our analytical approachreveals that strengthening the migration from populous areas contains the infection at the earlystage of the epidemic. Moreover, depending on the exponent of the nonlinear diffusion rate, epi-demic outbreaks do not always occur in the most populated areas as one might expect.PACS numbers: 89.75.Fb 87.23.Cc
I. INTRODUCTION
Density dependence is a crucial factor in dispersalmovements of individuals in many species, be they an-imal or human. In turn, these movements are a key fac-tor in the fate of epidemic outbreaks occurring withinspatially distributed populations [1].Migration or dispersal, as the movement of individualsfrom one place to another, are affected by the presenceof other conspecific members in two different ways. Apositive density dependence reflects the fact that highpopulation densities increase competition effects amongtheir individuals and induce emigration from heavily pop-ulated locations. On the contrary, a negative density de-pendence corresponds to conspecific attraction, i.e., thetendency for individuals of a species to settle near con-specifics, and it has been claimed to be one of the reasonsfor leaving low-density areas [2] and for population ag-gregation patterns [3]. In general, there is empirical evi-dence that diffusion rates can both increase and decreasewith population density. For instance, in some insectgroups, crowded conditions lead to the appearance of agreater fraction of long-winged adults [4] and, in aphids,the appearance of winged individuals within populationsof wingless adults [5]. On the other hand, negative den-sity dependence has been reported, for instance, in gulls,voles, and deers [2]. In humans, both positive and neg-ative dependences have been considered for prehistoricand historical human population dispersals [6].Traditionally, populations are assumed to be con-tinuously distributed over space and, hence, nonlinearreaction-diffusion equations have been used to study theeffect of the previous density dependencies on the spa-tial population dynamics [3, 7]. In this context, lineardependence of the diffusion coefficient on the populationdensity has been proposed by several authors for animaldispersal [8, 9], whereas a power law with positive expo-nent has been used to describe the relationship betweenthe insect dispersal rate and the population density [7].Alternatively, the increasing fragmentation of habitats of many species, as well as the fact that, at a large scale,human travel can be described by flows among a setof discrete locations, make metapopulation models veryuseful for studying the dynamics of spatially subdividedpopulations and, also, for the analysis of spreading pro-cesses on top of these populations. Under this modelingapproach, the spatial distribution of the whole popula-tion is described by a network of patches inhabited bylocal populations, as cities or habitats in a patchy land-scape, with migratory flows connecting them.The nature of these flows is, in fact, a key ingredient ofmetapopulation models. On the one hand, flows can bedue to uniformly random migration, which assumes thatany neighbouring patch of a “source” patch is reachedwith the same probability. In this case, the migrationprobability only depends on features of the source patchas, for instance, its population density or its degree orconnectivity. On the other hand, migratory flows be-tween patches can depend on features of both source anddestination patches. Such a non-uniformly random mi-gration is, in fact, assumed by the so-called gravity modelof movements of people and goods traditionally used insocial sciences [10, 11], and by the radiation model intro-duced in [12] to overcome some of its limitations. In anepidemic context, migration towards other patches canbe even more related to the healthy conditions in thedestination patches than to the conditions in the originpatch. In [13], for instance, a constant diffusion rate isassumed for each class of individuals, but the migrationprobability between two patches depends on the suscep-tible population in the destination patch.The goal of this paper is to deal with positive andnegative density-dependent diffusion processes extendingthe equations for continuous-time epidemic dynamics onmetapopulations derived in [14], which were based on theformalism of complex networks introduced in [15, 16]. Inparticular, we study the impact of density-dependent dif-fusion coefficients on the population distribution amongheavily populated and lightly populated areas, how thesedistributions affect the epidemic growth, and the contri-bution of each local population to the epidemic spread-ing.
II. THE MODEL
The spatial arrangement of patches (areas) is describedin terms of the connectivity (degree) distribution p ( k )and the conditional probability P ( k ′ | k ) for a patch (node)of degree k to be connected to a patch of degree k ′ . More-over, within a patch of degree k , any individual has thesame probability of leaving it through any of its k links,namely, 1 /k , which implies that the strength of the con-nections is independent of the travel distance betweenpatches.According to these assumptions and denoting by ρ S,k ( t ) and ρ I,t ( t ) the average number of suscepti-ble and infectious individuals in a patch of degree k ( k = 1 , . . . , k max ), respectively, the following equationsdescribe the epidemic dynamics over a metapopulation ρ ′ S,k ( t ) = (cid:16) µ − βc ( ρ k ) ρ S,k ρ k (cid:17) ρ I,k + δ ( ρ k − ρ S,k ) − D S ( ρ k ) ρ S,k + k X k ′ D S ( ρ k ′ ) P ( k ′ | k ) ρ S,k ′ k ′ ,ρ ′ I,k ( t ) = (cid:16) βc ( ρ k ) ρ S,k ρ k − µ (cid:17) ρ I,k − δρ I,k − D I ( ρ k ) ρ I,k + k X k ′ D I ( ρ k ′ ) P ( k ′ | k ) ρ I,k ′ k ′ , (1)where β denotes the infection transmission probabilitythrough an infectious contact, µ is the recovery rate, andequal birth and death rates δ are assumed. The aver-age population size in a patch of degree k is ρ k ( t ) = ρ S,k ( t ) + ρ I,k ( t ), and the contact rate c ( ρ k ) is a non-decreasing density-dependent function, generalizing thetwo cases considered in [14, 17]: c ( ρ ) = ρ (fully-mixedpopulation) and c ( ρ ) = 1 (limited homogeneous mixing).Finally, the density-dependent diffusion rates of suscep-tible and infectious individuals are denoted by D S ( ρ k )and D I ( ρ k ), respectively. Note that it is natural to as-sume that the total outflow of individuals of each typein a patch, D S ( ρ ) ρ and D I ( ρ ) ρ , are zero when ρ = 0(see Discussion for a comment on the implications of theviolation of this hypothesis).The first term in each equation of (1) corresponds tothe infection and recovery processes. The second term isthe neutral demographic turnover. The last one corre-sponds to the migration/diffusion process which can besplit into negative and positive terms. The former countthe number of individuals leaving a patch of degree k ,whereas the latter are the sum of the flows of individualsarriving at this patch from patches of degree k ′ , providedthat patches of degree k are connected to patches of de-gree k ′ , i.e., for those k ′ such that P ( k ′ | k ) > S ( t ) = N P k p ( k ) ρ S,k ( t ) and I ( t ) = N P k p ( k ) ρ I,k ( t ), respectively, where N is the numberof nodes of the network. From Eq. (1) and assuming theconsistency condition kP ( k ′ | k ) p ( k ) = k ′ P ( k | k ′ ) p ( k ′ ), itfollows that the total number of individuals is conservedin the metapopulation, i.e., S ( t ) + I ( t ) = N ρ , with ρ being the average number of individuals per patch.Denoting by ~ρ S and ~ρ I the vectors of the populationdistribution of each type in the metapopulation, system(1) can be written as ~ρ ′ S = diag (cid:0) ˜ µ − β c ( ρ k ) ρ S,k ρ k (cid:1) ~ρ I +( C − Id ) diag( D S ( ρ k )) ~ρ S ,~ρ ′ I = diag (cid:0) β c ( ρ k ) ρ S,k ρ k − ˜ µ (cid:1) ~ρ I +( C − Id ) diag( D I ( ρ k )) ~ρ I , (2)where ˜ µ = µ + δ , C denotes the connectivity matrixwith C kk ′ = kP ( k ′ | k ) /k ′ , and diag( · ) denotes a diago-nal matrix. C is non-negative and assumed to be irre-ducible. For the case of uncorrelated networks, P ( k ′ | k ) = k ′ p ( k ′ ) / h k i and, hence, C kk ′ = kp ( k ′ ) / h k i , where thebrackets stand for the mean value. III. DIFFUSION WITHOUT EPIDEMICS
In the absence of epidemic ( ~ρ I = ~ ρ k , weconsider the equation ~ρ ′ S = ( C − Id ) diag( D S ( ρ k )) ~ρ S andstudy the disease-free equilibrium ~ρ ∗ S . As C is irreducible, v k = k , k = 1 , . . . , k max , is the only positive eigenvectorof C− Id associated to the dominant eigenvalue λ = 0. So, ρ ∗ k satisfies D S ( ρ ∗ k ) ρ ∗ k = M k , with M being a suitableconstant. Assuming that the total outflow of individualsper patch F ( ρ ) := D S ( ρ ) ρ is continuous, strictly increas-ing, and satisfies the natural hypothesis F (0) = 0, theexistence and uniqueness of a disease-free equilibrium ρ ∗ k = F − ( M k ) (3)is guaranteed and, hence, M is computed from the nor-malizing condition P k p ( k ) ρ ∗ k = ρ . It is worth notingthat ρ ∗ k does not depend on the conditional probability P ( k ′ | k ), so it is independent of the degree-degree corre-lations, and, moreover, that it increases with k since thetotal outflow per patch F ( ρ ) is an increasing function.To deal with both positive and negative density depen-dencies, we will assume the following typical nonlinearform for the diffusion rate D S (see [3] and the referencestherein) D S ( ρ ) = D S (cid:18) ρρ (cid:19) α , α > − , (4)with D S > − ), which guarantees theprevious hypotheses on the total outflow per patch F ( ρ ). Regarding to the dimensionless parameter α , for −1 0 1 2 3 4 52030405060708090100 α % o f popu l a t i on i n H P a r ea s .
19% 36 . very large scale−free networks FIG. 1. (Color online) Percentage of population in heavilypopulated (HP) patches, (( ω − /ω ) ω − × ω =( γ − α ) >
1, as a function of the dimensionless migrationexponent α in (4), and degree distribution p ( k ) ∼ k − γ with k min = 3 and γ = 3. Solid dots correspond to the profiles for α = − / α = 0 when the ratioof the individuals living in heavily/lightly populated patchesis one, and α = 1 (lightly populated scenario). − < α < ρ < ρ ). On the other hand, an exponent α > ρ > ρ ), and it hasbeen used for modelling animal dispersal in [7] (Section11.3). Finally, α = 0 recovers the constant diffusion rateconsidered in [14, 17].Taking (4), the expression (3) for the disease-free equi-librium reads ρ ∗ k = k / (1+ α ) h k / (1+ α ) i ρ . (5)Note that this power-law distribution arises indepen-dently of the network topology. Moreover, other typesof profiles are possible by just changing the form in (4).It is worth mentioning that, in terms of our notation,the same profile as (5) is found in [18] by assuming aconstant overall diffusion rate D from a patch, but dif-fusion rates between pairs of patches depending on thedegrees of the involved nodes. Therefore, more than onemechanism of dispersal can explain power-law profiles of ρ ∗ k (see Discussion).Beyond the richer variety of population profiles foundso far, it is interesting to know to what extent these sta-tionary profiles correspond to metapopulations whose in-dividuals live mostly in heavily populated patches or, onthe contrary, in lightly populated ones. Precisely, we say that a patch of connectivity k is heavily populated (HP)if its population is greater than the average number ofindividuals per patch in the metapopulation, i.e., when ρ k > ρ . Otherwise, a patch is said to be lightly popu-lated ( ρ k ≤ ρ ).When the disease-free equilibrium is computed from(3), HP patches turn out to be those with connectivitylarger than k ∗ := D S ( ρ ) ρ /M , with M the normalizingconstant for the profile. To illustrate heavily/lightly pop-ulated scenarios, we assume the form (4) for the diffusionrate D S and a scale-free network with degree distribution p ( k ) = γ − k min ( k/k min ) − γ , k ≥ k min , γ >
2. From (5), HPpatches are those with connectivity k > k ∗ = (cid:18) ωω − (cid:19) α k min where ω = ( γ − α ), and the total population in HPpatches, in the limit of very large networks, is N Z ∞ k ∗ p ( k ) ρ ∗ k dk = N ρ (cid:18) ω − ω (cid:19) ω − , (6)with ω > ω = 2, the percentages ofindividuals living in both types of patches are equal to50%. That would correspond, e.g., to the case of constantdiffusion rate α = 0 and γ = 3 as in [17]. However, othermigration patterns allow us to describe a wider spectrumof population distributions, specifically, those with a per-centage of population living in HP patches higher ( ω < ω >
2) than 50%, see Fig. 1 for the case γ = 3and varying α . For any exponent γ > γ > α can be used as a tuning parame-ter to shape the profile to a metapopulation with a givenpercentage of population in HP patches.However, the percentage of population living in HPpatches, namely, (( ω − /ω ) ω − × . , ω , with the limit values 100 e − %and 100% obtained when ω → ∞ and ω →
1, respec-tively. Notice that it is possible to obtain 100% of thetotal population in highly populated patches because thecomputation in (6) assumes infinite network sizes.
IV. EARLY STAGE OF THE EPIDEMIC
Now, let us focus on the effect of the migration patternson the onset of the epidemic. The initial epidemic growthrate is given by the dominant eigenvalue (i.e. spectralbound) λ of the Jacobian matrix of (2) at the disease-free equilibrium. This matrix can be written in blocksas: −0.5 0 0.5 1 1.5 2 2.5 3050100150200250 migration exponent α i n i t i a l g r o w t h r a t e o f t he ep i de m i c nonlimited contact rateestimationactual value −0.5 0 0.5 1 1.5 2 2.5 3678910111213141516 migration exponent α i n i t i a l g r o w t h r a t e o f t he ep i de m i c saturated contact rateestimationactual value β c − ˜ µ FIG. 2. (Color online) Comparison of the initial epidemic growth rate (time − ) for two contact patterns: c ( ρ ) = 100 ρ/ρ (left) and c ( ρ ) = 1000 ρ/ ( ρ + ρ ) (right). Growth rate computed as its actual value (dots) given by the dominant eigenvalue λ of the Jacobian matrix ( k min = 3 , k max = 220), and its estimation (solid line) given by left-hand side of (8). Resultsobtained for uncorrelated scale-free networks ( γ = 3 , N = 5000 nodes) and parameters: β = 0 .
015 (dimensionless); ˜ µ = 1, D S ( ρ ) = D I ( ρ ) = 0 . ρ/ρ ) α (units of time − ), and ρ = 100 individuals. (cid:18) ( C − Id ) · diag ( F ′ ( ρ ∗ k )) − diag ( βc ( ρ ∗ k ) − ˜ µ )0 ( C − Id ) · diag ( D I ( ρ ∗ k )) + diag ( βc ( ρ ∗ k ) − ˜ µ ) (cid:19) , (7)where the total outflow per patch F ( ρ ) = D S ( ρ ) ρ is as-sumed to be differentiable. An exponential initial growthof the infectious population occurs when the disease-freeequilibrium becomes unstable, i.e., when λ >
0. Theblock structure of the Jacobian matrix and the fact thatthe dominant eigenvalue of the first block is 0, imply that λ is given by the dominant eigenvalue of the fourth blockwhen the latter is positive.For a constant contact rate, c ( ρ ) = c , the dominanteigenvalue is equal to λ = max { , βc − ˜ µ } and it turnsout to be independent of the migration pattern. On thecontrary, when the contact rate is density dependent, asufficient condition for λ > k { βc ( ρ ∗ k ) − ˜ µ − (1 − P ( k | k )) D I ( ρ ∗ k ) } > λ since it ne-glects off-diagonal terms of the connectivity matrix.To have a global description of the role of the migra-tion, we have computed the initial growth of the epidemic λ as the largest eigenvalue of the fourth block of the Ja-cobian matrix, for different migration exponents α andwe have checked the goodness of its estimation given bythe left-hand side of (8). Fig. 2 shows the fit for two con-tact patterns, namely, nonlimited and saturated numberof contacts per unit of time. The first is the standardmass action c ( ρ ) = c ρ and in the second we consider a saturated contact rate of the form c ( ρ ) = c ρ/ ( ρ + ρ ), c >
0, which comes from the assumption that the timean individual is available for contacts is limited [19]. Inboth cases shown in Fig. 2, the initial growth rate of theepidemic decreases as the migration exponent increases.Recalling the definition of the basic reproduction num-ber of a local population as the number of secondaryinfections produced by a typical individual in a whollysusceptible population, and computed as infection prob-ability × contact rate × average infectious period , we canrewrite (8) asmax k βc ( ρ ∗ k )˜ µ + (1 − P ( k | k )) D I ( ρ ∗ k ) > . (9)So, for a given k , the ratio in (9) is a lower bound of thebasic reproduction number of the populations living inpatches with connectivity k , since the immigration frompatches with connectivities k ′ = k is neglected. Note thatthe factor 1 − P ( k | k ) = P k ′ = k P ( k ′ | k ) in the denominatoraccounts for those individuals who emigrate to patcheswith other connectivities. Therefore, local epidemic out-breaks will undoubtedly take place in those patches withconnectivity k where the ratio in (9) is larger than one.In contrast to the case of constant diffusion [17], fordensity-dependent contact rates c ( ρ ) and for D I ( ρ ) = τ D S ( ρ ), τ >
0, with D S of the form (4), numerical com-putations show that the maximum in (9) is not necessar-ily attained at populations with the largest connectivity. k k k α = 1 α = 3 α = 1 . patch connectivity k α
50 100 150 2001.41.61.822.22.42.62.83 0.50.60.70.80.911.11.21.3
FIG. 3. (Color online) Dimensionless ratio in (9) for uncorrelated scale-free networks ( γ = 3 and k min = 3) and diffusion D S ( ρ ) = D I ( ρ ) = D · ( ρ/ρ ) α . Left panel, from left to right: migration exponent α = 1, 1.7, and 3. Right panel: filled contourplot of the ratio in (9) showing that the maximum value decreases as α increases for the considered range of the parametervalues. Parameters: β = 0 .
015 (dimensionless); ˜ µ = 1, c ( ρ ) = 100 ρ/ρ , D = 0 . − ), and ρ = 100 individuals. In some cases, this maximum is attained at the minimumdegree k min , and, in others, at an intermediate value of k (see left panel in Fig. 3 for an illustration of these cases).In addition, the right panel of Fig. 3 is a plot of theratio in (9) as a function of both the degree k and themigration exponent α , and shows that the higher the mi-gration exponent α , the lower the maximum value (9),for the considered range of parameters values. V. DISCUSSION
Using the modelling framework introduced in [14], adensity-dependent diffusion rate D ( ρ ) has been consid-ered to model positive and negative density dependenciesof individuals’ migratory movements within a metapop-ulation [2]. In order to analyse the effect of these de-pendencies on the dispersal process, D ( ρ ) has been as-sumed to be a power-law function of ρ with an exponent α > −
1. This lower bound for the allowed values of α guarantees both that the flow of individuals leavinga patch, D ( ρ ) ρ , strictly increases with its population ρ ,and that D ( ρ ) ρ → ρ → D ( ρ ) ρ , thedisease-free equilibrium is determined by the diffusionprocess (which does not need to be necessarily describedin terms of a power-law diffusion rate) and it is computedfrom the eigenvector associated to the largest eigenvalueof the connectivity matrix. With this respect, it is in-teresting to observe that if these hypotheses are violatedas by assuming that the number of individuals travelingout of a patch, D ( ρ ) ρ , is independent from its populationsize ρ , the diffusion process does not determine the pro- file of the disease-free equilibrium. Such an assumption,indeed, leads to a diffusion rate D inversely proportionalto the patch population ρ , and to an equilibrium whichis uniquely determined by the initial distribution ρ k (0)(see Sect. 3.2 in [18]).On the other hand, when D ( ρ ) ∼ ρ α ( α > − ρ ∗ k ∼ k / (1+ α ) independentlyof the network topology. So, a power-law profile for ρ ∗ k has been obtained under the only assumption of nonlin-ear diffusion. Remarkably, the same type of profile hasbeen obtained for uncorrelated networks in [18] under adifferent physical principle, namely, assuming that themigratory flow ω kk ′ between an origin patch of degree k and a destination patch of degree k ′ is w kk ′ = w ( kk ′ ) θ , θ >
0. From these flows, the diffusion rate between nodesof degree k and k ′ is taken to be equal to D w kk ′ /T k with T k = k P k ′ P ( k ′ | k ) w kk ′ , which ensures a constant diffu-sion rate D from any patch in the metapopulation [18].Which of these physical assumptions is mainly at workin a random dispersal process depends on the species weare interested in. In general, it is expectable that manyinsects and other animal species experience only localenvironmental conditions and, hence, density dependen-cies at the origin location would be the main factor inrandom dispersal within a metapopulation. On the con-trary, humans migrate with knowledge of the political,economic and environmental conditions in the potentialarrival patches. Under this circumstance, it is suitableto assume that the migratory flow between two patchesdepends on features of the involved patches, and, in fact,it is what is assumed in many human mobility models[11, 12]. So, distinct physical mechanisms can be behindthe observed scaling laws in the topological features ofmobility networks [10].A way to quantify the effect of different density depen-dencies is by computing the percentage of individuals liv-ing in heavily populated patches (i.e., above the averagepopulation per patch). This aggregate measure of thepopulation over the whole metapopulation is a comple-mentary description to that given by ρ ∗ k . In particular,for infinite scale-free networks and assuming D ( ρ ) ∼ ρ α ,this percentage turns out to be bounded from below by36 . α is large enough, i.e., when there is a high propensity toemigrate from heavily populated patches.From an epidemiological point of view, the main nov-elty of our results is to show that migration patterns de-termine where epidemic outbreaks can take place withhigher probability. In particular, for large enough mi-gration exponents, the onset of an epidemic may not betriggered by infectious individuals living in large size pop-ulations, as it could be expected from the fact that ρ ∗ k increases with the connectivity of the sites. This hap-pens when individuals, also the infectious ones, have ahigh propensity to move from heavily populated areas.Such a higher number of infectious individuals arrivingat patches with lower population sizes makes more likelythe occurrence of an outbreak in these patches. But, atthe same time, their lower sizes lead to lower values of thebasic reproduction number because the contact rate c ( ρ )increases with ρ (see Fig. 3). Therefore, in contrast towhat is usually claimed about the role of the long-rangemixing in exacerbating epidemics [20], we have seen thatstrengthening the migration from heavily populated ar-eas can help to contain an epidemic by causing its ap-pearance in less populated areas, for which the value ofthe basic reproduction number is significantly smaller.This work is partially supported by grants MINECOMTM2011-27739-C04-03 and MTM2014-52402-C3-3-P(Spain). The authors are members of the Catalan re-search group 2014 SGR 1083. [1] S. Meloni, N. Perra, A. Arenas, S. Gomez, Y. Moreno,and A. Vespignani, Sci. Rep. (2011).[2] E. Matthysen, Ecography , 403 (2005).[3] V. M´endez, D. Campos, I. Pagonabarraga, and S. Fedo-tov, J. Theor. Biol. , 113 (2012).[4] R. Kisimoto, Nature , 641 (1956).[5] R. Harrison, Annual Review of Ecology and Systematics , 95 (1980).[6] J. Steele, Human Biology , 121 (2009).[7] J. D. Murray, Mathematical Biology , 3rd ed., Interdisci-plinary Applied Mathematics (Springer, 2007).[8] W. Gurney and R. Nisbet, J. Theor. Biol. , 441 (1975).[9] N. Shigesada, J. Math. Biol. , 85 (1980).[10] D. Brockmann, “Human mobility and spatial disease dy-namics,” in Reviews of Nonlinear Dynamics and Com-plexity (Wiley-VCH Verlag GmbH & Co. KGaA, 2010)Chap. 1, pp. 1–24.[11] A. G. Wilson,
Complex Spatial Systems: The ModellingFoundations of Urban and Regional Analysis. (Prentice Hall, Harlow, UK, 2000) pp. viii+174.[12] F. Simini, M. C. Gonzalez, A. Maritan, and A.-L.Barabasi, Nature , 96 (2012).[13] B. Wang, L. Cao, H. Suzuki, and K. Aihara, Sci. Rep. (2012).[14] J. Salda˜na, Phys. Rev. E , 01290 (2008).[15] V. Colizza, R. Pastor-Satorras, and A. Vespignani, Nat.Phys. , 276 (2007).[16] V. Colizza and A. Vespignani, Phys. Rev. Lett. ,148701 (2007).[17] D. Juher, J. Ripoll, and J. Salda˜na, Phys. Rev. E ,041920 (2009).[18] V. Colizza and A. Vespignani, J. Theor. Biol. , 450(2008).[19] H. R. Thieme, Mathematics in population biology (Princeton University Press, Princeton, NJ, 2003) pp.xx+543.[20] J. Epstein, J. Parker, D. Cummings, and R. Hammond,PLoS ONE3