On the critical behavior of the Susceptible-Infected-Recovered (SIR) model on a square lattice
OOn the critical behavior of the Susceptible-Infected-Recovered (SIR) model on asquare lattice
Tˆania Tom´e and Robert M. Ziff Instituto de F´ısica, Universidade de S˜ao PauloCaixa Postal 66138, 05315-970 S˜ao Paulo SP, Braziland Michigan Center for Theoretical Physicsand Department of Chemical Engineering, University of MichiganAnn Arbor, MI 48109-2136, USA
By means of numerical simulations and epidemic analysis, the transition point of the stochastic,asynchronous Susceptible-Infected-Recovered (SIR) model on a square lattice is found to be c =0 . c is the probability a chosen infected site spontaneously recovers rather thantries to infect one neighbor. This point corresponds to an infection/recovery rate of λ c = (1 − c ) /c = 4 . − c ) / (1 + 3 c ) = 0 . I. INTRODUCTION
The Susceptible-Infected-Recovered (SIR) model is afundamental system in epidemiological modeling [1–8]and has been studied extensively within the context ofnon-equilibrium phase transitions and critical phenom-ena (i.e., [9–22, 24]). The model was developed to de-scribe the propagation of an epidemic that occurs duringa period of time much smaller than the lifetime of indi-viduals of a given population. It is assumed that the pop-ulation can be divided into three categories: Susceptible(S), Infected (I) and Recovered (R) individuals. Suscep-tible individuals become infected at a given rate throughcontact with infected individuals. Infected individualsrecover with a given rate and become immune and recov-ered. The model is capable of showing a threshold of theepidemic spreading as one increases the infection rate.The SIR process has been studied using different ap-proaches and contexts. Originally it was defined in 1927by Kermack and McKendrick [1] as a deterministic pro-cess by means of a set of ordinary differential equations;they showed that epidemics disappear before all the sus-ceptible individuals contract the disease. Afterwards themodel was given a stochastic description by means ofbirth and death processes [2], and later Grassberger [9]introduced a cellular automaton implementation, that is,a synchronous-update Markovian process on a lattice. Inthis paper we consider a stochastic, asynchronous-updatelattice version of the SIR model [19], in which lattice sitesare updated one at a time. This model is a special caseof the predator-prey stochastic lattice-gas model intro-duced by Satulovky and Tom´e [25, 26], and also consid-ered by Antal et al. [12, 27]. It is also a special caseof the Susceptible-Infected-Removed-Susceptible (SIRS)stochastic lattice gas model [19]. For the synchronous orasynchronous versions of the SIR process, one observesthat, as the model parameters are varied, a phase transi- tion takes place. This is a continuous phase transition be-tween two distinct regimes: one in which the populationremains susceptible (non-spreading regime) and anotherin which the epidemic spreads over the lattice (spread-ing regime), where a significant portion of the populationbecomes infected and eventually immune. At the transi-tion the system becomes critical and corresponds to theepidemic threshold. Cardy and Grassberger [10] arguedthat the transition found in the SIR cellular automaton isof the percolation universality class. This has been con-firmed in various ways by numerous studies in stochasticmodels with synchronous as well as asynchronous update[11, 12, 15, 16, 18, 19].In recent years there has been a great deal of interest inthe SIR model on networks and other systems (i.e., [13–19, 21–24, 28–35]). In this paper we consider the SIRmodel on the square lattice, which has wide applicationsto many physical problems such as the spread of diseasein plants, and which has been studied only to a limitedextent [9, 14, 19, 36]. One question that has come up inthe context of mean field and network studies is the ef-fect of heterogeneity in the infectious period [29, 30, 32].When the period of infection is fixed, there is a directmapping of the model to bond percolation on the lattice,and thus the transition point can be determined exactly[9, 14, 36, 37]. However, when the infection period isheterogeneous, such as the exponentially distributed in-fection period inherent to the original SIR model, thereis no exact solution. In this paper, we carry out a carefulnumerical study of the transition point for the exponen-tially distributed case, and then investigate the correla-tions in the transmission of the infection to the nearestneighbors. We work out the correlations explicitly, andfind that at the critical point there is a higher probabilityto infect more rather than fewer of the neighbors, unlikethe case of bond percolation where the distribution ofinfected neighbors is simply binomial with p = 1 / a r X i v : . [ c ond - m a t . d i s - nn ] O c t II we define the asynchronous SIR model in terms of themaster equation. In section III we describe the Monte-Carlo algorithm that we use, and in section IV we carryout a numerical analysis of the cluster size distribution ofrecovered individuals to precisely determine the criticalpoint of the model. In section V we show in detail howthe asynchronous SIR model differs from percolation ona local scale. In section VI we discuss the relation of ourresults to some other population biology models, and insection VII we give our conclusions.
II. THE SIR LATTICE MODEL
To model the dynamics of an epidemic with immu-nization we consider a stochastic lattice-gas model withasynchronous dynamics. The lattice plays the role of thespatial region occupied by the individuals and the latticesites are the possible locations for the individuals. Eachsite can be occupied by just one individual that can beeither a susceptible, an infected, or a recovered individ-ual, called, respectively, an S site, I site and R site. Ateach time step a site is randomly chosen and the follow-ing rules are applied: (i) If the chosen site is in state Sor R it remains unchanged. (ii) If it is in state I then(a) with probability c the chosen site becomes R and (b)with the complementary probability b = 1 − c a neighbor-ing site is chosen at random. If the chosen neighboringsite is in state S it becomes I; otherwise it remains un-changed. Notice that a site in state R remains foreverin this state so that the allowed transitions are S → I → R.From this set of dynamic rules it follows that the stateof the system will change only when the chosen site isin state I, a feature that will be used to speed up thesimulation as explained in section III. This algorithm isequivalent to making S → I with probability bn I /
4, where n I is the number of nearest-neighbor I sites, and I → Rwith probability c .The system evolves in time according to a master equa-tion. To each site i of a two-dimensional lattice, we as-sociate a stochastic variable η i that takes the values 0,1, or 2 according to whether the site i is occupied byan R or an S or an I individual, respectively. A micro-scopic configuration of the entire system is denoted bythe stochastic vector η = ( η , . . . , η i , . . . , η N ) where N isthe total number of sites. Because the possible transi-tions are the cyclic ones (1 → → η i obtained by a cyclic permutationof the state of site i , that is, η i = ( η , . . . , η (cid:48) i , . . . , η N )where η (cid:48) i is 1, 2, or 0 according to whether η i is 0, 1,or 2, respectively. The master equation for the probabil-ity distribution P ( η, t ) associated with the microscopicconfiguration η at time t , is given by ddt P ( η, t ) = (cid:88) j (cid:88) i (cid:48) { w Bij ( i η ) P ( i η ) − w Bij ( η ) P ( η ) } + + (cid:88) j { w Cj ( j η ) P ( j η ) − w Cj ( η ) P ( η ) } , (1)where the summation on i extends over the nearest-neighbor sites of site j and j η denotes the state obtainedfrom η by an anticyclic permutation of state (0 → → w Bij is the transition rate associated withthe infection process, given by w Bij ( η ) = βz δ ( η i , δ ( η j , , (2)where z is the number of nearest-neighboring sites of site j (the lattice coordination number) and δ ( x, y ) is theKronecker delta; and w Cj is the transition rate for therecovery process, given by w Cj ( η ) = γδ ( η j , . (3)Two external parameters are associated to these pro-cesses: the infection rate β and the recovery rate γ .The probabilities b and c are related to these rates by b = β/ ( β + γ ) and c = γ/ ( β + γ ) so that b + c = 1 , (4)as it should.From the master equation (1) we can derive the timeevolution equations for the densities of recovered ρ , sus-ceptible ρ , and infected ρ . The connection of thepresent stochastic lattice model with the approach de-veloped by Kermack and McKendrick [1] can be revealedby using a simple mean-field approximation. Within thisapproximation [19] the following set of ordinary differen-tial equations for the densities can be derived ddt ρ = − βρ ρ , (5) ddt ρ = βρ ρ − γρ , (6) ddt ρ = γ ρ . (7)These equations are essentially the equations introducedby Kermack and McKendrick [1] in their deterministicapproach for the spreading of a disease with immuniza-tion.To analyze an epidemic, one can begin from an initialcondition at time t = 0 where all the individuals are sus-ceptible, with the exception of a very small number ofinfected individuals; that is, ρ = 1 − ρ ∗ , ρ = ρ ∗ (cid:28) ρ = 0. Using this initial condition and Eqs. (5-7),one finds that the system evolves in time and reaches twotypes of states: one where the epidemic spreads, that is, ρ (cid:54) = 0, ρ (cid:54) = 0, ρ = 0, when t → ∞ , which occursfor sufficiently large values of the infection probability b ,and another where the epidemic does not spread, that is, ρ = 1 when t → ∞ , which occurs for small values of b . As one varies the parameter b , there is a continuousphase transition at b = 1 /
2. In stochastic lattice modelsone expects similar behavior but a distinct value for thecritical parameters. In that case one starts with a lat-tice full of susceptible individuals with the exception of asingle infected individual, and studies the epidemics thatensue.
III. SIMULATION ALGORITHM
The asynchronous SIR model may be simulated by akinetic Monte-Carlo process by following procedure: • Pick an I site randomly from a list of all I sites. • Generate a random number X in (0 , X ≤ c then let I → R. • Otherwise (if
X > c ), pick one nearest neighbor ofthe I site randomly. If the neighbor is S, then makeit I and add to list of I sites. • Repeat as long as there are available I sites.When we remove an I site from the list in the first step,we swap the empty location with the site at the top ofthe list and decrease the list length by one in order tokeep the list compact. New I sites in the third step areadded to the top of the list.The Monte-Carlo time t is determined by increment-ing t by 1 /N I , where N I is the current number of I siteson the list, each time an I is picked from the list. Thisreflects an effective increase of one Monte-Carlo step ifevery site on the lattice were chosen once, on the aver-age. However, as explained below, we don’t actually keeptrack of time in our simulations, but instead monitor justthe size (number of R sites) of the clusters.This model is closely related to standard percolationgrowth, in which an active growth site spreads to an un-occupied nearest neighbor with the given bond probabil-ity. However, in the percolation case, as described below,the growth site becomes inactive after the four nearestneighbors are checked (for both asynchronous and syn-chronous updating). In the SIR model, the I site re-mains active an exponentially distributed length of time,and can thus attempt to infect neighboring sites multipletimes. That behavior leads to a kind of correlation in thespreading in the SIR model, as we shall see.In order to study properties of this model in detail, it isnecessary to know the transition point to high accuracy.Recent work [19] has shown that c = 0 . c to about 200 times theaccuracy, and indeed, in so doing, we find nearly completeagreement with the previous central value. IV. NUMERICAL RESULTS
To find the transition point, we consider the statisticsof the cluster size distribution, and determine c as thepoint where the distribution follows a power law. Ac-cording to standard percolation theory, the probability P s that a point belongs to a cluster containing s sites isgiven by sn s , where n s is the number of clusters of size s (per site) on the lattice. At criticality, n s ∼ s − τ , so that P s ∼ s − τ . Integrating from s to ∞ we find the proba-bility P ≥ s that a point belongs to a cluster of size greateror equal to s . At criticality, P ≥ s ∼ s − τ , and within thescaling region one expects [38] P ≥ s ∼ s − τ F ( εs σ ) ≈ s − τ ( A + Bεs σ ) , (8)where ε = c − c and for the last term we made a Taylor-series expansion of the scaling function F for small val-ues of its argument [39]. In 2-d percolation, we have τ = 187 /
91 and σ = 36 /
91. While there is no “ n s ”for dynamic percolation (because we don’t have a latticefully populated with clusters), we expect P ≥ s for dynamicpercolation and the SIR model at criticality to show thesame power-law behavior as static percolation.To calculate P ≥ s , we ran simulations with the latticecovered entirely by S sites, except for a single I site inthe center, using the algorithm above. The lattice sizewas 16384 × s of the cluster was charac-terized by the number of R sites. When s went beyonda cut-off of 2 = 2097152, the growth was stopped. Thecluster sizes were binned to make a histogram. Variousruns were made at each value of c as listed in the captionof Fig. 1. For a pair of values of c which bracketed theapparent transition point, we generated about 10 clus-ters, requiring several weeks of computer time each. Therandom number generator we used was R(9689) of [40].While the time t might perhaps be more of a naturalvariable to consider for dynamic percolation, we chose toconsider the survival probability as a function of s . Oneadvantage of using s instead of t is that the percolationexponents for P ≥ s are known exactly in two dimensions,while those for P ≥ t ∼ t (2 − τ ) D/d min = t − . ... (at crit-icality) are related to d min which is known only approx-imately: d min ≈ . D = 91 /
48 is thefractal dimension.) Also, it is somewhat easier to keeptrack of the dependence upon a discrete variable s ratherthan a real variable t . In Fig. 1 we plot s τ − P ≥ s vs. s σ for various c , using the 2-d percolation values of τ and σ .According to (8), this plot should be a linear function,with a slope proportional to c − c , and slope equal tozero for c = c . Indeed, this linear behavior is followedvery well, except for smaller s where (8) is not valid dueto finite-size effects. in Fig. 2, we plot the slopes of thelinear parts of the three central values as a function of c . The intercept where the slope is zero gives the criticalpoint: c = 0 . s τ - P ≥ s s σ FIG. 1. (Color online) Plot of s τ − P ≥ s vs. s σ for theasynchronous SIR model, using dynamic percolation values τ = 187 /
91 and σ = 36 /
91, for (top to bottom for large s ): c = 0 . . , . . where the error is based upon the statistical error of thedata. The linear behavior of the curves in Fig. 1 forlarger s shows that the critical behavior is consistent withpercolation scaling. In Fig. 3, we plot ln P ≥ s vs. ln s at c = 0 . − . − τ = − / ≈ . c ≈ . → Rwith probability c and S → I with probability n I (1 − c ) / → Rwith probability c (cid:48) = c/ (1 − c ) and S → I with probability n I /
4. We also confirmed by simulation that this proce-dure yields the correct critical point c (cid:48) = c / (1 − c ) ≈ . transmissibility ) our result implies by (12),prob(infection) = 0 . p (bond) c = 1 / p (site) c ≈ . V. RELATION TO PERCOLATION GROWTH
While the asynchronous SIR model is clearly in thedynamic percolation universality class, it is not strictlyidentical to percolation. In this section we show explicitlyhow the two differ locally by considering the probabili-ties of an infection spreading from a single I site. Thisdifference between the SIR model and bond percolation -8 0 8 0.176494 0.1765 0.176506 S l o p e c FIG. 2. Slope of the linear part of the three central curves ofFig. 1, measured for s σ ≥ . The intercept of the line onthe horizontal axis gives an estimate for c . was noted by Kuulasmaa [36] and by Grassberger [9], andmore recently discussed in general by Kenah and Robins[29]. Here we carry out a brief analysis directly related tothe computer algorithm we developed, and give explicitnumerical results for the infectivity or transmissibility atthe critical point.First, for comparison, we formulate epidemic bond per-colation growth in the SIR language. As in the algorithmin section III, we start with a system of all S, with one Isite at the center. The algorithm we follow for percola-tion is: • Pick an I site randomly from the list of I sites. Setthis I to R and remove it from the list. • Consider the I site’s four nearest-neighbor sites,and for each do the following: • If the neighbor is an S site, then generate a randomnumber X in (0 , • If X ≤ p then let the S become I, and add to list;otherwise, do nothing. • Repeat as long as there are available I sites.This is equivalent to bond percolation because eachpossible bond is considered only once and with a fixedprobability p . Note that the algorithm does not generateall the bonds of a cluster, as no internal bonds are consid-ered, but it finds all the “wetted sites” with the correctprobability. The bonds that are occupied form a mini-mum spanning tree on the cluster. A similar algorithmwas given by Grassberger [9].If a given I site has m nearest-neighbor S sites ( m =0 , . . . , P ( m ) k that exactly k ofthem become infected as a result of the central I sitefor percolation is just P ( m ) k = (cid:18) mk (cid:19) p k (1 − p ) m − k . (10)At the threshold p c = 1 /
2, for m = 4, for example, wehave P (4)0 = 1 / P (4)1 = 4 / P (4)2 = 6 / P (4)3 = 4 / P (4)4 = 1 /
16, with an average occupancy of 2.Next we calculate P ( m ) k for the SIR model. First, wederive the probability of spreading for “transmissivity”for the SIR model. Consider an I site surrounded by atleast one S site. We desire the probability that a givenone of those S sites will become infected by the I site.The I site will remain infectious for an exponentially dis-tributed number of trials. The probability p ( n ) that itremains infectious for n trials is simply p ( n ) = (1 − c ) n c n = 0 , , , . . . , (11)because c is the probability that it recovers in a giventrial. This implies that, on the average, an I site willbe considered (cid:80) ∞ n =0 np ( n ) = (1 − c ) /c times before itrecovers. y = -0.05524x - 0.20471 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0 5 10 15 l n P ≥ s ln s FIG. 3. Plot of ln P ≥ s vs. ln s for the SIR model, simulatedon a 16384 × s = 2097152.The equation represents a linear fit of the points for s > x represents ln s and y represents ln P ≥ s . For each trial that occurs on the given I site, we checkone of its four neighbors randomly to see if any is S, notremembering if we have already checked that neighborbefore. The probability that a given nearest neighborS site is chosen in at least one of the n trials is givenby 1 − (3 / n . Multiplying by p ( n ) from Eq. (11) andsumming, we find the net probability a given nearest-neighbor S site becomes infected:Prob(infected) = ∞ (cid:88) n =0 (1 − c ) n c (cid:20) − (cid:18) (cid:19) n (cid:21) = 1 − c c (12)If we set this probability equal to 1 / c = 1 /
5, which is abovethe observed value c = 0 . / not give the correct threshold. This is because that bond is notoccupied independently, but is correlated with its neigh-bors, and so the effective bond probability (probability ofspreading to a neighbor) is not 1/2. Putting c = 0 . . − p (0) = 0 . p (site) c , or c = 1 − p (site) c = 0 . p (site) c is the site perco-lation threshold on the square lattice [42–44], and haveused p (0) from Eq. (11). This also implies Prob(infected)= 0 . z .Then an extension of equation (12) to arbitrary z givesProb(infected) = b/ [ z − ( z − b ]. Setting this equal to1 / ( z − z −
1) nearest-neighborS sites and the critical point is when, on the average, onenearest-neighbor site is infected, we find [23] b = z z − . (13)Thus, for z = 3, the transition is at c = 1 /
4, and for z = 4, the transition is at c = 1 /
3. The latter is muchgreater than the critical c we found for the square lattice.This is expected, because for the Bethe lattice the netbond probability should be less than for a regular latticewith the same z .Now, we turn to the probabilities P ( m ) k for the SIRmodel. From a recursion relation analysis, and using thegenerating function method, we find the coefficients foreach m = 1 , , ,
4. The derivation is given in the Ap-pendix. The resulting P ( m ) k calculated from Eqs. (28-29)in the Appendix at c = 0 . k for a given m .The probabilities add to 1, and the mean number of in-fected neighbors is b/ (4 − b ) m = 0 . m , consistentwith Eq. (12).For m = 4, the results give the outgoing probabil-ities in the construction of Kuulasmaa [36], in whicha quenched double-directed lattice showing all possiblespreading probabilities is considered. In this construc-tion, one does not consider whether the neighbor is S ornot, but draws outgoing directed bonds (with the properprobabilities) to all neighboring sites. Then an epidemicis identified by following all directed paths emanatingfrom a given seed.Thus, we see in general that in the SIR model at crit-icality there is an enhancement in the number of neigh- TABLE I. P ( m ) k for the SIR model at c = 0 . p c = 1 /
2, which is also at a criticalpoint but with much different P (4) k . m \ k p c = 1 /
2. The red (dark) dot marks the place where thecluster growth began.FIG. 5. (Color online) Example of a completed SIR cluster of51034 recovered sites, generated at the threshold c = 0 . bors that are infected, compared to the case for dynamicpercolation. Because of the correlations in the infectedneighbors in the SIR model, it is necessary for the netneighbor infection probability to be higher than the ran-dom percolation value of 1/2.We note that for c = 1 /
5, the value obtained by settingthe probability of spreading (12) equal to the bond perco-lation value 1 / P ( m ) k is simply a uniformly distribution P ( m ) k = 1 / ( m + 1). This behavior is in sharp contrast to(10) for p = 1 / p .For the z = 4 Bethe lattice at its critical point c = 1 / P (4) k are 5 /
15, 4 /
15, 3 /
15, 2 /
15 and 1 / k = 0, 1, 2, 3, and 4, respectively—trending in theopposite direction as for the square lattice.In Figs. 4 and 5 we show actual pictures of critical clus-ters containing about 50,000 sites each from the percola-tion and SIR models, respectively. On a local scale theSIR clusters appear slightly denser but otherwise there isno apparent difference between the two critical clusters. VI. RELATION TO OTHER POPULATIONBIOLOGY MODELS
The SIR model can be considered as a particular case ofother population biology models such as the Susceptible-Infected-Recovered-Susceptible (SIRS) model [17, 19]and the Predator-Prey model [18]. The SIRS model de-scribes an epidemic process without permanent immu-nization and is defined by the following three processes:S → I, I → R, and R → S. The first two processes, S → I andI → R, are the same as those for the SIR model, as de-scribed in section 2, and occur with rates β and γ , re-spectively. The third process R → S is spontaneous andoccurs at rate α . The Predator-Prey model, in the epi-demic language, is similar to the SIRS model except thatthe third process, R → S occurs with rate αn S /
4, where n S is the number of nearest-neighbor S sites. In the epidemiclanguage the following correspondence is used: prey as S,predator as I, and an empty site as R.The SIRS and the Predator-Prey models exhibit non-equilibrium phase transitions between an absorbing sus-ceptible phase and active phase where the individuals arecontinuously being infected. Their critical behavior be-longs to the universality class of directed percolation for α (cid:54) = 0. When α = 0 one recovers [18, 19], from thesetwo models the SIR model, which belongs to the univer-sality class of dynamic isotropic percolation, as we haveconfirmed.In the opposite regime, namely, when α is large enoughcompared to β and γ , both the SIRS and the Predator-Prey models map [19, 45] into the contact process (CP)[46] with a creation rate λ = β/γ . The CP can be identi-fied as the Susceptible-Infected-Susceptible (SIS) modelof epidemiological modeling [8]. The SIS model describesthe dynamics of infection with no immunity and has twoprocesses: S → I with a rate βn I /
4, and I → S with rate γ .Changing the time scale, these two processes can be de-scribed by the infection rate λ = β/γ and recovery rateequal to 1. A site occupied by a particle in the CP modelcorresponds to an I site in the SIS model and an emptysite to an S site.The SIR and the SIS models correspond then to twoextremal behaviors of the SIRS model with respect toimmunity. Let us consider the periods of time spent byan individual in the R state. These periods of time aredistributed around a mean value T , which is proportionalto the inverse of the rate α . In the SIRS model T is finitewhich means that individuals have a partial immuniza-tion. The SIR model can be understood as the SIRSmodel in which T → ∞ , meaning that an R individ-ual has a lifelong immunity. The SIS model, on the otherhand, can be regarded as the SIRS model in which T → b = β/ ( α + β + γ ), c = γ/ ( α + β + γ ) and a = α/ ( α + β + γ ), so that b , c and a areinterpreted as the probabilities of infection, recovery andre-infection, respectively. In this formulation there arejust two independent parameters, as a + b + c = 1.For small values of β and γ compared to α , the crit-ical line of the SIRS model can be obtained by meansof the mapping into the SIS or CP [19, 45]. In two di-mensions the critical value for the creation rate for theCP is λ c = 1 . λ = β/γ and b/c = β/γ , it follows that the critical line is given by b/c = λ c = 1 . a →
1. If one extrapolates thisline to a = 0, b + c = 1, which corresponds to the SIRmodel (that is, suppressing the parameter a ), one gets(1 − c ) /c = λ c or c = 1 / ( λ c + 1) ≈ . → S rather than tries to infect one neigh-bor. This value should be compared to the critical valuefor the SIR process c ≈ . b = 1 − c ≈ . b = 1 − c ≈ . VII. CONCLUSIONS
We have provided further evidence that the asyn-chronous SIR model is in the universality class of stan-dard percolation through the behavior of the cluster sizedistribution. We showed that local correlations differfrom those of standard percolation. By extensive numer-ical simulation of the cluster size distribution, we haveshown that the transition in the SIR model defined on the square lattice occurs at c = 0 . c compares with c = 1 / . . c is useful for otherstudies of this critical state, such as studying the scalingof the average cluster size with L for finite systems [50]. Note added in revision . While this paper was underreview, a paper appeared online [51] in which heterogene-ity of the transmissibility in a square-lattice SIR modelwas also considered. In that paper, the heterogeneityis formed by having classes of individuals each with dif-ferent fixed transmissibilities ψ i corresponding to havingdifferent fixed infections times. In terms of the corre-lations we studied, each of these classes of individualsproduce a binomial distribution for P ( m ) k given by Eq.(10) with p replaced by ψ i , and by taking a set of classesit is possible for example to reproduce the exponential in-fectivity distribution exactly. Our critical transmissibil-ity 0 . σ = 0 . VIII. ACKNOWLEDGMENTS
This work was supported by the Brazilian agencyCNPq and INCT of Complex Fluids (CNPq andFAPESP) (T.T.), and supported in part by the U. S.National Science Foundation Grant No. DMS-0553487(R.M.Z.).
IX. APPENDIX
Consider that there are m nearest-neighbor S sites toa given I site, and four nearest neighbors total. Then wedefine: • P ( m ) n,k = the probability that exactly k distinct Ssites are visited after n trials on the I site, k =0 , . . . m , for n = 0 , , , . . . ∞ .with P ( m )0 , = 1; P ( m )0 ,k = 0, k = 1 , . . . , m . Given p ( n )defined in (11), the net probability that k sites of type Sare visited is given by P ( m ) k = ∞ (cid:88) n = k p ( n ) P ( m ) n,k = ∞ (cid:88) n = k c (1 − c ) n P ( m ) n,k , (14)We can write recursion relations to find the P ( m ) n,k . Forexample, for m = 1 we have P (1) n, = 34 P (1) n − , (15) P (1) n, = 14 P (1) n − , + P (1) n − , . (16)The first expression states that the probability that the Ssite was not visited at the n -th trial is 3 / n -th trial, while if the S site wasvisited in some previous trial, then it will surely still havebeen visited in this trial.To solve the recursion relations, we use the generatingfunctions, defined for general m by G ( m ) k ( x ) = ∞ (cid:88) n = k P ( m ) n,k x n . (17)A straightforward application to (15-16) yields G (1)0 ( x ) = 11 − (3 / x (18) G (1)1 ( x ) = x − x )(1 − (3 / x ) . (19)According to (14), we have simply P ( m ) k = cG ( m ) k (1 − c ) , (20)and we thus find P (1)0 = 4 c c (21) P (1)1 = 1 − c c , (22)which satisfies P (1)0 + P (1)1 = 1. Likewise, for m = 2, we have P (2) n, = 24 P (2) n − , (23) P (2) n, = 24 P (2) n − , + 34 P (2) n − , (24) P (2) n, = 14 P (2) n − , + P (2) n − , (25)and similarly one can write the recursions for m = 3and 4. In fact one can summarize them for all m by thegeneral formulas P ( m ) n, = 4 − m P ( m ) n − , (26) P ( m ) n,k = m + 1 − k P ( m ) n − ,k − + 4 − m + k P ( m ) n − ,k (27)for k = 1 , . . . , m , and derive a general recursion relationfor the averaged quantities for all m = 1 , . . . , P ( m )0 = 4 cm + (4 − m ) c (28) P ( m ) k = ( m − k + 1)(1 − c ) m − k + (4 + k − m ) c P ( m ) k − k = 1 , ..., m (29)Using these formulas, we obtain the numerical resultsgiven in Table I. For m = 4 (the last row in that table),the explicit formulas are: P (4)0 = c (30) P (4)1 = 4(1 − c ) c c (31) P (4)2 = 12(1 − c ) c (3 + c )(2 + 2 c ) (32) P (4)3 = 24(1 − c ) c (3 + c )(2 + 2 c )(1 + 3 c ) (33) P (4)4 = 6(1 − c ) (3 + c )(2 + 2 c )(1 + 3 c ) (34)Note that it is possible to write explicit formulas forall the P ( m ) n,k , such as P (3) n, = 4 n − n ) + 3(2 n ) − n , (35)where the coefficients in some cases are related to Stirlingnumbers of the second kind. Also, we can generalize (28-29) to all coordination numbers z simply by replacing all4’s in those equations by z ’s. Then one can show that P ( m ) m + P ( m ) m − + . . . + P ( m ) m − n +1 = n (1 − c ) zc P ( m ) m − n , (36)for n = 1 , . . . , m , and taking n = m , one can verifydirectly that (cid:80) mk =0 P ( m ) k = 1. [1] W. O. Kermack and A. G. McKendrick, Proc. Royal Soc.London A , 700 (1927). [2] N. T. J. Bailey, Biometrika , 177 (1953); The Mathe-matical Theory of Epidemics (Hafner, New York, 1957). [3] K. Dietz, J. Roy. Stat. Soc. A , 505 (1967).[4] D. Mollison, J. Roy. Stat. Soc. B , 283 (1977).[5] E. Renshaw, Modelling Biological Populations in Spaceand Time (Cambridge University Press, Cambridge,1991).[6] A. Hastings,
Population Biology: Concepts and Models (Springer, New York, 1996).[7] J. D. Murray,
Mathematical Biology (Springer, New York,2003).[8] M. J. Keeling and P. Rohani,
Modeling Infectious Dis-eases in Human and Animals (Princeton UniversityPress, Princeton, 2008).[9] P. Grassberger, Math. Biosc. , 157 (1983).[10] J. Cardy and P. Grassberger, J. Phys. A , L267 (1985).[11] H. Hinrichsen, Adv. Phys. , 815 (2000).[12] T. Antal, M. Droz, A. Lipowski and G. ´Odor, Phys. Rev.E , 036118 (2001).[13] M. E. J. Newman, Phys. Rev. E , 016128 (2002).[14] L. M. Sander, C. P. Warren, and I. M. Sokolov, PhysicaA , 1 (2003).[15] S. M. Dammer and H. Hinrichsen, Phys. Rev. E ,016114 (2003).[16] S. M. Dammer and H. Hinrichsen, J. Stat. Mech.: Theor.Exp. P07011 (2004).[17] J. Joo and J. L. Lebowitz, Phys. Rev. E , 036114(2004).[18] E. Arashiro and T. Tom´e, J. Phys. A , 887 (2007).[19] D. R. de Souza and T. Tom´e, Physica A , 1142 (2010).[20] M. Henkel, H. Hinrichsen, and S. L¨ubeck, NonequilibriumPhase Transitions , vol. 1 (Springer, Berlin, 2008).[21] G. M. Sch¨utz, M. Brandau, and S. Trimper, Phys. Rev.E , 061132 (2008).[22] A. J. Black, A. J. McKane, A. Nunes, and A. Parisi,Phys. Rev. E , 021922 (2009).[23] T. Tom´e and M. de Oliveira, to be published.[24] V. R. V. Assis and M. Copelli, Phys. Rev. E , 061105(2009).[25] J. E. Satulovsky and T. Tom´e, Phys. Rev. E , 5073(1994).[26] J. E. Satulovsky and T. Tom´e, J. Math. Biol. , 344(1997). [27] T. Antal and M. Droz, Phys. Rev. E , 056119 (2001).[28] M. B. Hastings, Phys. Rev. Lett. , 148701 (2006).[29] E. Kenah and J. M. Robins, Phys. Rev. E , 036113(2007).[30] J. C. Miller, Phys. Rev. E , 010101(R) (2007).[31] D. Alonso, A. J. McKane, and M. Pascual, J. Roy. Soc.Interface, , 575 (2007).[32] E. Volz, J. Math. Biol. , 293 (2008).[33] E. Ben-Naim and P. L. Krapivsky, Phys. Rev. E ,050901(R) (2004).[34] D. A. Kessler and N. M. Shnerb, Phys. Rev. E ,010901(R) (2007).[35] M. A. Serrano and M. Bogu˜n´a, Phys. Rev. Lett. ,088701 (2006).[36] K. Kuulasmaa, J. Appl. Probab. , 745 (1982).[37] K. Kuulasmaa and S. Zachary, J. Appl. Probab. , 911(1984).[38] D. Stauffer and A. Aharony, Introduction to PercolationTheory , Taylor and Francis, 1994.[39] C. D. Lorenz and R. M. Ziff, Phys. Rev. E , 230 (1998).[40] R. M. Ziff, Computers Phys. , 385 (1998).[41] P. Grassberger, J. Phys. A , 6233 (1999).[42] M. E. J. Newman and R. M. Ziff, Phys. Rev. E , 016706(2001).[43] M. J. Lee, Phys. Rev. E , 031131 (2008).[44] X. Feng, Y. Deng, and H. W. J. Bl¨ote, Phys. Rev. E ,031136 (2008).[45] K. C. de Carvalho and T. Tom´e, Int. J. Mod. Phys. C , 1647 (2006).[46] T. E. Harris, Ann. Probab. , 969 (1974).[47] R. Dickman, Phys. Rev. E , R2441 (1999).[48] M. M. S. Sabag and M. J. de Oliveira, Phys. Rev. E ,036115 (2002).[49] T. Vojta, A. Farquhar and J. Mast, Phys. Rev. E79