Spatial extent of an outbreak in animal epidemics
Eric Dumonteil, Satya N. Majumdar, Alberto Rosso, Andrea Zoia
SSpatial extent of an outbreak in animal epidemics
Eric Dumonteil ∗ , Satya N. Majumdar † , Alberto Rosso † , and Andrea Zoia ∗ , ∗ CEA/Saclay, DEN/DM2S/SERMA/LTSD, 91191 Gif-sur-Yvette Cedex, France, and † CNRS - Universit´e Paris-Sud, LPTMS, UMR8626, 91405 Orsay Cedex, FranceSubmitted to Proceedings of the National Academy of Sciences of the United States of America
Characterizing the spatial extent of epidemics at the outbreak stageis key to controlling the evolution of the disease. At the outbreak,the number of infected individuals is typically small, so that fluctu-ations around their average are important: then, it is commonly as-sumed that the susceptible-infected-recovered (SIR) mechanism canbe described by a stochastic birth-death process of Galton-Watsontype. The displacements of the infected individuals can be modelledby resorting to Brownian motion, which is applicable when long-range movements and complex network interactions can be safelyneglected, as in case of animal epidemics. In this context, the spa-tial extent of an epidemic can be assessed by computing the convexhull enclosing the infected individuals at a given time. We derive theexact evolution equations for the mean perimeter and the mean areaof the convex hull, and compare them with Monte Carlo simulations.
Epidemics | Branching Brownian motion | Convex hull M odels of epidemics traditionally consider three classes ofpopulations, namely, the susceptibles (S), the infected(I), and the recovered (R). This provides the basis of the so-called SIR model [1, 2], a fully connected mean-field modelwhere the population sizes of the three species evolve withtime t via the coupled nonlinear equations: dS/dt = − βIS ; dI/dt = βIS − γI and dR/dt = γI . Here γ is the rate at whichan infected individual recovers and β denotes the rate at whichit transmits the disease to a susceptible [4, 5, 6]. In the sim-plest version of these models, the recovered can not be infectedagain. These rate equations conserve the total population size I ( t ) + S ( t ) + R ( t ) = N , and one assumes that initially thereis only one infected individual and the rest of the populationis susceptible: I (0) = 1, S (0) = N −
1, and R (0) = 0. Ofparticular interest is the outbreak stage, i.e., the early timesof the epidemic process, when the susceptible population ismuch larger than the number of infected or recovered. Dur-ing this regime, for large N , the susceptible population hardlyevolves and stays S ( t ) ≈ N , so that nonlinear effects can besafely neglected and one can just monitor the evolution of the infected population alone : dI/dt ≈ ( βN − γ ) I ( t ). Thus, the ul-timate fate of the epidemics depends on the key dimensionlessparameter R = βN/γ , which is called the reproduction rate.If R > R < R = 1 the infected population remainsconstant [7, 8, 9].This basic deterministic SIR has been generalized to a va-riety of both deterministic as well as stochastic models, whosedistinct advantages and shortcomings are discussed at lengthin [3, 10, 11]. Generally speaking, stochastic models are moresuitable in presence of a small number of infected individuals,when fluctuations around the average may be relevant [3, 10].During the outbreak of epidemics, the infected population istypically small: in this regime, the evolution can be modeledby resorting to a stochastic birth-death branching process ofthe Galton-Watson type for the number of infected [3, 10, 11],where each infected individual transmits the disease to an-other individual at rate βN and recovers at rate γ . The epi-demic may become endemic for R >
1, becomes extinct for R <
1, whereas for R = 1 fluctuations are typically longlived and completely control the time evolution of the infectedpopulation [4, 5, 6]. How far in space can an epidemic spread? Branching pro-cesses alone are not sufficient to describe an outbreak, andspatial effects must necessarily be considered [1, 5, 12, 13, 14].Quantifying the geographical spread of an epidemic is closelyrelated to the modelling of the population displacements.Brownian motion is often considered as a paradigm for de-scribing the migration of individuals, despite some well-knownshortcomings: for instance, finite speed effects and preferen-tial displacements are neglected. Most importantly, a num-ber of recent studies have clearly shown that individuals ge-ographically far apart can actually be closely related to eachother through the so-called small-world connections, such asair traffic, public transportation and so on: then, the spreadof epidemics among humans can not be realistically modelledwithout considering these complex networks of interconnec-tions [15, 16, 17, 18]. Nonetheless, Brownian motion providesa reasonable basis for studying disease propagation in animalsand possibly in plants (here, pathogen vectors are insects) [4].While theoretical models based on branching Brownianmotion have provided important insights on how the pop-ulation size grows and fluctuates with time in a given do-main [1, 5, 13, 14], another fundamental question is how thespatial extension of the infected population evolves with time .Assessing the geographical area travelled by a disease is keyto the control of epidemics, and this is especially true at theoutbreak, when confinement and vaccination could be mosteffective. One major challenge in this very practical field ofdisease control is how to quantify the area that needs to bequarantined during the outbreak. For animal epidemics, thisissue has been investigated experimentally, for instance in thecase of equine influenza [19]. The most popular and widelyused method for this consists in recording the set of positionsof the infected animals and, at each time instant, construct a convex hull, i.e., a minimum convex polygon surrounding thepositions (Fig. 1; for a precise definition of the convex hull,see below). The convex hull at time t then provides a roughmeasure of the area over which the infections have spread upto time t . The convex hull method is also used to estimatethe home range of animals, i.e., the territory explored by aherd of animals during their daily search for food [20, 21].In this paper, we model the outbreak of an epidemic asa Galton-Watson branching process in presence of Brownianspatial diffusion. Despite infection dynamics being relativelysimple, the corresponding convex hull is a rather complexfunction of the trajectories of the infected individuals up totime t , whose statistical properties seem to be a formidableproblem. Our main goal is to characterize the time evolution Reserved for Publication Footnotes
Issue Date
Volume Issue Number 1 – a r X i v : . [ q - b i o . P E ] M a r f the convex hull associated to this process, in particular itsmean perimeter and area.The rest of the paper is organized as follows. We first de-scribe precisely the model and summarize our main results.Then, we provide a derivation of our analytical findings, sup-ported by extensive numerical simulations. We conclude withperspectives and discussions. Some details of the computa-tions are relegated to the Supplementary Material. The model and the main results
Consider a population of N individuals, uniformly distributedin a two dimensional plane, with a single infected at the ori-gin at the initial time. At the outbreak, it is sufficient to keeptrack of the positions of the infected, which will be marked as‘particles’. The dynamics of the infected individuals is gov-erned by the following stochastic rules. In a small time interval dt , each infected alternatively (i) recovers with probability γ dt . This corresponds to thedeath of a particle with rate γ . (ii) infects, via local contact , a new susceptible individual fromthe background with probability b dt . This corresponds to thebirth of a new particle that can subsequently diffuse. Theoriginally infected particle still remains infected, which meansthat the trajectory of the originally infected particle branchesinto two new trajectories. The rate b replaces the rate βN inthe SIR or the Galton-Watson process mentioned before. (iii) diffuses with diffusion constant D with probability 1 − ( γ + b ) dt . The coordinates { x ( t ) , y ( t ) } of the particle getupdated to the new values { x ( t ) + η x ( t ) dt, y ( t ) + η y ( t ) dt } ,where η x ( t ) and η y ( t ) are independent Gaussian white noiseswith zero mean and correlators (cid:104) η x ( t ) η x ( t (cid:48) ) (cid:105) = 2 Dδ ( t − t (cid:48) ), (cid:104) η y ( t ) η y ( t (cid:48) ) (cid:105) = 2 Dδ ( t − t (cid:48) ) and (cid:104) η x ( t ) η y ( t (cid:48) ) (cid:105) = 0.The only dimensionless parameter in the model is the ratio R = b/γ , i.e, the basic reproduction number.Consider now a particular history of the assembly of thetrajectories of all the infected individuals up to time t , start-ing from a single infected initially at the origin (see Fig. 1).For every realization of the process, we construct the asso-ciated convex hull C . To visualize the convex hull, imaginestretching a rubber band so that it includes all the points ofthe set at time t inside it and then releasing the rubber band.It shrinks and finally gets stuck when it touches some pointsof the set, so that it can not shrink any further. This finalshape is precisely the convex hull associated to this set.In this paper, we show that the mean perimeter (cid:104) L ( t ) (cid:105) andthe mean area (cid:104) A ( t ) (cid:105) of the convex hull are ruled by two cou-pled nonlinear partial differential equations that can be solvednumerically for all t (see Fig. 2). The asymptotic behavior forlarge t can be determined analytically for the critical ( R = 1),subcritical ( R <
1) and supercritical ( R >
1) regimes. Inparticular, in the critical regime the mean perimeter saturates to a finite value as t → ∞ , while the mean area diverges log-arithmically for large t (cid:104) L ( t → ∞ ) (cid:105) = 2 π (cid:114) Dγ + O ( t − / ) [1] (cid:104) A ( t → ∞ ) (cid:105) = 24 πD γ ln t + O (1) . [2] This prediction seems rather paradoxical at a first glance.How can the perimeter of a polygon be finite while its areais divergent? The resolution to this paradox owes its originprecisely to statistical fluctuations. The results in Eqs. [1]and [2] are true only on average. Of course, for each sample,the convex hull has a finite perimeter and a finite area. How- ever, as we later show, the probability distributions of theserandom variables have power-law tails at long time limits. Forinstance, while Prob(
L, t → ∞ ) ∼ L − for large L (thus lead-ing to a finite first moment), the area distribution behaves asProb( A, t → ∞ ) ∼ A − for large A . Hence the mean area isdivergent as t → ∞ (see Fig. 2).When R (cid:54) = 1, the evolution of the epidemic is controlledby a characteristic time t ∗ , which scales like t ∗ ∼ | R − | − .For times t < t ∗ the epidemic behaves as in the critical regime.In the subcritical regime, for t > t ∗ the quantities (cid:104) L ( t ) (cid:105) and (cid:104) A ( t ) (cid:105) rapidly saturate and the epidemic goes eventually to ex-tinction. In contrast, in the supercritical regime (which is themost relevant for virulent epidemics that spread fast), a newtime-dependent behavior emerges when t > t ∗ , since there ex-ists a finite probability (namely 1 − /R ) that epidemic nevergoes to extinction (Fig. 5). More precisely, we later show that (cid:104) L ( t (cid:29) t ∗ ) (cid:105) = 4 π (cid:18) − R (cid:19) (cid:112) D γ ( R − t [3] (cid:104) A ( t (cid:29) t ∗ ) (cid:105) = 4 π (cid:18) − R (cid:19) D γ ( R − t . [4] The ballistic growth of the convex hull stems from an un-derlying travelling front solution of the non-linear equationgoverning the convex hull behavior. Indeed, the prefactorof the perimeter growth is proportional to the front veloc-ity v ∗ = 2 (cid:112) D γ ( R − The statistics of the convex hull
Characterizing the fluctuating geometry of C is a formidabletask even in absence of branching ( b = 0) and death ( γ = 0),i.e., purely for diffusion process in two dimensions. Majorrecent breakthroughs have nonetheless been obtained for dif-fusion processes [23, 24] by a clever adaptation of the Cauchy’sintegral geometric formulae for the perimeter and area of anyclosed convex curve in two dimensions. In fact, the prob-lem of computing the mean perimeter and area of the convexhull of any generic two dimensional stochastic process can bemapped, using Cauchy’s formulae, to the problem of comput-ing the moments of the maximum and the time at which themaximum occurs for the associated one dimensional compo-nent stochastic process [23, 24]. This was used for computing,e.g., the mean perimeter and area of the convex hull of a twodimensional regular Brownian motion [23, 24] and of a twodimensional random acceleration process [25].Our main idea here is to extend this method to computethe convex hull statistics for the two dimensional branchingBrownian motion. Following this general mapping and usingisotropy in space (see Supplementary Materials), the averageperimeter and area of the convex hull are given by (cid:104) L ( t ) (cid:105) = 2 π (cid:104) x m ( t ) (cid:105) [5] (cid:104) A ( t ) (cid:105) = π (cid:2) (cid:104) x m ( t ) (cid:105) − (cid:104) y ( t m ) (cid:105) (cid:3) , [6] where x m is the maximum displacement of our two-dimensional stochastic process in the x direction up to time t , t m is the time at which the maximum displacement along x direction occurs and y ( t m ) is the ordinate of the process at t m , i.e., when the displacement along the x direction is maxi-mal. A schematic representation is provided in Fig. 4, wherethe global maximum x m is achieved by one single infected in-dividual, whose path is marked in red. A crucial observationis that the y component of the trajectory connecting O to his red path is a regular one dimensional Brownian motion.Hence, given t m and t , clearly (cid:104) y ( t m ) (cid:105) = 2 D (cid:104) t m (cid:105) . Therefore, (cid:104) A ( t ) (cid:105) = π (cid:2) (cid:104) x m ( t ) (cid:105) − D (cid:104) t m ( t ) (cid:105) (cid:3) . [7] Equations [5] and [7] thus show that the mean perimeter andarea of the epidemics outbreak are related to the extremestatistics of a one dimensional branching Brownian motionwith death. Indeed, if we can compute the joint distribution P t ( x m , t m ), we can in turn compute the three moments (cid:104) x m (cid:105) , (cid:104) x m (cid:105) and (cid:104) t m (cid:105) that are needed in Eqs. [5] and [7]. We showbelow that this can be performed exactly. The convex hull perimeter and the maximum x m . For theaverage perimeter, we just need the first moment (cid:104) x m ( t ) (cid:105) = (cid:82) ∞ x m q t ( x m ) dx m , where q t ( x m ) denotes the probability den-sity of the of the maximum of the one dimensional com-ponent process. It is convenient to consider the cumula-tive distribution Q t ( x m ) i.e., the probability that the max-imum x -displacement stays below a given value x m up totime t . Then, q t ( x m ) = dQ t ( x m ) /dx m and (cid:104) x m ( t ) (cid:105) = (cid:82) ∞ [1 − Q t ( x m )] dx m . Since the process starts at the origin,its maximum x -displacement, for any time t , is necessarilynonnegative, i.e., x m ≥
0. We next write down a backwardFokker-Planck equation describing the evolution of Q t ( x m ) byconsidering the three mutually exclusive stochastic moves ina small time interval dt : starting at the origin at t = 0, thewalker during the subsequent interval [0 , dt ] dies with prob-ability γdt , infects another individual (i.e., branches) withprobability b dt = R γdt , or diffuses by a random displace-ment ∆ x = η x (0) dt with probability 1 − γ (1 + R ) dt . In thelast case, its new starting position is ∆ x for the subsequentevolution. Hence, for all x m ≥
0, one can write Q t + dt ( x m ) = γdt + R γdtQ t ( x m )+[1 − γ ( R + 1)] dt (cid:104) Q t ( x m − ∆ x ) (cid:105) , [8] where the expectation (cid:104)(cid:105) is taken with respect to the randomdisplacements ∆ x . The first term means that if the processdies right at the start, its maximum up to t is clearly 0 andhence is necessarily less than x m . The second term denotesthe fact that in case of branching the maximum of each branchstays below x m : since the branches are independent, one getsa square. The third term corresponds to diffusion. By using (cid:104) ∆ x (cid:105) = 0 and (cid:104) ∆ x (cid:105) = 2 Ddt and expanding Eq. [8] to thefirst order in dt and second order in ∆ x we obtain ∂∂t Q = D ∂ ∂x m Q − γ ( R + 1) Q + γR Q + γ [9] for x m ≥
0, satisfying the boundary conditions Q t (0) = 0 and Q t ( ∞ ) = 1, and the initial condition Q ( x m ) = Θ( x m ), whereΘ is the Heaviside step function. Hence from Eq. [5] (cid:104) L ( t ) (cid:105) = 2 π (cid:90) ∞ [1 − Q t ( x m )] dx m . [10] Equation [9] can be solved numerically for all t and all R ,which allows subsequently computing (cid:104) L ( t ) (cid:105) in Eq. [10] (de-tails and figures are provided in the Supplementary Material). The convex hull area.
To compute the average area in Eq. [7],we need to evaluate (cid:104) x m ( t ) (cid:105) as well as (cid:104) t m (cid:105) . Once the cumula-tive distribution Q t ( x m ) is known, the second moment (cid:104) x m ( t ) (cid:105) can be directly computed by integration, namely, (cid:104) x m ( t ) (cid:105) = (cid:82) ∞ dx m x m (1 − Q t ( x m )). To determine (cid:104) t m (cid:105) , we need to alsocompute the probability density p t ( t m ) of the random vari-able t m . Unfortunately, writing down a closed equation for p t ( t m ) is hardly feasible. Instead, we first define P t ( x m , t m )as the joint probability density that the maximum of the x component achieves the value x m at time t m , when the fullprocess is observed up to time t . Then, we derive a backwardevolution equation for P t ( x m , t m ) and then integrate out x m to derive the marginal density p t ( t m ) = (cid:82) ∞ P t ( x m , t m ) dx m .Following the same arguments as those used for Q t ( x m ) yieldsa backward equation for P t ( x m , t m ): P t + dt ( x m , t m ) = [1 − γ ( R + 1) dt ] (cid:104) P t ( x m − ∆ x, t m − dt ) (cid:105) +2 γR dtQ t ( x m ) P t ( x m , t m − dt ) . [11] The first term at the right hand side represents the contribu-tion from diffusion. The second term represents the contribu-tion from branching: we require that one of them attains themaximum x m at the time t m − dt , whereas the other staysbelow x m ( Q t ( x m ) being the probability that this conditionis satisfied). The factor 2 comes from the interchangeabilityof the particles. Developing Eq. [11] to leading order gives (cid:20) ∂∂t + ∂∂t m (cid:21) P t = (cid:20) D ∂ ∂x m − γ ( R + 1) + 2 γ R Q t (cid:21) P t . [12] This equation describes the time evolution of P t ( x m , t m ) inthe region x m ≥ ≤ t m ≤ t . It starts from the ini-tial condition P ( x m , t m ) = δ ( x m ) δ ( t m ) (since the processbegins with a single infected with x component located at x = 0, it implies that at t = 0 the maximum x m = 0 and also t m = 0). For any t > x m >
0, we have the condition P t ( x m ,
0) = 0. We need to also specify the boundary condi-tions at x m = 0 and x m → ∞ , which read (i) P t ( ∞ , t m ) = 0(since for finite t the maximum is necessarily finite) and (ii) P t (0 , t m ) = δ ( t m ) q t ( x m ) | x m =0 . The latter condition comesfrom the fact that, if x m = 0, this corresponds to the eventthat the x component of the entire process, starting at 0 ini-tially, stays below 0 in the time interval [0 , t ], which happenswith probability q t ( x m ) | x m =0 : consequently, t m must neces-sarily be 0. Furthermore, by integrating P t ( x m , t m ) with re-spect to t m we recover the marginal density q t ( x m ).The numerical integration of the full Eq. [12] would berather cumbersome. Fortunately, we do not need this. Sincewe are only interested in (cid:104) t m (cid:105) , it is convenient to introduce T t ( x m ) = (cid:90) t t m P t ( x m , t m ) dt m , [13] from which the average follows as (cid:104) t m (cid:105) = (cid:82) dx m T t ( x m ). Mul-tiplying Eq. [12] by t m and integrating by parts we get ∂∂t T t − q t ( x m ) = (cid:20) D ∂ ∂x m + 2 γR Q t − γ ( R + 1) (cid:21) T t , [14] with the initial condition T ( x m ) = 0, and the boundary con-ditions T t (0) = 0 and T t ( ∞ ) = 0. Eq. [14] can be integratednumerically, together with Eq. [9] (details are provided in theSupplementary Material), and the behavior of (cid:104) A ( t ) (cid:105) = π (cid:90) ∞ dx m [2 x m (1 − Q t ( x m )) − T t ( x m )] [15] as a function of time is illustrated in Fig. 2. The critical regime.
We now focus on the critical regime R = 1. We begin with the average perimeter: for R = 1,Eq. [9] admits a stationary solution as t → ∞ , which canbe obtained by setting ∂Q/∂t = 0 and solving the resultingdifferential equation. In fact, this stationary solution was al-ready known in the context of the genetic propagation of amutant allele [22]. Taking the derivative of this solution with Footline Author PNAS
Issue Date
Volume Issue Number espect to x m , we get the stationary probability density of themaximum x m q ∞ ( x m ) = ∂ x m Q ∞ ( x m ) = 2 (cid:112) γ D (cid:0) (cid:112) γ D x m (cid:1) . [16] The average is (cid:104) x m (cid:105) = (cid:82) ∞ x m q ∞ ( x m ) dx m = (cid:112) D/γ , whichyields then Eq. [1] for the average perimeter of the convexhull at late times.To compute the average area in Eq. [7], we need toalso evaluate the second moment (cid:104) x m ( t ) (cid:105) , which diverges as t → ∞ , due to the power-law tail of the stationary probabil-ity density q ∞ ( x m ) ∝ x − m for large x m . Hence, we need toconsider large but finite t . In this case, the time dependentprobability density q t ( x m ) displays a scaling form which canbe conveniently written as q t ( x m ) (cid:39) q ∞ ( x m ) f (cid:18) x m √ Dt (cid:19) , [17] where f ( z ) is a rapidly decaying function with f ( z (cid:28) (cid:39) f ( z (cid:29) (cid:39)
0. Using the scaling form of Eq. [17] andEq. [9] one can derive a differential equation for f ( z ). But itturns out that we do not really need the solution of f ( z ).From Eq. [17] we see that the asymptotic power-law decayof q t ( x m ) for large x m has a cut-off around x ∗ m ∼ √ Dt and f ( z ) is the cut-off function. The second moment at finite butlarge times t is given by (cid:104) x m ( t ) (cid:105) = (cid:82) ∞ x m q t ( x m ) dx m . Sub-stituting the scaling form and cutting off the integral over x m at x ∗ m = c √ t (where the constant c depends on the preciseform of f ( z )) we get, to leading order for large t , (cid:104) x m ( t ) (cid:105) (cid:39) (cid:90) x ∗ m x m q ∞ ( x m ) dx m (cid:39) Dγ ln t . [18] Thus, interestingly the leading order result is universal, i.e,independent of the details of the cut-off function f ( z ) (the c -dependence is only in the subleading term). To complete thecharacterization of (cid:104) A ( t ) (cid:105) in Eq. [7], we still need to determine (cid:104) t m (cid:105) : in the Supplementary Material we explicitly determinethe stationary solution P ∞ ( x m , t m ) for R = 1. By followingthe same arguments as for (cid:104) x m ( t ) (cid:105) , we show that (cid:104) t m (cid:105) (cid:39) γ ln t [19] for large t , which leads again to a logarithmic divergence intime. Finally, substituting Eqs. [18] and [19] in Eq. [7] givesthe result announced in Eq. [2].A deeper understanding of the statistical properties of theprocess would demand knowing the full distribution Prob( L, t )and Prob(
A, t ) of the perimeter and area. These seem ratherhard to compute, but one can obtain the asymptotic tailsof the distributions by resorting to scaling arguments. Fol-lowing the lines of Cauchy’s formula (see the SupplementaryMaterial), it is reasonable to assume that for each sample theperimeter scales as L ( t ) ∼ x m ( t ). We have seen that the distri-bution of x m ( t ) has a power-law tail for large t : q ∞ ( x m ) ∼ x − m for large x m . Then, assuming the scaling L ( t ) ∼ x m ( t ) andusing Prob( L, t → ∞ ) dL ∼ q ∞ ( x m ) dx m , it follows that atlate times the perimeter distribution also has a power-lawtail: Prob( L, t → ∞ ) ∼ L − for large L . Similarly, usingthe Cauchy formula for the area, we can reasonably assumethat for each sample A ( t ) ∼ x m ( t ) in the scaling regime. Onceagain, using Prob( A, t → ∞ ) dA = q ∞ ( x m ) dx m , we find thatthe area distribution also converges, for large t , to a stationarydistribution with a power-law tail: Prob( A, t → ∞ ) ∼ A − for large A . Moreover, the logarithmic divergence of the meanarea calls for a precise ansatz on the tail of the area distribu-tion, namely,Prob( A, t ) −−−→ A (cid:29) πD γ A − h (cid:18) ADt (cid:19) , [20] where the scaling function h ( z ) satisfies the conditions h ( z (cid:28)
1) = 1, and h ( z (cid:29) (cid:39)
0. It is not difficult to verify that thisis the only scaling compatible with Eq. [2]. These two resultsare consistent with the fact that for each sample typically A ( t ) ∼ L ( t ) at late times in the scaling regime. Our scalingpredictions are in agreement with our Monte Carlo simula-tions (see Fig. 2). The power-law behavior of Prob( A, t ) im-plies that the average area is not representative of the typicalbehavior of the epidemic area, which is actually dominated byfluctuations and rare events, with likelihood given by Eq. [20].
The supercritical regime.
When R >
1, it is convenient torewrite Eq. [9] in terms of W ( x m , t ) = 1 − Q ( x m , t ): ∂∂t W = D ∂ ∂x m W + γ ( R − W − γR W [21] starting from the initial condition W ( x m ,
0) = 0 for all x m > (cid:104) L ( t ) (cid:105) = 2 π (cid:82) ∞ W ( x m , t ) dx m isjust the area under the curve W ( x m , t ) vs. x m , up to a factor2 π . As t → ∞ , the system approaches a stationary state for all R ≥
1, which can be obtained by setting ∂ t W = 0 in Eq. [21].For R > W ( x m , ∞ ) approaches theconstant 1 − /R exponentially fast as x m → ∞ , namely, W ( x m , ∞ ) − R − → exp[ − x m /ξ ], with a characteristiclength scale ξ = (cid:112) D/γ ( R − t , W ( x m , t ) as a function of x m has a two-step form: it firstdecreases from 1 to its asymptotic stationary value 1 − /R over the length scale ξ , and then decreases rather sharply from1 − /R to 0. The frontier between the stationary asymptoticvalue 1 − /R (stable) and 0 (unstable) moves forward withtime at constant velocity, thus creating a travelling front at theright end, which separates the stationary value 1 − /R to theleft of the front and 0 to the right. This front advances witha constant velocity v ∗ that can be estimated using the stan-dard velocity selection principle [27, 28, 29]. Near the frontwhere the nonlinear term is negligible, the equation admitsa travelling front solution: W ( x m , t ) ∼ exp[ − λ ( x m − v t )],with a one parameter family of possible velocities v ( λ ) = Dλ + γ ( R − /λ , parametrized by λ . This dispersion rela-tion v ( λ ) has a minimum at λ = λ ∗ = (cid:112) γ ( R − /D , where v ∗ = v ( λ ∗ ) = 2 (cid:112) Dγ ( R − v ∗ . The width of the front remains of ∼ O (1) at large t . Thus,due to this sharpness of the front, to leading order for large t one can approximate W ( x m , t ) (cid:39) (1 − /R )Θ( v ∗ t − x m )near the front. Hence, to leading order for large t one gets (cid:104) x m ( t ) (cid:105) (cid:39) (1 − /R ) v ∗ t and (cid:104) x m (cid:105) (cid:39) (1 − /R ) ( v ∗ t ) . Theformer gives, from Eq. [5], the result announced in Eq. [3].For the mean area in Eq. [7], the term (cid:104) x m (cid:105) ∼ t for large t dominates over (cid:104) t m (cid:105) ∼ t (which can be neglected), and we getthe result announced in Eq. [3]. Conclusions
In this paper, we have developed a general procedure for as-sessing the time evolution of the convex hull associated to theoutbreak of an epidemic. We find it extremely appealing thatone can successfully use mathematical formulae (Cauchy’s) rom two dimensional integral geometry to describe the spa-tial extent of an epidemic outbreak in relatively realistic sit-uations. Admittedly, there are many assumptions in this epi-demic model that are not quite realistic. For instance, we haveignored the fluctuations of the susceptible populations duringthe early stages of the epidemic: this hypotheses clearly breaksdown at later times, when depletion effects begin to appear,due to the epidemic invading a thermodynamical fraction ofthe total population. In addition, we have assumed that thesusceptibles are homogeneously distributed in space, which isnot the case in reality. Nonetheless, it must be noticed thatin practical applications whenever strong heterogeneities ap-pear, such as mountains, deserts or oceans, one can split theanalysis of the evolving phenomena by conveniently resortingto several distinct convex hulls, one for each separate region.For analogous reasons, the convex hull approach would not besuitable to characterize birth-death processes with long rangedisplacements, such as for instance branching L´evy flights.The model discussed in this paper based on branchingBrownian motion is amenable to exact results. More gener-ally, realistic models could be taken into account by resortingto cumbersome Monte Carlo simulations. The approach pro-posed in this paper paves the way for assessing the spatialdynamics of the epidemic by more conveniently solving twocoupled nonlinear equations, under the assumption that theunderlying process be rotationally invariant.We conclude with an additional remark. In our computa-tions of the mean perimeter and area, we have averaged overall realizations of the epidemics up to time t , including thosewhich are already extinct at time t . It would also be interest-ing to consider averages only over the ensemble of epidemicsthat are still active at time t . In this case we expect differ-ent scaling laws for the mean perimeter and the mean area ofthe convex hull. In particular, in the critical case, we believethat the behavior would be much closer to that of a regularBrownian motion. Appendix: Supplementary MaterialsCauchy’s formula.
The problem of determining the perime-ter and the area of the convex hull of any two dimensionalstochastic process [ x ( τ ) , y ( τ )] with 0 ≤ τ ≤ t can be mappedto that of computing the statistics of the maximum and thetime of occurrence of the maximum of the one dimensionalcomponent process x ( τ ) [23, 24]. This is achieved by resort-ing to a formula due to Cauchy, which applies to any closedconvex curve C .A sketch of the method is illustrated in Fig. 5. Choosethe coordinates system such that the origin is inside the curve C and take a given direction θ . For fixed θ , consider a stickperpendicular to this direction and imagine bringing the stickfrom infinity and stop upon first touching the curve C . Atthis point, the distance M ( θ ) of the stick from the origin iscalled the support function in the direction θ . Intuitively, thesupport function measures how close can one get to the curve C in the direction θ , coming from infinity. Once the supportfunction M ( θ ) is known, then Cauchy’s formulas [30] give theperimeter L and the area A enclosed by C , namely L = (cid:82) π M ( θ ) dθA = (cid:82) π (cid:2) M ( θ ) − ( M (cid:48) ( θ )) (cid:3) dθ, [22] where M (cid:48) ( θ ) = dM/dθ . For example, for a circle ofradius R = r , M ( θ ) = r , and one recovers the stan-dard formulae: L = 2 πr and A = πr . When C isthe convex hull of associated with the process at time t , we first need to compute its associated support function M ( θ ). A crucial point is to realize that actually M ( θ ) =max ≤ τ ≤ t [ x ( τ ) cos( θ ) + y ( τ ) sin( θ )] [23, 24]. Furthermore, ifthe process is rotationally invariant any average is indepen-dent of the angle θ . Hence for the average perimeter wecan simply set θ = 0 and write (cid:104) L ( t ) (cid:105) = 2 π (cid:104) M (0) (cid:105) , wherebrackets denote the ensemble average over realizations. Sim-ilarly for the average area, (cid:104) A ( t ) (cid:105) = π (cid:2) (cid:104) M (0) (cid:105) − (cid:104) M (cid:48) (0) (cid:105) (cid:3) .Clearly, M (0) = max ≤ τ ≤ t [ x ( τ )] is then the maximum of theone dimensional component process x ( τ ) for τ ∈ [0 , t ]. As-suming that x ( τ ) takes its maximum value x ( t m ) at time τ = t m (see Fig. 4). Then, M (0) = x ( t m ) = x m ( t ), and M (cid:48) (0) = y ( t m ) [31]. Now, by taking the average over Cauchy’sformulas, and using isotropy, we simply have Eqs. [5] and [6]for the mean perimeter and the mean area of the convex hull C at time t . Note that this argument is very general and is appli-cable to any rotationally invariant two dimensional stochasticprocess. Since the branching Brownian motion with death isrotationally invariant we can use these formulae. Numerical methods.
Numerical integration.
Equations [9] and[14] have been integrated numerically by finite differencesin the following way. Time has been discretized by setting t = ndt , and space by setting x = idx , where dt and dx aresmall constants. For the sake of simplicity, here we considerthe case R = 1. We thus have Q n +1 ( i ) == Q n ( i ) + γ dt [1 − Q n ( i )] + D dt ( dx ) [ Q n ( i + 1) − Q n ( i ) + Q n ( i − [23] and T n +1 ( i ) == T n ( i ) + 2 γ dt T n ( i ) [ Q n ( i ) −
1] +
D dt ( dx ) [ T n ( i + 1) − T n ( i ) + T n ( i − dtdx [ T n ( i ) − T n ( i − . [24] As for the initial conditions, Q (0) = 0 and Q ( i >
0) = 1,and T ( i ) = 0 ∀ i . The boundary conditions at the originare Q n (0) = 0 and T n (0) = 0. In order to implement theboundary condition at infinity, we impose Q n ( i max ) = 1 and T n ( i max ) = 0 ∀ n , where the large value i max is chosen so that T n ( i max ) − T n ( i max − < − . We have verified that numeri-cal results do not change when passing to the tighter condition T n ( i max ) − T n ( i max − < − .Once Q n ( i ) and T n ( i ) are known, we use Eqs. [10] and [15]to determine the average perimeter and area, respectively. Monte Carlo simulations.
The results of numerical in-tegrations have been confirmed by running extensive MonteCarlo simulations. Branching Brownian motion with deathhas been simulated by discretizing time with a small dt : ineach interval dt , with probability bdt the walker branches andthe current walker coordinates are copied to create a new ini-tial point, which is then stored for being simulated in thenext dt ; with probability γdt the walker dies and is removed;with probability 1 − ( b + γ ) dt the walker diffuses: the x and y displacements are sampled from Gaussian densities of zeromean and standard deviation √ Ddt and the particle posi-tion is updated. The positions of all the random walkers arerecorded as a function of time and the corresponding convexhull is then computed by resorting to the algorithm proposedin [32].
Footline Author PNAS
Issue Date
Volume Issue Number erimeter statistics. In order the complete the analysisof the convex hull statistics, in Fig 6 and Fig 7 we show theresults for the perimeter.
Analysis of t m . In the critical case R = 1, the stationaryjoint probability density P ∞ ( x m , t m ) satisfies (upon setting ∂P t /∂t = 0 in Eq. [12]) ∂∂t m P ∞ ( x m , t m ) == (cid:34) D ∂ ∂x m − γ (cid:2) (cid:112) γ D x m (cid:3) (cid:35) P ∞ ( x m , t m ) . [25] For any x m >
0, we have the condition P ∞ ( x m ,
0) = 0. Theboundary conditions for Eq. (25) are P ∞ ( x m → ∞ , t m ) = 0and P ∞ (0 , t m ) = q ∞ (0) δ ( t m ) = 2 (cid:112) γ/ (6 D ) δ ( t m ). We firsttake the Laplace transform of (25), namely,˜ P ∞ ( x m , s ) = (cid:90) ∞ e − st m P ∞ ( x m , t m ) dt m . [26] This gives for all x m > Ds ∂ ∂x m ˜ P ∞ ( x m , s ) = sD ( (cid:113) Dγ + x m ) ˜ P ∞ ( x m , s ) , [27] where we have used the condition P ∞ ( x m ,
0) = 0 for any x m >
0. This second order differential equation satisfiestwo boundary conditions: ˜ P ∞ ( ∞ , s ) = 0 and ˜ P ∞ (0 , s ) =2 (cid:112) γ/ (6 D ). The latter condition is obtained by Laplace trans-forming P ∞ (0 , t m ) = 2 (cid:112) γ/ (6 D ) δ ( t m ). By setting z = (cid:18)(cid:114) Dγ + x m (cid:19) (cid:114) sD , [28] we rewrite the equation as ∂ ∂z ˜ P ∞ − ˜ P ∞ − z ˜ P ∞ = 0 . [29] Upon making the transformation ˜ P ∞ ( z ) = √ z F ( z ), the func-tion F ( z ) then satisfies the Bessel differential equation d dz F ( z ) + 1 z ddz F ( z ) − (cid:20) z (cid:21) F ( z ) = 0 . [30] The general solution of this differential equation can be ex-pressed as a linear combination of two independent solutions: F ( z ) = A I / ( z )+ B K / ( z ) where I ν ( z ) and K ν ( z ) are mod-ified Bessel functions. Since, I ν ( z ) ∼ e z for large z , it isclear that to satisfy the boundary condition ˜ P ∞ ( ∞ , s ) = 0(which mean F ( z → ∞ ) = 0), we need to choose A = 0.Hence we are left with F ( z ) = BK / ( z ), where the con-stant B is determined from the second boundary condition ˜ P ∞ (0 , s ) = 2 (cid:112) γ/ (6 D ). By reverting to the variable x m , wefinally get˜ P ∞ ( x m , s ) = 2 (cid:114) γ D (cid:114) γ D x m K / (cid:104)(cid:16)(cid:113) Dγ + x m (cid:17) (cid:112) sD (cid:105) K / (cid:104)(cid:113) sγ (cid:105) . [31] Now, we are interested in determining the Laplace trans-form of the marginal density ˜ p ∞ ( s ) = (cid:82) ∞ e − s t m p ∞ ( t m ) dt m where p ∞ ( t m ) = (cid:82) ∞ P ∞ ( x m , t m ) dx m . Taking Laplace trans-form of this last relation with respect to t m gives ˜ p ∞ ( s ) = (cid:82) ∞ ˜ P ∞ ( x m , s ) dx m . Once we know ˜ p ∞ ( s ), we can invert it toobtain p ∞ ( t m ). Since we are interested only in the asymptotictail of p ∞ ( t m ), it suffices to investigate the small s behaviorof ˜ p ∞ ( s ). Integrating Eq. (31) over x m and taking the s → p ∞ ( s ) = 1 + 35 γ s ln( s ) + · · · . [32] We further note that (cid:90) ∞ e − st m t m p ∞ ( t m ) dt m = d ds ˜ p ∞ ( s ) (cid:39) γs , [33] which can then be inverted to give the following asymptoticbehavior for large t m p ∞ ( t m ) (cid:39) γt m . [34] Analogously as for (cid:104) x m (cid:105) , the moment (cid:104) t m (cid:105) → ∞ , due to thepower-law tail p ∞ ( t m ) ∝ t − m . Hence we need to compute (cid:104) t m (cid:105) for large but finite t : in this case, the time-dependent solutiondisplays a scaling behavior p t ( t m ) (cid:39) p ∞ ( t m ) g (cid:18) t m t (cid:19) , [35] where the scaling function g ( z ) satisfies the conditions g ( z (cid:28) (cid:39) g ( z (cid:29)
1) = 0. Much like in Eq. [17] for themarginal density q t ( x m ), we have a power-law tail of p t ( t m )for large t m that has a cut-off at a scale t ∗ m ∼ t , and g ( z ) isthe cut-off function. As in the case of x m , we do not need theprecise form of g ( z ) to compute the leading term of the firstmoment (cid:104) t m (cid:105) = (cid:82) ∞ p t ( t m ) t m dt m for large t . Cutting off theintegral at t ∗ m = c t (where c depends on the precise form of g ( z )) and performing the integration gives (cid:104) t m (cid:105) (cid:39) (cid:90) t t m p ∞ ( t m ) dt m (cid:39) γ ln t, [36] which is precisely the result announced in Eq. [19]. ACKNOWLEDGMENTS.
S.N.M. acknowledges support from the ANR grant 2011-BS04-013-01 WALK-MAT. S.N.M and A.R. acknowledge support from the Indo-French Centre for thePromotion of Advanced Research under Project 4604-3.
1. Bailey N T J (1975) The Mathematical Theory of Infectious Diseases and its Appli-cations, Griffin, London2. McKendrick A G (1925) Applications of mathematics to medical problems. Kapil Pro-ceedings of the Edinburgh Mathematical Society 44:1-343. P. Whittle (1955) The outcome of a stochastic epidemic - a note on Baileys paper.Biometrika 42:1161224. Murray J D (1989) Mathematical Biology, Springer-Verlag, Berlin5. Bartlett M S (1960) Stochastic Population Models in Ecology and Epidemiology 6. Andersson H, Britton T (2000) Stochastic Epidemic Models and their Statistical Anal-ysis Lecture Notes in Statistics 151, Springer-Verlag, New York7. Antal T, Krapivsky P L (2012) Outbreak size distributions in epidemics with multiplestages. arXiv:1204.42148. Anderson R, May R (1991) Infectious Diseases: Dynamics and Control, Oxford Uni-versity Press, Oxford9. Antia R, Regoes R R, Koella J C, Bergstrom C T (2003) The role of evolution in theemergence of infectious diseases. Nature 426:658-661
0. Kendall D G (1956) Deterministic and stochastic epidemics in closed populations. InProc. 3rd Berkeley Symp. Math. Statist. Prob. 4:149-16511. Bartlett M S (1956) An introduction to stochastic processes, Cambridge UniversityPress12. Elliott P, Wakefield J C, Best N G, Briggs D J (2000) Spatial Epidemiology: Methodsand Applications, Oxford University Press13. Radcliffe J (1976) The Convergence of a Position-Dependent Branching Process Usedas an Approximation to a Model Describing the Spread of an Epidemic. Journal ofApplied Probability 13:338-34414. Wang J S (1980) The Convergence of a Branching Brownian Motion Used as a ModelDescribing the Spread of an Epidemic. Journal of Applied Probability 17:301-31215. Riley S, et al (2007) Large-Scale Spatial-Transmission Models of Infectious Disease.Science 316:1298-130116. Fraser C, et al (2009) Pandemic Potential of a Strain of Influenza A (H1N1): EarlyFindings Science 324:1557-156117. Colizza V, Barrat A, Barth´elemy M, Vespignani A (2006) The role of the airline trans-portation network in the prediction and predictability of global epidemics. Proc. Natl.Acad. Sci. USA 103:201518. Brockmann D, Hufnagel L, Geisel T (2006) The scaling laws of human travel. Nature439:462-46519. Cowled B, Ward M P, Hamilton S, Garner G (2009) The equine influenza epidemic inAustralia: Spatial and temporal descriptive analyses of a large propagating epidemic.Preventive Veterinary Medicine 92:607020. Worton B J (1995) A convex hull-based estimator of home-range size. Biometrics 51:1206-121521. Giuggioli L, Abramson G, Kenkre V M, Parmenter R R, Yates T L (2006) Theoryof home range estimation from displacement measurements of animal populations. J.Theor. Biol. 240: 126-135. 22. Sawyer S and Fleischman J (1979) Maximum geographic range of a mutant alleleconsidered as a subtype of a Brownian branching random field. Proc Natl Acad SciUSA 76:872875.23. Randon-Furling J, Majumdar S N, Comtet A (2009) Convex Hull of N planar BrownianMotions: Application to Ecology. Phys. Rev. Lett. 103:14060224. Majumdar S N, Comtet A, Randon-Furling J (2010) Random Convex Hulls and Ex-treme Value Statistics. J. Stat. Phys. 138:955-100925. Reymbaut A, Majumdar S N, Rosso A (2011) The Convex Hull for a Random Accel-eration Process in Two Dimensions. J. Phys. A-Math. & Theor. 44:41500126. Cauchy A (1850) Mem. Acad. Sci. Inst. Fr. 22: 3 ; see also the book by L. A. Santal´o,Integral Geometry and Geometric Probability (Addison-Wesley, Reading, MA, 1976)27. van Saarloos W (2003) Front propagation into unstable states. Phys. Rep. 386: 29-222.28. Brunet E, Derrida B (2009) Statistics at the tip of a branching random walk and thedelay of traveling waves Europhys. Lett. 87:60010.29. Majumdar S N, Krapivsky P L (2003) Extreme value statistics and traveling fronts:various applications. Physica A 318:161-170.30. A. Cauchy (1850), Mem. Acad. Sci. Inst. Fr. 22, 3; see also the book by L. A. Santal´o,Integral Geometry and Geometric Probability (Addison-Wesley, Reading, MA, 1976).31. Actually, tm implicitly depends on θ , hence formally M (cid:48) ( θ ) = − x ( tm ) sin( θ ) + y ( tm ) cos( θ ) + dtmdθ dzθ ( t ) dt (cid:12)(cid:12)(cid:12) t = tm . Nonetheless, since zθ ( t ) is maximum at t = tm , by definition dzθ ( t ) /dt | t = tm = 0 .32. Module: Finding the convex hull of a set of 2D points, by Gehrman D C, PythonCookbook, edited by Martelli A and Ascher D Footline Author PNAS
Issue Date
Volume Issue Number O O
Day 1 Day 2 Day 3 xy xy xy
Fig. 1.
The snapshots of the trajectories of an assembly of infected individuals at the epidemics outbreak at three different times (schematic), starting from a single infectedat the origin O at time t = 0 . Individuals that are still infected at a given time t are displayed as red dots, while those already recovered are shown as black dots. The convexhull enclosing the trajectories (shown as a dashed line) is a measure of geographical area covered by the spreading epidemic. As the epidemic grows in space, the associatedconvex hull also grows in time. -9 -8 -7 -6 -5 -4 -3 -2 P r ob ( A ,t ) A 0 500 1000 1500 2000 2500 3000 3500 10 < A ( t ) > t R = . R = R = . R = . R = . Fig. 2.
Left. The average area (cid:104) A ( t ) (cid:105) of the convex hull as a function of the observation time. For the parameter values, we have chosen D = 1 / and b = R γ = 0 . .We considered five different values of R . We have obtained these results by two different methods: (i) via the numerical integration of Eqs. [9] and [14] and using Eq. [15].These results are displayed as solid lines. (ii) by Monte Carlo simulations of the two-dimensional branching Brownian motion with death with the same parameters, averagedover samples. Monte Carlo are displayed as symbols. The dashed lines represent the asymptotic limits as given in Eq. [2] for the critical case R = 1 . Further details ofthe numerical simulations are provided in the Supplementary Material. Right. Distribution of the area of the convex hull for the critical case R = 1 , with γ = 0 . and D = 1 / , as obtained by Monte Carlo simulations with · realizations. The dashed line corresponds to the power-law (24 πD/ γ ) A − as predicted by Eq [20]. t = t = t = t = ∞
0 10 30 40x t = t = t = ∞ < A ( t ) > t R = . R = . R = . R = Fig. 3.
Left. The time behavior of the average area in the supercritical regime for different values of R > . Dashed lines represent the asymptotic scaling as in Eq.[4]. The red curve corresponds to the critical regime. Right. The behavior of W t ( x ) = 1 − Q t ( x ) for R = 1 . at different times, as in Eq. [21]. When t → ∞ , W t ( x ) → − R − , and for large but finite times the travelling front behavior is clearly visible. The inset displays the exponential convergence of W t ( x ) to the asymptoticlimit. The dashed line represents ξ = (cid:112) D/γ ( R − . xy x t t yx m y( t m ) x m y( t m ) t m t m Fig. 4.
Left. A branching random walk composed of five individuals. At time t = 0 , a single infected is at the origin O , and starts diffusing (blue line). At later times, thisindividual branches and gives rise to other infected individuals. Among these, the red path reaches the maximum x m along the x component up to the final time t . Infectedindividuals at a given time t are displayed as red dots, whereas recovered as black dots. Center. The displacement along the x direction as a function of time. The red pathreaches the global maximum x m at time t m . Right. The displacement along the y direction as a function of time. When the red path reaches the global maximum x m at time t m , its y coordinate attains the value y ( t m ) . A crucial observation is that the y component of the trajectory connecting O to the red path is a regular Brownianmotion. This is not the case for the x component, which is constrained to reach the global maximum of the branching process. O y xM( θ ) C θ Fig. 5.
Cauchy’s construction of the two-dimensional convex hull, with support function M ( θ ) representing the distance along the direction θ .Footline Author PNAS Issue Date
Volume Issue Number < L ( t ) > t 10 -6 -5 -4 -3 -2 P r ob ( L ,t ) L R = . R = R = . R = . R = . Fig. 6.
Left. The average perimeter (cid:104) L ( t ) (cid:105) of the convex hull as a function of the observation time. For the parameter values, we have chosen D = 1 / and b = R γ = 0 . . We considered five different values of R . We have obtained these results by two different methods: (i) via the numerical integration of Eq. [9] and usingEq. [10] (with the choices dt = 0 . and dx = 0 . ). These results are displayed as solid lines. (ii) by Monte Carlo simulations of the two-dimensional branchingBrownian motion with death with the same parameters and with the choice of the Monte Carlo time step dt = 0 . with the results averaged over samples. MonteCarlo are displayed as symbols. The dashed lines represent the asymptotic limits as given in Eq. [1] for the critical case R = 1 . Right. Distribution of the perimeter of theconvex hull for the critical case R = 1 , with γ = 0 . and D = 1 / , as obtained by Monte Carlo simulations with time step dt = 1 and t = 4 · . The number ofrealizations is · . The dashed line of the left panel corresponds to the power-law L − (up to an arbitrary prefactor). < L ( t ) > t R = . R = . R = . R = Fig. 7.
The time behavior of the average perimeter in the supercritical regime for different values of R > . Dashed lines represent the asymptotic scaling as in Eq. [3].The red curve corresponds to the critical regime.10