An infection process near criticality: Influence of the initial condition
IInfection process near criticality: Influence of the initial condition
P. L. Krapivsky Department of Physics, Boston University, Boston, Massachusetts 02215, USA
We investigate how the initial number of infected individuals affects the behavior of the criticalsusceptible-infected-recovered process. We analyze the outbreak size distribution, duration of theoutbreaks, and the role of fluctuations.
I. INTRODUCTION
Epidemics greatly affect Homo sapiens [1–3] and otheranimals. Epidemic modeling goes back to Bernoulli whomodeled the spread of smallpox [4, 5]. The methodsand concepts developed in this field have been applied tomodeling recent human epidemics like HIV, H1N1 SwineFlu, Zika Virus, and COVID-19, and also to modelingnon-biological processes like the spread of rumors or cul-tural fads, diffusion of innovations, computer viruses, etc.(see [6–9] and references therein).Epidemics often operate close to the critical regime.Subcritical epidemics are the most numerous, but theyquickly die out. Supercritical epidemics can kill a finitefraction of the population thereby destroying the environ-ment that the virus needs to thrive and multiply. Epi-demics close to the critical regime are optimal in a viewof never-ending competition between the hosts and theviruses. Non-biological factors play an increasingly im-portant role in the spread of epidemics in human society.On one side, modern human society is much more in-terconnected than ever — this facilitates the spread ofinfections and leaves no hope for their total eradication.On the other side, advances in technology and wealthhelp to devise and employ the containment measures thatsuppress supercritical epidemics.Mathematical models of epidemics are over-simplified.In physics, for instance, we know that electrons are iden-tical. In contrast, organisms are different, even twins aredifferent. This heterogeneity is rarely taken into account,and it is far from clear how to model it in a reasonableway. In physics, we usually rely on binary and symmet-ric interactions. Both these features are questionable inthe realm of epidemics. Other realistic features are alsomostly ignored. However, the populations where the epi-demics spread are usually very large and the lore fromstatistical physics tells us that in large systems qualita-tive behaviors can be predicted even if one greatly simpli-fies the model. This is especially true in critical regimes.The very concept of the critical regime comes fromepidemic modeling. This concept clearly emerges fromthe well-known susceptible-infected-recovered (SIR) pro-cess [10–14], a toy model that mimics the spread of in-fection. According to the rules of the SIR process, in-fected individuals recover (become immune or die) withequal rates and every infected individual transmits a dis-ease to every susceptible individual with the rate R /N ,where N is the population size. Thus on average, each in-fected individual spreads the infection to R individuals before recovery. Therefore the behavior of the SIR pro-cess greatly depends on whether the reproduction num-ber R is smaller or larger than the recovery rate whichwe set to unity. When R <
1, i.e., for subcritical SIRprocesses, outbreaks quickly end, namely just a few indi-viduals catch the disease. For supercritical SIR processes( R > R − . With com-plementary probability, rN + O ( √ N ) individuals catchthe disease before the outbreaks dies out; the fraction r = r ( R ) is implicitly determined by r + e − R r = 1 (1)Huge outbreaks killing finite fractions of the popula-tion continue to devastate animal species. They also usedto decimate human societies [1–3], e.g., the Black Deathkilled about 50% of the European population [3]. Pre-ventive and containment measures such as quarantine,improved hygiene, etc. suppress supercritical infectiousdiseases and often drive them to a critical situation. Thiscritical state is effectively self-organized. Indeed, sup-pressing the disease to the subcritical regime may be pos-sible but costly and psychologically difficult to maintainwhen the number of newly infected starts to decreaseexponentially. Therefore, if the outbreak is not quicklyand completely eradicated, the containment measures arerelaxed and the system may return to the supercriticalstage, the disease gets again out of control, so the con-tainment measures are tightened driving the system backto the subcritical state. It would be interesting to devisea self-organized process of the spread of infection withdynamics similar to the critical SIR process. In this pa-per, however, we merely consider the critical SIR processwith many initially infected individuals.The SIR processes are often treated using a determin-istic framework [10–14]. This framework can be appliedto the supercritical regime where it gives e.g. the sim-plest derivation of Eq. (1). Stochastic effects are un-avoidable, however, for the critical and subcritical SIRprocesses (see [15–18]). When the population of suscep-tible is finite, finite-size corrections become important,particularly for the critical SIR process [19–25].In this paper, we study the critical SIR process andhence we employ stochastic methods. We consider fi-nite populations. The population size is assumed to belarge, N (cid:29)
1. We focus on the situation when the ini-tial number of infected individuals is also large: k (cid:29) a r X i v : . [ q - b i o . P E ] D ec single infected individual. Our goal, however, is to under-stand the behavior of the SIR process that has begun atthe supercritical regime, exhibited an exponential growthregime, and was subsequently suppressed (by preventiveand containment measures) to the critical SIR. Ignoringthe earlier regime yields the critical SIR process with acertain large number k of initially infected individuals.The same critical SIR process with a large numberof initially infected individuals has been recently stud-ied, using large-scale simulations and scaling arguments,by Radicchi and Bianconi [26]. We rely on exact cal-culations and asymptotic analyses. Our analytical andasymptotic predictions qualitatively agree with simula-tion results [26]. Some quantitive discrepancies suggesttrying slightly different scaling fits that are a little sim-pler than the fits used in [26]. The chief reason for subtlebehaviors are algebraic tails; the average size and dura-tion of outbreaks are especially sensitive to these tails.The chief outcome of the epidemic is the outbreak size n . The full description of this random quantity is pro-vided by the probability A n ( k, N ) that starting with k infected individuals in the population of size N , exactly n individuals catch the infection before the epidemic stops.We examine this probability distribution. In particular,we study the behavior of the average and variance E k ( N ) = (cid:104) n (cid:105) , V k ( N ) = (cid:104) n (cid:105) − (cid:104) n (cid:105) (2)We argue that in the N → ∞ limit, these quantitiesdepend on a single variable κ = k/N / . More precisely,they exhibit the following scaling behaviors E k ( N ) = N / E ( κ ) , V k ( N ) = N / V ( κ ) (3)We shall show that the scaled distributions have simpleextremal behaviors E ( κ ) = (cid:40) C κ when κ (cid:28) √ κ when κ (cid:29) V ( κ ) = (cid:40) C κ when κ (cid:28) (cid:112) /κ when κ (cid:29) k of initiallyinfected individuals and derive several exact results inthe infinite-population limit.The outline of this paper is as follows. In Sec. II, weconsider the infinite-population limit, present exact re-sults for the outbreak size distribution, and show thatexact results approach a simple scaling form in the mostinteresting situation when k (cid:29)
1. We then analyze the critical SIR process in a population with N (cid:29) II. INFINITE-POPULATION LIMIT:OUTBREAK SIZE DISTRIBUTION
In the infinite-population limit, the SIR process re-duces to the branching process [27–31]. Branching pro-cesses involve duplication and death. Symbolically, P + PP (cid:79) (cid:79) (cid:47) (cid:47) ∅ For the critical branching process, the rates of duplicationand death are equal.Branching processes have numerous applications, e.g.,they mimic cell division and death; for adult organisms,the critical branching process is appropriate as the num-ber of cells remains (approximately) constant.We begin with a classical setting when one individualwas initially infected. Let A n be the probability thatexactly n individuals catch the infection before the epi-demic is over. With probability , the initially infectedindividual joins the population of recovered before infect-ing anyone else, so A = . Further, A = A since atthe first step a new individual must get infected, and thenboth must recover without spreading infection. Proceed-ing along these lines we arrive at the recurrence A n = 12 (cid:88) i + j = n A i A j + 12 δ n, (5)reflecting that the first infection event creates two in-dependent infection processes [28]. A solution to (5) isfound by introducing the generating function A ( z ) = (cid:88) n ≥ A n z n (6)converting the recurrence (5) into a quadratic equation2 A ( z ) = [ A ( z )] + z (7)whose solution reads A ( z ) = 1 − √ − z (8)Expanding A ( z ) in powers of z we find A n = 1 √ π Γ (cid:0) n − (cid:1) Γ( n + 1) (cid:39) √ π n − / (9)In particular, the probabilities A n are given by , , , , , , , , , , for n = 1 , . . . , k initially infected individuals, infection processesoriginated with each individual are independent. Hencethe probability A ( k ) n that exactly n individuals catch theinfection before the epidemic is over can be expressedvia the probabilities A m ≡ A (1) m corresponding to theclassical situation with one initially infected individual: A ( k ) n = (cid:88) i + ... + i k = n A i . . . A i k (10)The generating function A ( k ) ( z ) = (cid:88) n ≥ k A ( k ) n z n (11)is therefore A ( k ) ( z ) = [ A ( z )] k = (cid:8) − √ − z (cid:9) k (12)This generating function encapsulates all A ( k ) n .Let us first look at the probabilities A ( k ) n for small k .Needless to say, A ( k ) n = 0 when n < k . When n ≥ k , onecan express A ( k ) n through the probabilities (9). Here area first few explicit formulas A (2) n = 2 A n A (3) n = 4 A n − A n − A (4) n = 8 A n − A n − A (5) n = 16 A n − A n − + A n − A (6) n = 32 A n − A n − + 6 A n − A (7) n = 64 A n − A n − + 24 A n − − A n − A (8) n = 128 A n − A n − + 80 A n − − A n − A (9) n = 256 A n − A n − + 240 A n − − A n − + A n − suggesting the general representation of A ( k ) n as a sum A ( k ) n = (cid:98) k − (cid:99) (cid:88) p =0 ( − p (cid:18) k − − pp (cid:19) k − − p A n − p (13)Here (cid:98) x (cid:99) is the largest integer ≤ x and (cid:0) ab (cid:1) are binomialcoefficients. Massaging the sum in Eq. (13) one obtainsa neat final formula A ( k ) n = k n − k Γ(2 n − k )Γ( n + 1 − k )Γ( n + 1) (14)The derivation of Eq. (14) from (13) is outlined in Ap-pendix A, where we also present a more straightforwardderivation of Eq. (14) from the generating function (12). / k0.0050.0100.0150.020 k A n ( k ) FIG. 1: Top to bottom: √ k A ( k ) n versus n/k for k = 1 , , π ) − / ( n/k ) − / in agreementwith Eq. (18). Note that substituting (9) into the recurrence (10) onecan directly compute k ≥ A ( k ) k = 12 k k ≥ A ( k ) k +1 = k k +2 k ≥ A ( k ) k +2 = k k +3 + k ( k − k +5 (15)These results are recovered from the general solution (14)thereby providing a consistency check. As another consis-tency check we note that (14) agrees with normalization (cid:88) n ≥ k A ( k ) n = 1 (16)The sequence A ( k ) n has a single peak located at n = k when k ≤
4, while for k ≥ n = ν ( k ) > k ,see Fig. 1. The sequence A ( k ) n grows from A ( k ) k = 2 − k tothe maximum b ( k ) := max { A ( k ) n | n ≥ k } (17)at n = ν ( k ), then decays and eventually approaches A ( k ) n (cid:39) k √ π n − / (18)This asymptotic behavior is straightforwardly deducedfrom (14). Using the general solution (14) together withthe Stirling formula one deduces the behaviors of b ( k )and ν ( k ) in the k → ∞ limit: ν ( k ) (cid:39) k , b ( k ) (cid:39) Bk − (19)with B = 3 e − / (cid:114) π = 0 . . . . (20) μ Φ FIG. 2: The scaled distribution Φ( µ ) versus the scaled out-break size µ given by Eq. (21). The general solution (14) approaches the scaling form A ( k ) n = 4 k Φ( µ ) , Φ( µ ) = π − / µ − / e − /µ (21)in the scaling limit n → ∞ , k → ∞ , µ = 4 nk = finite (22)The scaled distribution Φ( µ ) has a single peak (Fig. 2)and it vanishes faster than any power of µ when µ → (cid:104) n a (cid:105) = (cid:80) n ≥ k n a A ( k ) n diverge when a ≥ .To perform the summation in the a < range where themoments converge, let us first consider modified momentswith n a replaced by Γ( n +1) / Γ( n +1 − a ). Such momentsadmit an analytical expression (cid:88) n ≥ k Γ( n + 1)Γ( n + 1 − a ) A ( k ) n = Γ(1 − a )Γ(1 − a ) Γ(1 + k )Γ(1 − a + k ) (23)This result has been established by substituting A ( k ) n given by (14) into Eq. (23) and computing the sum usingidentities that can be found, e.g., in Ref. [32].When k (cid:29)
1, the modified moments are asymptoti-cally equal to the regular moments since Γ( n +1)Γ( n +1 − a ) → n a for n ≥ k (cid:29)
1. Therefore Eq. (23) simplifies to (cid:104) n a (cid:105) = (cid:88) n ≥ k n a A ( k ) n (cid:39) Γ(1 − a )Γ(1 − a ) k a (24)when a < . In the complementary a ≥ , the moments (cid:104) n a (cid:105) remain finite for finite populations, but diverge as N → ∞ . This divergent leading behavior of the momentsin finite populations is computed in the next section; inthe simplest case of k = 1, it is given by Eq. (27). III. OUTBREAK SIZE DISTRIBUTION:SCALING ANALYSIS
Consider the critical SIR process with k initially in-fected individuals, but in a finite population. Denote by N the size of the population and by A n ( k ; N ) the prob-ability that the size of the outbreak is n . In the infinite-population limit, A n ( k ; ∞ ) ≡ A ( k ) n is given by (14); itapproaches the scaling form (21)–(22) when k (cid:29) A n ( N ) ≡ A n (1; N ) forthe critical SIR process starting with a single infectedindividual has been investigated in Refs. [20–25]. Theprobability distribution A n ( N ) acquires a scaling form A n ( N ) = A n A ( ν ) (25)in the scaling limit n → ∞ , N → ∞ , ν = nN / = finite (26)This scaling was proposed in [21]. Simulations of verylarge systems with N comparable to the current humanpopulation [21, 25] are well fitted by the scaling form(25)–(26). Kessler and Shnerb [22] derived the scalingfunction A ( ν ); similar scaling functions have been rigor-ously studied in Refs. [23, 24].Recall that the moments (cid:104) n a (cid:105) with a ≥ diverge in theinfinite-population limit. For finite populations, thesemoments are finite and the scaling (25)–(26) implies thefollowing leading asymptotic behavior (cid:104) n a (cid:105) (cid:39) (cid:40) (9 π ) − / ln N a = C a N a − a > (27)with C a = (cid:90) ∞ dν √ πν ν a A ( ν ) (28)Two important properties of the epidemics are the av-erage and variance of the size of outbreaks. These quan-tities are defined by Eq. (2). When k = 1, the averageand the variance scale with population size as E ( N ) (cid:39) C N / , V ( N ) (cid:39) C N (29)Indeed, Eq. (27) shows that (cid:104) n (cid:105) (cid:29) (cid:104) n (cid:105) when N (cid:29)
1, soEq. (29) follows from (27), and the amplitudes are givenby (28). We have computed the amplitudes using theexact expression [22] for the scaling function to give C = 1 . . . . , C ≈ .
99 (30)In the general case, the distribution A n ( k, N ) dependson three variables. The interesting range is n ∼ N / , k ∼ N / (31)The first scaling law follows from Eq. (26), the second isan outcome of Eqs. (22) and (26). In the scaling region(31), the distribution A n ( k, N ) is expected to acquire ascaling form A n ( k, N ) = N − / A ( κ, ν ) (32)with ν = n/N / , see Eq. (26), and the scaled initialnumber of infected individuals κ = k/N / . More pre-cisely, (32) should hold in the scaling limit (26) and k → ∞ , N → ∞ , κ = kN / = finite (33)The normalization condition, (cid:80) n ≥ k A n ( k, N ) = 1, gives (cid:82) ∞ dν A ( κ, ν ) = 1 and explains the pre-factor in (32).The two-variable distribution A ( κ, ν ) is unknown, solet us discuss the average outbreak size which is the ba-sic quantity with simpler scaling behavior. The averageoutbreak size E k ( N ) depends on two variables and itsconjectural scaling behavior is E k ( N ) = N / E ( κ ) (34)The two scaling behaviors, (32) and (34), are compatibleif E ( κ ) = (cid:82) ∞ dν ν A ( κ, ν ).Let us probe the extremal behaviors of the scaled dis-tribution E ( κ ). To establish the small κ asymptotic, wenote that in the k (cid:28) N / regime, the infectious pro-cesses generated by each initially infected individual aremutually independent. Thus E k ( N ) (cid:39) C k N / when k (cid:28) N / (35)Comparing (34) and (35) we obtain E ( κ ) (cid:39) C κ as κ → C appearing in Eq. (4a), it is given by (30).To appreciate the large κ behavior of E ( κ ) we mentionthat E N ( N ) = N . This suggests that E ( N / ) ∼ N / ,from which E ( κ ) ∼ √ κ as κ → ∞ . This argument isheuristic. In Sec. IV, we derive the large κ asymptoticstated in Eq. (4a). In Appendix B, we present anotherelementary derivation based on the observation that thebehavior in the κ → ∞ limit is essentially deterministic.The deterministic analysis is significantly simpler thanthe exact approach, but it is not suitable for studyingfluctuations.Simulation results [26] are in a fairly good agreementwith the linear behavior of the scaled average outbreaksize in the small κ limit: E ( κ ) = C κ with C ≈ . C = 1 . . . . .In the large κ limit, numerical data in [26] were fitted to √ κ up to a logarithmic correction.We now turn to the variance of the size of the out-breaks. For sufficiently small k , we rely again on themutual independence of k infectious processes to deduce V k ( N ) (cid:39) C k N when k (cid:28) N / (36)The scaling region is given by (33). The natural scalingbehavior of the variance compatible with (36) is V k ( N ) = N / V ( κ ) (37) where V (( κ ) (cid:39) C κ when κ → κ → ∞ limit stated in Eq. (4b) is established in Sec. IV.Simulation results [26] support an algebraic decay inthe κ → ∞ limit: V ∼ κ − γ . The uncertainty [26] in themagnitude of the exponent is significant: γ = 0 . ± . γ = , andwe have also derived the amplitude: V (cid:39) (cid:112) /κ as statedin Eq. (4b). IV. OUTBREAK SIZE DISTRIBUTION: EXACTTREATMENT
The critical SIR process admits an exact treatment.Denote by s, i and r the population sizes in the suscep-tible, infected and recovered individuals. The entire sizeof the population is N , so s + i + r = N (38)Due to the constraint (38), the state of the process canbe described by any pair of variables s, i, r . We choose( i, x ) with x = N − s . For the critical SIR process, inthe interesting regime s is close to N , viz. N − s (cid:28) N ,and hence x is a more convenient variable than s . Theconstraint (38) shows that x = i + r , so x ≥ i .Infection and recovery events are symbolically( i, x ) → ( i + 1 , x + 1) rate i ( N − x ) /N (39a)( i, x ) → ( i − , x ) rate i (39b)Denote by t ( i, x ) the number of transitions from thestate ( i, x ) to termination. We are mostly interested in t ( k, k ), i.e., starting with k infected and no recovered.The process terminates at some state (0 , n ), where n isthe size of the outbreak. The rules (39a) and (39b) showthat the quantity i − x decreases by 1 in each transition.Thus starting at ( i, x ) = ( k, k ) gives i − x = − k − T after T transitions, and in particular n = k + t ( k, k )2 (40)The rates of the processes (39a) and (39b) imply thatthey occur with probabilities p + ( x ) = N − x N − x , p − ( x ) = N N − x (41)The stochastic transition time t ( i, x ) evolves accordingto the rules t ( i, x ) = (cid:40) t ( i + 1 , x + 1) prob p + ( x )1 + t ( i − , x ) prob p − ( x ) (42) A. Average number of transitions: Exact solution
Averaging (42) we find that T ( i, x ) = (cid:104) t ( i, x ) (cid:105) satisfies T = 1 + p + T ( i + 1 , x + 1) + p − T ( i − , x ) (43)To avoid cluttering of formulae, we write p ± ≡ p ± ( x ) and T ≡ T ( i, x ) when there is no confusion. The recurrence(43) should be solved subject to the boundary condition T (0 , x ) = 0 (44)The boundary-value problem (43)–(44) admits an ex-act solution T = i + 2( N − x ) − N − x (cid:88) j =1 (cid:18) NN + j (cid:19) i + N − x − j B ( N − x ) j ( N ) (45)with B pj ( N ) determined recurrently from B ( p ) j ( N ) = pp − j B ( p − j ( N ) , j = 1 , . . . , p − B ( p ) p ( N ) = 2 p − p − (cid:88) j =1 (cid:18) NN + j (cid:19) p − j B ( p ) j ( N ) (46)One can verify by direct substitution that Eq. (45) withamplitudes determined from the recurrent relations (46)satisfies (43)–(44). In Appendix C, we describe the ra-tionale behind the ansatz (45).Specializing Eq. (45) to i = x = k we obtain T ( k, k ) = 2 N − k − N − k (cid:88) j =1 (cid:18) NN + j (cid:19) N − j B ( N − k ) j ( N ) (47)The average size of an outbreak, E k ( N ) = [ k + T ( k, k )] / E k ( N ) = N − N − k (cid:88) j =1 (cid:18) NN + j (cid:19) N − j B ( N − k ) j ( N ) (48)An exact formula (48) for the average size of the out-break valid for arbitrary k and N was a pleasant out-come, yet we have not succeeded so far in extracting anasymptotic behavior of the sum in Eq. (48). There aretwo technical challenges: (i) the amplitudes determinedby the recurrence relations (46), are unwieldy; (ii) thesum in Eq. (48) involves N − k terms. One can find com-pact exact results when N − k = O (1), see Appendix C.Since k = O ( N / ) in the interesting range, the numberof terms in the sum in Eq. (48) is close to N (cid:29) E ( κ ) from the exact solution (46)–(48),and perhaps even sub-leading corrections to the leadingasymptotic behavior, Eq. (34). A numerical integrationof (46)–(48) may prove advantageous to simulations sinceone does not need to perform the averaging. Here wemerely demonstrate how to obtain explicit analytical re-sults for sufficiently small N . Solving the recurrence inthe top line Eq. (46) gives B ( p ) j = (cid:18) pj (cid:19) B j , B j ≡ B ( j ) j (49) N E ( N ) E ( N )2 TABLE I: The average size of an outbreak starting with oneinfected individual, E ( N ), and with two infected individuals, E ( N ), for N = 2 , . . . , E , E FIG. 3: Bottom: The average size of an outbreak for epi-demics starting with a single infected individual, E ( N ). Top:The average size of an outbreak for epidemics starting withtwo infected individuals, E ( N ). and therefore the bottom in Eq. (46) turns into a systemof linear equations p − (cid:88) j =1 (cid:18) NN + j (cid:19) p − j (cid:18) pj (cid:19) B j + B p = 2 p (50)The average size of an outbreak simplifies to E k = N − N − k (cid:88) j =1 (cid:18) NN + j (cid:19) N − j (cid:18) N − kj (cid:19) B j (51)with B j = B j ( N ) found by solving Eqs. (50).Exact results for E ( N ) and E ( N ) are given in Table Ifor N ≤
9, and plotted in Fig. 3 for N ≤
20. Even forsuch small populations, the exact results are rather closeto the asymptotic behavior (35).
B. Average number of transitions: Continuumtreatment
An often potent line of attack on asymptotic behav-ior of solution to discrete problems relies on employingcontinuum methods. In the present case we assume that T ( i, x ) is a smooth function of i and x , and we expand T ( i + 1 , x + 1) and T ( i − , x ) appearing in Eq. (43) upto the second order T ( i − , x ) = T − ∂ i T + ∂ i T T ( i + 1 , x + 1) = T + ∂ i T + ∂ x T + ∂ i T + ∂ x T + ∂ i ∂ x T (52)Here we use shorthand notation ∂ i = ∂/∂i , ∂ x = ∂/∂x ,etc. for the partial derivatives. Inserting (52) into (43)and keeping only dominant terms we obtain ∂ i T + ∂ x T + 2 = xN ∂ i T (53)Suppose that the scaling in the interesting range is i ∼ N α , x ∼ N β , T ∼ N γ (54)Plugging (54) into (53) we find that the terms in arecomparable only when α = and β = γ = . Thus were-scale the variables i = N / I, x = N / X (55)and the average number of transitions T ( i, x ) = N / T ( I, X ) (56)One can verify that terms not included in (53) are sub-dominant. For instance, computing the second derivatesgives ∂ i ∂ x T = O ( N − / ) and ∂ x T = O ( N − / ), sothese derivatives can indeed be dropped.The transformation (55)–(56) turns (53) into a partialdifferential equation (PDE) ∂ T ∂I + ∂ T ∂X + 2 = X ∂ T ∂I (57)for the re-scaled transition time T ( I, X ).We must solve Eq. (57) in the quadrant I ≥ X ≥
0. The boundary condition, Eq. (44), yields T (0 , X ) = 0 (58)Solving the boundary-value problem (57)–(58) is anintriguing challenge that we leave for the future. Herewe limit ourselves by a simpler problem of computingthe asymptotic behavior of T ( k, k ) when k (cid:29) N / .In the realm of the framework (55)–(58), we shouldlearn how to extract the large I behavior of T ( I, X ). Thiscan be done by noting that when I (cid:29)
1, the diffusionterm can be dropped from Eq. (57). Thus we arrive atthe first order PDE ∂ T ∂I − X ∂ T ∂X = 2 X (59)Introducing new variables u = I + X , v = I − X (60) we recast (59) into ∂ T ∂v = 1 √ u − v (61)The solution is T = − √ u − v + f ( u )with an arbitrary function f ( u ). The boundary condition(58) gives T = 0 when v = − u . This fixes f ( u ) = 2 √ u .Combining T = 2 √ u − √ u − v and (60) we obtain T ( I, X ) = 2 (cid:112) I + X − X (62)We want to determine T ( k, k ) = N / T ( κ, N − / κ ).Thus I = κ (cid:29) X = 0 as we always considerlarge populations, N (cid:29)
1. More precisely, setting X = 0amounts for a tacit assumption κ (cid:28) N / . Summarizing,our asymptotic results are valid in the range N / (cid:28) k (cid:28) N / (63)The upper and lower bounds are well separated when N / (cid:29)
1. This is satisfied for large populations, yet theconvergence may be slow as the effective small parameteris N − / .Thus T ( k, k ) = N / T ( κ,
0) when the bounds (63) areobeyed. Using Eq. (62) we get T ( k, k ) = √ kN whichwe insert into E k ( N ) = [ k + T ( k, k )] / E k ( N ) = √ kN . This completes the derivation of thelarge κ behavior announced in Eq. (4a). C. Variance
To derive the governing equations for the variance wefirst take the square of Eq. (42). Performing averagingwe find that T ( i, x ) = (cid:104) t ( i, x ) (cid:105) satisfies T = 1 + p + T ( i + 1 , x + 1) + p − T ( i − , x )+ 2 p + T ( i + 1 , x + 1) + 2 p − T ( i − , x ) (64)Again we shortly write p ± ≡ p ± ( x ) and T ≡ T ( i, x ).We now subtract the square of Eq. (43) from Eq. (64)and find that the variance V ( i, x ) = (cid:104) t ( i, x ) (cid:105) − (cid:104) t ( i, x ) (cid:105) (65)satisfies V ( i, x ) = p + V ( i + 1 , x + 1) + p − V ( i − , x )+ p + p − [ T ( i + 1 , x + 1) − T ( i − , x )] (66)We do not attempt to solve (66) and switch to the contin-uum treatment. Similar to (52) we expand the variance V ( i − , x ) = V − ∂ i V + ∂ i VV ( i + 1 , x + 1) = V + ∂ i V + ∂ x V + ∂ i V + ∂ x V + ∂ i ∂ x V (67)Plugging the expansions (52) and (67) into (66) and keep-ing only dominant terms we obtain ∂ i V + ∂ x V − xN ∂ i V + 2( ∂ i T ) = 0 (68)We use the same rescaled variables (55) as before, andseek the variance in the scaling form V ( i, x ) = N / V ( I, X ) (69)The transformation (55) and (69) turns (68) into ∂ V ∂I + 2 (cid:18) ∂ T ∂I (cid:19) = X ∂ V ∂I − ∂ V ∂X (70)When I (cid:29)
1, we can drop again the diffusion termfrom (70). We also use the asymptotic expression (62)for T ( I, X ) and arrive at the first order PDE
X ∂ V ∂I − ∂ V ∂X = 82 I + X (71)Using the variables (60) we recast (71) into ∂ V ∂v = 2 u √ u − v (72)which is integrated to give V = 4 u − (cid:2) √ u − √ u − v (cid:3) , or V ( I, X ) = 8 √ I + X − X I + X (73)Setting I = κ and X = 0 in Eq. (73) gives the large κ behavior: V ( κ,
0) = 4 (cid:112) /κ . Using Eq. (40) we find V k ( N ) = N / V ( κ, κ behavior an-nounced in Eq. (4b). V. DURATION OF OUTBREAKS
Some basic features of the duration of the outbreaks inthe critical SIR process in a finite system can be extractedfrom the temporal behaviors in the infinite-populationlimit. We first recall these infinite-population results inthe simplest case with one initially infected individual[16, 25]. The probability P i ( t ) to have i infected individ-uals at time t satisfies˙ P i = ( i − P i − − iP i + ( i + 1) P i +1 , i ≥ P = P (74b)where dot denotes the derivative with respect to time.A solution of an infinite set of equations (74a)–(74b)subject to the initial condition P i (0) = δ i, reads P i ( t ) = (1 + t ) − τ i − , i ≥ P ( t ) = τ ≡ t t (75b) This soultion can be verified by a direct substitution, orderived using e.g. generating function techniques [see Ap-pendix D for details]. The probability that the outbreakis still alive at time t , is P ( t ) = (cid:88) i ≥ P i ( t ) = 1 − P ( t ) = 11 + t (76)The average number of infected individuals in outbreakswhich are still alive at time t is therefore (cid:104) i (cid:105) = (cid:80) iP i ( t ) (cid:80) P i ( t ) = 1 + t (77)A general solution of Eqs. (74a)–(74b) describing aninfinite-population limit subject to an arbitrary numberof initially infected individuals is somewhat cumbersome,it is presented in Appendix D.In a finite population, the infection process eventu-ally comes to an end. To estimate heuristically this finaltime t f one uses (77) to express [21] the final size of theoutbreak through the final time: n f ∼ (cid:82) t f dt (cid:104) i (cid:105) ∼ t .The maximal outbreak size scales as n ∗ ∼ N / (see e.g.[21, 22, 25]) and hence the maximal duration is t ∗ ∼ n / ∗ ∼ N / (78)The average duration of the outbreak is formally E [ t ] = (cid:90) ∞ dt t (cid:18) − dPdt (cid:19) = (cid:90) ∞ dt tP ( t ) (79)Recalling that P ( t ) = (1+ t ) − in the infinite-populationlimit (equivalently, for the critical branching process), wenotice that the integral in (79) diverges. We should use,however, the finite upper limit given by (78). This leadsto an estimate E [ t ] (cid:39) (cid:90) √ N dt t (1 + t ) (cid:39)
13 ln N (80)where the subscript reminds that the process begins witha single infected individual. The logarithmic growth ofthe average duration time was predicted by Ridler-Rowemany years ago [19], albeit with incorrect amplitude; thecorrect amplitude 1 / E [ t ] = ln N + c + o (1 /N ) (81)Since the logarithm is a slowly growing function, the sub-leading constant term c significantly contributes to theaverage duration. The computation of the sub-leadingterm requires much more comprehensive analysis thanwhat we have used so far.If the number of initially infected individuals is suffi-ciently small, k (cid:28) N / , one can generalize the predic-tion (80) for the average duration of an outbreak with-out using a complete solution of the infinite-populationlimit (Appendix D), it suffices to rely on the indepen-dence of infection processes generated by each initiallyinfected individual. The probability that the infectionis over at time t is P k , with P given by (75b). Thus − dP k /dt = kP k − / (1 + t ) is the probability densitythat the infection is eradicated at time t , from which E k [ t ] (cid:39) (cid:90) √ N dt kt k (1 + t ) k +1 (cid:39) k N (82)implying that the average duration of the outbreak ex-hibits a simple logarithmic scaling with amplitude pro-portional to the initial number of infected individuals.Equation (82) suggests a plausible scaling form E k [ t ] = N / ln N κ ) (83)with Θ( κ ) (cid:39) κ when κ → κ = k/N / . The N − dependent pre-factor, how-ever, contains an additional logarithmic factor. The un-expected logarithmic factor was originally observed insimulations [26]. In Appendix B, we probe the large κ using a deterministic approach. This analysis supportsthe scaling form (83). Overall, the extremal behaviors ofthe scaled average time areΘ( κ ) = (cid:40) κ when κ (cid:28) (cid:112) /κ when κ (cid:29) V [ t ] ∼ (cid:90) √ N dt t (1 + t ) ∼ N / (85)but not an amplitude. Independence gives V [ t ] ∼ kN / when k (cid:28) N / , and hence the hypothetical scaling formof the variance is V k [ t ] = N / Θ ( κ ) (86)with Θ ( κ ) ∼ κ when κ →
0. In contrast to the aver-age, the scaled form for the variance does not contain alogarithmic factor.Simulations [26] support the scaling law (86) and thelinear small κ behavior. Simulations also suggest [26] that the scaled distribution Θ ( κ ) is inversely propor-tional to κ in the large κ limit. ThusΘ ( κ ) ∼ (cid:40) κ κ → κ − κ → ∞ (87) VI. CONCLUSIONS
We have studied the critical SIR process starting witha large number of initially infected individuals, k (cid:29) k scalesas a cubic root of the population size, k ∼ N / . Thecritical SIR process exhibits large fluctuations, so we haverelied on the stochastic formulation. We have treated theproblem using a combination of exact calculations andasymptotic methods. We have focused on the size andduration of the outbreaks, more precisely on the averageand variance of these quantities. The analysis of the sizeof outbreaks is rather detailed, Secs. III–IV.Our analytical and asymptotic predictions qualita-tively agree with simulation results [26]. Whenever thereis a discrepancy between simulation results and theoret-ical predictions, it seems that data may be fitted usingslightly different scaling forms, viz. simpler than the fitsused in Ref. [26]. The chief reason for subtle behaviorsare algebraic tails, the average size and duration of out-breaks are especially sensitive to these tails.We have found an exact expression for the average sizeof outbreaks valid for arbitrary N and k . Extractingthe scaled size distribution E ( κ ) and sub-leading correc-tions to the leading behavior (34) is left for future work.The derivation of the exact average size in Sec. IV can beprobably generalized to establish the variance and cumu-lants (perhaps even the cumulant generating function).Continuum methods allow, in principle, to determine thescaling functions, one should be able to solve linear PDEswith non-constant coefficients. Acknowledgments.
I want to thank Ginestra Bianconiand Sid Redner for collaboration on similar problems,G. Bianconi and F. Radicchi for sending a preliminaryversion of Ref. [26], and Boston University Network groupfor discussions and encouragement. I am also grateful tothe referee who has given a simple derivation of Eq. (14)presented in Appendix A. [1] M. McNeill,
Plagues and People (Anchor Books, NewYork, 1989).[2] M. Oldstone,
Viruses, Plagues, and History (Oxford Uni-versity Press, Oxford, 1998).[3] O. J. Benedictow,
The Black Death 1346–1353: TheComplete History (Boydell Press, Wiltshire, 2012).[4] D. Bernoulli, M´em. Math. Phys. Acad. Roy. Sci., Paris, 1 (1760).[5] K. Dietza and J. A. P. Heesterbeek, Math. Biosciences , 1 (2002).[6] D. P. Maki and M. Thompson,
Mathematical Models andApplications, with emphasis on the social, life, and man-agement sciences (Englewood Cliffs, N.J., Prentice-Hall,1973). [7] M. A. Ludwig, The Giant Black Book of ComputerViruses (American Eagle Publications Inc., Show Low,1998).[8] C. Castellano, S. Fortunato, and V. Loreto, Rev. Mod.Phys. 81, 591 (2009).[9] P. L. Krapivsky, S Redner, and D. Volovik, J. Stat. Mech.P12003 (2011).[10] A. G. McKendrick, Proc. Edin. Math. Soc. , 98 (1926).[11] W. O. Kermack and A. G. McKendrick, Proc. Roy. Soc.A , 700 (1927).[12] R. Anderson and R. May, Infectious Diseases: Dynamicsand Control (Oxford University Press, Oxford, 1991).[13] H. W. Hethcote, SIAM Rev. , 599 (2000).[14] J. D. Murray, Mathematical Biology. I. An Introduction (Springer-Verlag, New York, 2002).[15] N. T. J. Bailey, Biometrika , 193 (1950); ibid , 177(1953).[16] N. T. J. Bailey, The Mathematical Theory of InfectiousDiseases (Oxford University Press, Oxford, 1987).[17] H. Andersson and T. Britton,
Stochastic Epidemic Mod-els and Their Statistical Analysis (New York, Springer,2000).[18] P. L. Krapivsky, S. Redner and E. Ben-Naim,
A KineticView of Statistical Physics (Cambridge, Cambridge Uni-versity Press, 2010).[19] C. J. Ridler-Rowe, J. Appl. Prob. , 19 (1967).[20] A. Martin-L¨of, J. Appl. Probab. , 671 (1998).[21] E. Ben-Naim and P. L. Krapivsky, Phys. Rev. E ,050901(R) (2004).[22] D. A. Kessler and N. M. Shnerb, Phys. Rev. E ,010901(R) (2007).[23] L. F. Gordillo, S. A. Marion, A. Martin-L¨of, and P. E.Greenwood, Bull. Math. Biol. , 589 (2008).[24] R. van der Hofstad, A. J. E. M. Janssen, and J. S. H.Leeuwaarden, Adv. Appl. Probab. , 706 (2010).[25] E. Ben-Naim and P. L. Krapivsky, Eur. Phys. J. B ,145 (2012).[26] F. Radicchi and G. Bianconi, arXiv:2007.15034.[27] W. Feller, An Introduction to Probability Theory and ItsApplications , Vol. I, 3rd edn. (John Wiley, New York,1968).[28] T. E. Harris,
The Theory of Branching Processes (Dover,New York, 1989).[29] K. B. Athreya and P. E. Ney,
Branching Processes (DoverPublications, Inc., Mineola, New York, 2004).[30] M. Kimmel and D. Axelrod,
Branching Processes in Bi-ology (Springer, New York, 2002).[31] P. Haccou, P. Jagers, and V. A. Vatutin,
Branching pro-cesses: variation, growth, and extinction of populations (Cambridge University Press, New York, 2005).[32] R. L. Graham, D. E. Knuth, and O. Patashnik,
Con-crete Mathematics: A Foundation for Computer Science (Reading, Mass.: Addison-Wesley, 1989).[33] The same leading behavior (80) is obtained if one woulduse the upper bound C √ N with arbitrary C ; due to thelogarithmic divergence of the integral, only the scaling ofthe upper bound with N is required for extracting theleading behavior.[34] S. N. Majumdar and R. M. Ziff, Phys. Rev. Lett. ,050601 (2008). Appendix A: Derivation of Eq. (14)
Plugging (9) into the sum in Eq. (13), one reduces thesum to a hypergeometric series [32]. Using identities in-volving hypergeometric series [32], as well as identitiesinvolving the gamma function (particularly, the duplica-tion formula), one computes the sum in (13) and arrivesat the announced expression (14) for A ( k ) n . This deriva-tion relies on identities for the hypergeometric series thatare little known. The intermediate result, Eq. (13), alsorequires a rather technical derivation.We now present an alternative derivation of Eq. (14)directly from the generating function (12). Employingthe Cauchy theorem, reduces the problem to computingthe integral A ( k ) n = (cid:73) C dz π i (cid:8) − √ − z (cid:9) k z n +1 (A1)over a small simple counter-clockwise contour C aroundthe origin in the complex z plane. To simplify the com-putation of the integral, let us try to change the variable z so that the square root √ − z would become a linearfunction of a novel complex variable. The change of vari-able, 1 − z = (2 v − , achieves this goal and reduces theintegral in (A1) to A ( k ) n = (cid:73) C dv π i (1 − v ) − ( n +1 − k ) n − k − (cid:20) v n − v n +1 (cid:21) (A2)over the appropriate contour C . The integral in (A2) isfound without a calculation, it suffices to use the residuetheorem and the binomial theorem [32](1 − v ) − m = (cid:88) (cid:96) ≥ Γ( m + (cid:96) )Γ( m ) Γ( (cid:96) + 1) v (cid:96) (A3)The final outcome is Eq. (14).The quantity A ( k ) n admits an interesting interpretationin terms of records of a one-dimensional random walkdiscrete in time and with an arbitrary symmetric andcontinuous jump distribution [34]. Namely, A ( k ) n is theprobability that such a random walk has k maxima in n time steps, with the last position being the maximum.Formulas similar to (12) and (14) appear in the contextof universal record statistics of random walks [34].It would be interesting to extend the relation betweenthe outbreak size distribution and the records of a one-dimensional discrete in time random walk to more com-plicated versions of the critical branching process, suchas the critical branching process with triplication ratherthan duplication or the discrete time critical branchingprocess. Appendix B: Deterministic treatment
For the SIR process starting with a single infected in-dividual, the deterministic framework is applicable only1in the super-critical regime if the process enters the run-away regime. In this situation, the epidemic ends only af-ter a finite fraction of the population [see Eq. (1)] catchesthe disease, and this fraction can be determined using thedeterministic framework. If the initial number of infectedindividuals is very large, the deterministic framework al-ways applies and one can describe the SIR process viathe system of differential equations for the densities S ( t )of susceptible, I ( t ) of infected, and R ( t ) of recovered in-dividuals.For the critical SIR process, R = 1, the deterministicrate equations read ˙ S = − SI (B1a)˙ I = − I + SI (B1b)˙ R = I (B1c)Equations (B1a)–(B1c) are consistent with the conserva-tion law S ( t ) + I ( t ) + R ( t ) = 1. The initial condition S (0) = 1 − (cid:15), I (0) = (cid:15), R (0) = 0 (B2)Below we always assume that N − (cid:28) (cid:15) (cid:28) (cid:15)N , is large and the deterministic framework should beasymptotically correct. Setting (cid:15) (cid:28) (cid:15) (cid:28) S as a time variable. Dividing (B1b) by (B1a)and integrating gives I = 1 − S + ln S − (cid:15) (B4)Similarly dividing (B1c) by (B1a) and integrating yields R = − ln S − (cid:15) (B5)The outbreak ends when I f = 0. Equation (B4) showsthat the final fraction of susceptible is implicitly deter-mined by 1 − S f = − ln S f − (cid:15) (B6)from which R f = 1 − S f (cid:39) √ (cid:15) (B7)when (cid:15) (cid:28)
1. Since R f = N − E k ( N ) and (cid:15) = k/N ,we can re-write (B7) as E k ( N ) = √ kN , leading to theasymptotic Ψ( κ ) = √ κ when κ (cid:29)
1. The applicability of the deterministic framework forsufficiently large k is intuitively clear. We have deter-mined the crossover k ∼ N / when the stochastic ef-fects become important using the stochastic framework.It would be interesting to deduce the crossover in therealm of the deterministic framework.The applicability of the deterministic framework onthe final stage with a few remaining infected individualsis questionable. We know, however, that for the criti-cal SIR process starting from a single infected individ-ual, large outbreaks are of the order of N / . Thus when κ (cid:29)
1, the outbreak size is N / √ κ + O ( N / ). The firstterm describes the deterministic stage, while the secondterm accounts for the stochastic final stage. The deter-ministic first term dominates when κ (cid:29)
1. Therefore inthis regime, one can use the deterministic framework.To estimate the average duration time using the de-terministic framework, we combine (B1a) and (B4) andobtain E [ t ] = (cid:90) S (0) S f dσσ (cid:104) − σ + ln σ − (cid:15) (cid:105) (B8)Writing σ = 1 − √ (cid:15) y , expanding the denominator inthe integral in Eq. (B8) in powers of (cid:15) (cid:28)
1, and recallingthat 1 − S (0) = (cid:15) and 1 − S f (cid:39) √ (cid:15) , we obtain E k [ t ] (cid:39) (cid:114) (cid:15) (cid:90) dy (cid:15) − y (cid:39) ln(8 /(cid:15) ) √ (cid:15) = N / ln(8 N / /κ ) √ κ (cid:39) N / ln N (cid:114) κ (B9)(We have taken into account that (cid:15) = k/N = κ/N / .)This asymptotic is consistent with the scaling form (83)and gives the announced large κ asymptotic in Eq. (84). Appendix C: Exact results
In Sec. IV, we have used asymptotic methods to deter-mine the leading behaviors. It is possible to derive ex-act results for the average outbreak size. Unfortunately,these results are expressed through the sum of a largenumber of variables with complicated individual terms.Here we explain how one could guess these results; theverification of the guess (45) is straightforward. We havemade the guess (45) by establishing a few explicit exactresults. These exact results describe a non-interesting re-gion where the initial number of susceptible individuals issmall, so their virtue is that they simplify the guesswork.The average number of transitions T ( i, x ) from thestate ( i, x ) to termination satisfies the recurrence (43)and the boundary condition (44). We start by solving(43)–(44) in the simplest cases when N − x = O (1). First,2we notice that T ( i, N ) = i (C1)Indeed, only recovery events are possible if there are nosusceptible. One can also formally derive (C1) by special-izing (43) to x = N . This gives T ( i, N ) = 1+ T ( i − , N )which in conjunction with T (0 , N ) = 0 lead to (C1).Similarly we specialize Eq. (43) to x = N − T ( i, N −
1) = 1 + NN + 1 T ( i − , N −
1) + i + 1 N + 1Solving this recurrence subject to T (0 ,
1) = 0 one obtains T ( i, N −
1) = i + 2 − (cid:18) NN + 1 (cid:19) i (C2)Specializing (43) to x = N − T ( i, N −
2) = 1 + NN + 2 T ( i − , N − N + 2 (cid:34) i + 3 − (cid:18) NN + 1 (cid:19) i +1 (cid:35) from which T ( i, N −
2) = i + 4 − (cid:18) NN + 1 (cid:19) i +1 − N + 1 (cid:18) NN + 2 (cid:19) i (C3)Similarly we compute T ( i, N −
3) = i + 6 − (cid:18) NN + 1 (cid:19) i +2 − N + 1 (cid:18) NN + 2 (cid:19) i +1 −
12 + 18 N ( N + 1) ( N + 2) (cid:18) NN + 3 (cid:19) i (C4)We are interested in T ( k, k ). Specializing Eq. (C1) to i = N yields T ( N, N ) = N . Similarly from (C2)–(C4)we extract T ( N − , N −
1) = N + 1 − (cid:18) NN + 1 (cid:19) N − T ( N − , N −
2) = N + 2 − (cid:18) NN + 1 (cid:19) N − − N + 1 (cid:18) NN + 2 (cid:19) N − T ( N − , N −
3) = N + 3 − (cid:18) NN + 1 (cid:19) N − − N + 1 (cid:18) NN + 2 (cid:19) N − −
12 + 18 N ( N + 1) ( N + 2) (cid:18) NN + 3 (cid:19) N − Looking at the above expressions for T ( N − p, N − p )with p = 0 , , ,
3, one guesses the general formula T ( N − p, N − p ) = N + p − p (cid:88) j =1 (cid:18) NN + j (cid:19) N − j B ( p ) j (C5)This is exactly (47) in different notation. More generally,Eqs. (C1)–(C4) suggest our chief ansatz, Eq. (45), whichis then straightforwardly verified.The amplitudes B ( p ) j ( N ) are rather simple for small p ,but quickly become cumbersome. Here are a few seriesof the amplitudes extracted from (46): B ( p )1 = 2 pB ( p )2 = 2 p ( p − N + 1 B ( p )3 = p ( p − p − N )( N + 1) ( N + 2) (C6)Therefore when p = O (1) is fixed T ( N − p, N − p ) = N + p (cid:18) − e (cid:19) − p ( p − e − N + 1 − p ( p − p − e − ( N + 1) + O ( N − )when N (cid:29)
1. The asymptotic behavior of the averagesize of the outbreak is therefore E N − p ( N ) = N − pe − p ( p − e N − + O ( N − ) (C7) Appendix D: The distribution P i ( t ) We want to solve Eqs. (74a)–(74b) subject to the initialcondition P i ( t = 0) = δ i,k . Using the generating function g ( z, t ) = (cid:88) i ≥ P i ( t ) z i (D1)we reduce the infinite system (74a)–(74b) of ordinary dif-ferential equations to a partial differential equation ∂ t g = (1 − z ) ∂ z g (D2)Introducing the auxiliary variable ζ = (1 − z ) − we recast(D2) into ( ∂ t − ∂ ζ ) g = 0 which is solved to yield g ( z, t ) = G ( t + ζ ) (D3)The function G is fixed by the initial condition G ( ζ ) = g ( z ) (D4)In terms of the original variable z , the solution (D3)–(D4)becomes g ( z, t ) = g (cid:18) − − z t − tz (cid:19) (D5)3This solution is valid for an arbitrary initial condition. Itis useful to rewrite this solution in terms of the reducedgenerating function P ( z, t ) = (cid:88) m ≥ P m ( t ) z m = g ( z, t ) − g (0 , t ) (D6)accounting only for the active part. One gets P ( z, t ) = g (cid:18) − − z t − tz (cid:19) − g ( τ ) (D7)In the classical case P i (0) = δ i, we get g ( z ) = z , so(D7) reduces to P ( z, t ) = 1 − τ − − z t − tz = 1(1 + t ) z − τ z (D8)Expanding (D8) we recover (75a).If P i (0) = δ i,k , we get g ( z ) = z k and P ( z, t ) = (cid:18) − − z t − tz (cid:19) k − τ k (D9)There are no infected with probability P = g ( τ ) = τ k , τ = t t (D10) Generally the probabilities P i ( t ) are obtained by expand-ing (D9) in powers of z . Explicit formulas(1 + t ) k +1 P = (cid:18) k (cid:19) t k − (1 + t ) k +2 P = (cid:18) k (cid:19) t k + (cid:18) k (cid:19) t k − (1 + t ) k +3 P = (cid:18) k (cid:19) t k +1 + 2 (cid:18) k (cid:19) t k − + (cid:18) k (cid:19) t k − for i = 1 , , P i = i (cid:88) a =1 (cid:18) ka (cid:19)(cid:18) i − a − (cid:19) t k + i − a (1 + t ) k + ii