From individual-based epidemic models to McKendrick-von Foerster PDEs: A guide to modeling and inferring COVID-19 dynamics
Félix Foutel-Rodier, François Blanquart, Philibert Courau, Peter Czuppon, Jean-Jil Duchamps, Jasmine Gamblin, Élise Kerdoncuff, Rob Kulathinal, Léo Régnier, Laura Vuduc, Amaury Lambert, Emmanuel Schertzer
FFrom individual-based epidemic models toMcKendrick-von Foerster PDEs: A guide to modeling and inferring COVID-19 dynamics
Félix Foutel-Rodier (cid:63), , François Blanquart , Philibert Courau Peter Czuppon , Jean-Jil Duchamps , Jasmine Gamblin Élise Kerdoncuff , Rob Kulathinal , Léo Régnier Laura Vuduc Amaury Lambert (cid:63)(cid:63), , Emmanuel Schertzer (cid:63)(cid:63), , July 21, 2020 (cid:63) first author (cid:63)(cid:63) co-last authors SMILE Group, Center for Interdisciplinary Research in Biology (CIRB), Collège deFrance, CNRS UMR 7241, INSERM U 1050, PSL Research University, Paris, France Laboratoire de Probabilités, Statistique & Modélisation (LPSM), Sorbonne Université,CNRS UMR 8001, Université de Paris, Paris, France Infection, Antimicrobials, Modeling, Evolution (IAME), Université de Paris, INSERMUMR 1137, Paris, France Institute of Ecology and Environmental Sciences of Paris (IEES), Sorbonne Université,UPEC, CNRS UMR 7618, IRD, INRA, Paris, France Institut de Systématique, Biodiversité, Évolution (ISYEB), Muséum National d’HistoireNaturelle (MNHN), CNRS UMR 7205, Paris, France Department of Biology, Temple University, Philadelphia, PA, USA1 a r X i v : . [ q - b i o . P E ] J u l bstract We present a unifying, tractable approach for studying the spread of virusescausing complex diseases that require to be modeled using a large number of types(e.g., infective stage, clinical state, risk factor class). We show that recording eachinfected individual’s infection age, i.e., the time elapsed since infection,1. The age distribution n ( t, a ) of the population at time t can be described bymeans of a first-order, one-dimensional partial differential equation (PDE)known as the McKendrick-von Foerster equation.2. The frequency of type i at time t is simply obtained by integrating the proba-bility p ( a, i ) of being in state i at age a against the age distribution n ( t, a ) .The advantage of this approach is three-fold. First, regardless of the number of types,macroscopic observables (e.g., incidence or prevalence of each type) only rely on aone-dimensional PDE “decorated” with types. This representation induces a simplemethodology based on the McKendrick-von Foerster PDE with Poisson sampling toinfer and forecast the epidemic. We illustrate this technique using a French datafrom the COVID-19 epidemic.Second, our approach generalizes and simplifies standard compartmental modelsusing high-dimensional systems of ordinary differential equations (ODEs) to accountfor disease complexity. We show that such models can always be rewritten in ourframework, thus, providing a low-dimensional yet equivalent representation of thesecomplex models.Third, beyond the simplicity of the approach, we show that our population modelnaturally appears as a universal scaling limit of a large class of fully stochasticindividual-based epidemic models, where the initial condition of the PDE emergesas the limiting age structure of an exponentially growing population starting froma single individual. Contents Inference 25
A Maximum likelihood method 37B Best fitting prevalence curves under admission model 37
The transmission of pathogens between species is a global concern [1, 2]. As such zoonoticepisodes are expected to become increasingly common in humans, it is critical to developanalytic tools that can quickly transform epidemiological observations into informed publicpolicy in order to mitigate and control epidemics.A novel coronavirus, SARS-CoV-2, has recently crossed the species barrier into humansand, within months, has rapidly spread to all corners of our planet [3]. The sheer scaleof this pandemic has overburdened our medical infrastructure, caused fatalities estimatedwell into the hundreds of thousands, and shut down entire economies. Remarkably, therapid spread of COVID-19 and its consequences can be attributed to the unique life cycleof a 30,000 base pair single-stranded virus. SARS-CoV-2 is an airborne pathogen trans-mitted by both symptomatic and asymptomatic carriers in close proximity to non-infectedindividuals. Milder COVID-19 symptoms include a dry cough, fever, and/or shortnessof breath while more serious cases include respiratory failure and eventual death. Withmillions of infections and hundreds of thousands of documented deaths and recoveries, theCOVID-19 pandemic is providing a wealth of independent estimates of important clinicalcharacteristics that can help predict health outcomes specific for a country or region.It quickly became understood that accurate descriptions of the life cycle of this diseaseneeded to distinguish between several stages of the disease, referred to as compartments,depending on whether an infected individual is infectious or not, symptomatic or not,hospitalized, etc . However it remains unclear to what extent making precise predictionsof the dynamics of such a complex disease requires to have a precise knowledge of clinicalfeatures such as incubation period, generation time, and duration times between infection,symptom establishment, hospitalization, recovery and death, to know how these durationscorrelate and what are the exact probabilities of transition between stages.In this work, we consider a fully stochastic, generic epidemiological model with anarbitrary number of compartments, that encompasses life cycles of most complex diseasesand that of COVID-19 in particular. We show how structuring the infected populationby its infection age, i.e., time elapsed since infection, allows us to decouple dependenciesbetween stages and to time. More specifically, when the population size is large enough,the joint evolution of all compartment sizes can be described by means of a linear, first-order partial differential equation (PDE) known as the McKendrick-von Foerster equationdescribing the number n ( t, a ) of infecteds of (infection) age a at time t . The boundarycondition at age 0 is driven by the infection rate from infecteds of age a averaged over all3ife cycles and the number of individuals of age a in compartment i at time t is obtainedby thinning n ( t, a ) by a factor p ( a, i ) which is the probability of being in compartment i conditional on having age a , averaged over all life cycles.In the case of COVID-19, we display a simple procedure to infer these parameters,some from the biological literature and most from time series of numbers of severe cases,hospitalized cases, discharged patients and deaths that can be applied easily to any re-gional or national dataset. We also allow for time inhomogeneity in the infection rate toaccount for temporary mitigation measures such as lockdowns or social distancing. Weapply this procedure to French COVID-19 data from March to May 2020 and estimatevarious parameters of interest including R in different phases of the epidemic (before,during, and after lockdown) and biological parameter values that we compare to empiricalestimates. We consider a population model of the SIR fashion where each individual goes throughsuccessive stages, starting from stage S (susceptible) and ending in one of two states: R (recovered) or D (dead). Depending on disease complexity, the number of stages inthis life cycle can vary. In the SARS-CoV-2 example, typical intermediate stages are A (asymptomatic) or P (presymptomatic), I (mild case) or C (severe case), H (hospitalized), U (intensive care unit). These stages are sometimes called compartments, types, classes,stages or simply states. See Figure 4 for an illustration.We assume that S individuals are always in excess (branching assumption) and thateach individual infects new S individuals at successive times of a random point process,one at a time. We further assume that upon infection, an S individual immediatelychanges state and never returns to state S (ruling out multiple infections in particular).More formally, • The set of all possible states is denoted S and a stochastic process ( X ( a ); a ≥ with values in S gives the state of a typical individual of age a , where here age means age of infection , i.e., time elapsed since infection. For the sake of simplicity,we will assume that S is a discrete space. • A random point measure P describes the times of secondary infections. Due to theprevious assumptions, atoms of P all have mass 1 and only charge (0 , ∞ ) . • We can (and will) superimpose time heterogeneity to this process by means of a suppression function ( c ( t ); t ≥ valued in [0 , thinning the infection process.More precisely, if t is a potential time of infection for individual x (i.e., t is anatom of its infection point process P x ), we ignore the event with probability − c ( t ) . This suppression function can model the effect of vaccination, of density-dependence (i.e., relaxing the branching assumption due to an excess of removedor of deceased individuals), or of governmental mitigation measures (i.e., socialdistancing, lockdown).The population is thus described by a Crump-Mode-Jagers branching process , where allindividuals x are characterized by independent copies ( P x , X x ) of the pair ( P , X ) describ-ing, respectively, the process of infection and the life cycle. In the branching processliterature, X is often referred to as a random characteristic of individual x [4, 5, 6].4 emark 1. A typical infection measure is a Poisson Point Process with intensity λ killedat an independent exponential random variable with parameter γ . (By killing we meanthat we erase all atoms of the PPP after the killing time.) This corresponds to the classicalSIR process with rate of infection λ and recovery rate γ . More generally, we can constructa SEIR process (E for exposed) as follows. Let L be a random variable and consider ξ be a Poisson Point Process with an intensity measure λ ( x )d x . Then define P ([ a, b ]) = ξ ([( a − L ) + , ( b − L ) + ]) so that no infection occurs during the incubation period [0 , L ] . Remark 2. P and X are generally not independent . As a simple example, since (most)diseases cannot be transmitted by deceased individuals, no atom of P is allowed once X hasreached the end-state D . In the same spirit, one could assume that the infection potentialof a given individual is reduced once in the hospital and that individuals with many atomsin their infection process P (high infectiosity) are identified and isolated. Remark 3.
Classes such as P , I or H must sometimes be refined to account for ad-ditional structuring variables like general health condition, (real) age, spatial position orprevious exposure to similar pathogen. Knowledge of such variables can help predict moreaccurately the outcome of the infection and parametrize more precisely the infection pro-cess. Regarding this last point, note that the assumptions in force here allow for anyimplicit or explicit structure provided that transmission from an individual of type i toan individual of type j does not depend on j (but may depend on i , as we have seen).Relaxing this assumption would result in describing the large population limit by a mul-tidimensional PDE instead of a one-dimensional PDE (see Section 1.4). Also, note thatignoring structuring by a hidden variable such as spatial position or health condition canlead to difficulties in estimating sojourn times in each compartment (such as P , I or H )from clinical data, due to over- or under-representation in this compartment of subsets ofindividuals carrying certain values of the hidden variable. The stochastic epidemic models we consider here are fairly general and can exhibit quitecomplex dependencies (i) between states and time, due to the lack of any Markov-typeassumption, (ii) between states, due to possibly hidden structuring variables impactingthe life cycle, (iii) between state and infection rate, and (iv) between past and futureinfection events. The main message of this note is that despite this apparent complexity,most of this complexity vanishes when the size of the population is large. More specifically,we show that in the limit of large populations (obtained by starting from a large initialpopulation or as a consequence of natural exponential growth), the population of infectedsstructured by age (of the infection) can be described by a one-dimensional PDE that onlydepends on(a) The average infection rate τ (d a ) := E ( P (d a )) , and(b) The 1d marginals of the life-cycle process p ( a, i ) := P ( X ( a ) = i ) . arge initial population. Let us start with N infected individuals and define theempirical measure µ Nt describing the ages and types of infected individuals present attime t µ Nt (d a × { i } ) := (cid:88) x σ x
0) = c ( t ) (cid:90) ∞ τ (d a ) n ( t, a ) , ∀ a ≥ , n (0 , a ) = g ( a ) . (2) After lockdown.
Our second result (Theorem 8) displays a similar, but more subtle,convergence in the case when the process is supercritical, where natural growth leadsby itself to large population sizes. Let Z ( t ) denote the total population size at time t and assume that Z (0) = 1 , i.e., we start from a single individual . By a slight abuse ofnotation, denote by µ t K t the empirical measure of ages and types as in (1), but under theassumption that the suppression function at time t is equal to c ( t − t K ) where c is equalto 1 for negative arguments, and t K is some large, random time. We are motivated bymodeling a situation where the infection is separated into two distinct phases:(Phase 1) We let the epidemic develop until a certain random time t K . For instance, t K couldbe the time at which the number of recorded deaths exceeds a large threshold K .We assume no suppression before t K ;(Phase 2) We let the suppression function vary after time t K , e.g., due to mitigation measuresand/or behavioral changes (i.e., lockdown phase).Conditional on non-extinction, letting ( t K ) be any sequence of stopping times such that t K → ∞ on the non-extinction event, we have the following convergence in probability lim K →∞ µ t K t K + t Z ( t K ) = n ( t, a ) p ( a, i )d a. Now n ( t, a ) solves the McKendrick-von Foerster PDE with the same boundary conditionas previously, but with initial condition n (0 , a ) = αe − αa , where α > is the exponentialgrowth rate of Z ( t ) for t ≤ t K , also called Malthusian parameter. This second result canbe seen as a refinement of limit theorems for exponentially growing populations countedwith random characteristics, where here the characteristic of a typical individual is thenumber of her descendants in class i of age a , born at least s time units after her birth(summed over s x = t K − σ x ). In particular, taking t = 0 in the statement yields the6onvergence to the exponential stable age distribution decorated by the 1d marginals p of X . The way we state the result nicely displays dependencies between characteristicscorresponding to different t ’s.To summarize: (1) the macroscopic infection process is characterized by the sole in-tensity measure τ and dictates an explicit age structure of the population, and (2) theclass structure is deduced by averaging the life-cycle process against the limiting age pro-file . In order to validate our approach, we use those deterministic approximations toinfer epidemiological parameters ( R before and during lockdown) from recent empiricalobservations, and show that our findings are in accordance with the current literature. Outline.
The paper is organized as follows. In Section 2, we study the McKendrick-von Foerster equation (2). After providing a precise description of the branching processthat we are considering in Section 2.1, we give the definition of a weak solution to (2) inSection 2.2. Then, we give two probabilistic representations of these weak solutions. Wefirst show in Section 2.3 that the weak solution to (2) corresponds to the first momentof the branching process that we are studying, when viewed as a random measure on theages of infection. Second, Section 2.4 provides a construction of the weak solution using adual genealogical process. The two laws of large numbers are proved in Section 3. Finallythe inference in carried out in Section 4. Section 4.2 and Section 4.3 describe our choiceof parametrization and the inference results are discussed in Section 4.4.
Some of the assumptions underlying the previous models can be relaxed and our generalframework can be adapted to more complex and realistic populations.
Contact matrix.
So far, infectious individuals infect new individuals uniformly at ran-dom. In general, a contact matrix specifies the contact rate, depending on contact location(household, school, work...) or individual types of source and target (real age, suscepti-bility...) [7, 8]. More precisely, each individual now belongs to one class, and we denoteby C the (finite) set of all classes. An individual in class j ∈ C is characterized by amulti-dimensional process ( X j , P j, , . . . , P j,C ) where the atoms of P j,k provide the ageat which this individual infects a new individual of type k ∈ C , and ( X j ( a ); a ≥ is a S -valued process whose distribution can depend on j .We define the mean contact matrix as ∀ j, k ∈ C , τ jk (d u ) := E ( P j,k (d u )) . The population is again described by a multi-type Crump-Mode-Jagers branching process.Analogously to (1), µ N ,jt denotes the empirical measure of j individuals starting with N = ( N j ) j ∈C infected individuals at time t = 0 . Assume that there exists a constant N such that for all j ∈ C , N j /N → y j as N → ∞ . Under the usual technical conditions(the matrix τ is irreducible and Malthusian, x log( x ) type condition, see [9]), ∀ j ∈ C , ∀ i ∈ S , N j µ N ,jt (d a × { i } ) = ⇒ n j ( t, a ) P ( X j ( a ) = i )d a ( n j ) j ∈C satisfies the multidimensional McKendrick-von Foerster equation ∂ t n j ( t, a ) + ∂ a n j ( t, a ) = 0 ∀ t > , n j ( t,
0) = (cid:88) k ∈C c kj ( t ) (cid:90) ∞ τ kj (d a ) n k ( t, a ) , ∀ a ≥ , n j (0 , a ) = y j g j ( a ) . (3)where g j describes the initial age profile of class j and c jk ( t ) is a matrix at time t general-izing the suppression function of the previous section. Following Theorem 8, it is naturalto consider the initial condition g j ( u ) = α exp( − αu ) φ j where α is the Malthusian parameter, i.e. the unique α such that the largest eigenvalueof the matrix (cid:90) exp( − αu ) τ (d u ) is equal to and ( φ j ) j ∈C is its Perron-Frobenius left eigenvector with (positive) entriessumming up to . As in Theorem 8, such initial condition can be justified by starting witha single individual and let the population grow up to a large random time t K conditionalon non-extinction [9]. Shortage of susceptibles.
Another natural extension consists in taking into accountsaturation of infected in the population. Start with a finite, but large population of size N with a fraction of infected individuals x ∈ (0 , with an age profile g at time t = 0 . Hereinfection is only effective if the target individual is susceptible, which thins infection rateby the fraction of susceptibles in the population. At the limit, this saturation translatesinto a non-linear McKendrick-von Foerster equation ∀ t, a > , ∂ t n ( t, a ) + ∂ a n ( t, a ) = 0 ∀ t > , n ( t,
0) = S ( t ) c ( t ) (cid:90) ∞ τ (d a ) n ( t, a ) , ∀ a ≥ , n (0 , a ) = xg ( a ) , (4)where we have defined S ( t ) := 1 − (cid:90) ∞ n ( t, a )d a, and the limiting empirical measure is given by n ( t, a ) p ( a, j )d a . Convergence results tothis limiting PDE are addressed by some of the present authors in [10]. An important special case of our model is when the process ( X ( a ); a ≥ is a Markovprocess and infections from individual x occur at a rate that only depends on the currentstate of x .Under these assumptions, the McKendrick-von Foerster PDE reduces to a finite setof ODEs. Similar sets of ODEs have been widely used to model the SARS-CoV-2 epi-demic [11, 12, 13, 14], and in that sense, taking into account explicitly the (infection) age tructure of the population allows us to incorporate all these models into the same generalframework. More precisely, for i ∈ S define ∀ t ≥ , N i ( t ) = (cid:90) ∞ n ( t, a ) p ( a, i )d a to be the number of individuals in state i . We will assume that ( X ( a ); a ≥ is aMarkov process with transitions ( λ ij ) i,j ∈S . Moreover, we suppose that conditional on ( X ( a ); a ≥ , the infection process P x from individual x is a Poisson point process withintensity rate τ i when x is in state i . Then a direct computation shows the followingresult. Proposition 1.
Suppose that ( X ( a ); a ≥ is a Markov process with transitions ( λ ij ) i,j ∈S ,and that conditional on ( X ( a ); a ≥ , P is a Poisson point process with intensity ( τ X ( a ) ; a ≥ . Then, if ( n ( t, a )) denotes the solution to (4) , ( N i ( t ); t ≥ i ∈S solvesthe following set of ODE: ∀ i ∈ S , ˙ N i = (cid:88) j ∈S λ ji N j − N i (cid:88) j ∈S λ ij + (cid:88) j ∈S a ji N j S, (5) where a ji := τ j p (0 , i ) and S := 1 − (cid:80) i ∈S N i .Proof. Recall that N i ( t ) = (cid:90) ∞ n ( t, a ) p ( a, i )d a. By differentiating both sides with respect to time we get ˙ N i ( t ) = (cid:90) ∞ ∂ t n ( t, a ) p ( a, i )d a = − (cid:90) ∞ ∂ a n ( t, a ) p ( a, i )d a = (cid:90) ∞ n ( t, a ) ∂ a p ( a, i )d a + n ( t, p (0 , i ) . By using the boundary condition and the fact that ( X ( a )) a ≥ is a Markov process, weobtain that ˙ N i ( t ) = (cid:90) ∞ n ( t, a ) (cid:16) (cid:88) j ∈S λ ji p ( a, j ) − λ ij p ( a, i ) (cid:17) d a + S ( t ) p (0 , i ) (cid:90) ∞ n ( t, a ) (cid:88) j ∈S τ j p ( a, j )d a. = (cid:88) j ∈S λ ji N j ( t ) − N j ( t ) (cid:88) j ∈S λ ij + (cid:88) j ∈S a ji S ( t ) N j ( t ) , which ends the proof.Conversely, let us consider from the start the system of ODEs (5). Here, λ ij is inter-preted as the rate at which a type i individual turns into type j and a ij as the rate atwhich a type i individual gives birth to a type j individual. If the contact matrix withgeneric entries a ij has rank 1 and non-negative entries, it can always be decomposed as a ij = τ i p (0 , j ) where τ i ≥ and ( p (0 , j )) is a probability vector. The vectors ( τ i ) and ( p (0 , j )) can be recovered by τ i = (cid:88) j ∈S a ij and p (0 , j ) = a ij τ i ( ∀ i, j ) . a ij = λ ( i ) p (0 , j ) , so that actually a type i individual gives birth at rate λ ( i ) and her offspring has type independently distributed according to p (0 , · ) . As a result,the contact matrix with generic entries a ij has rank 1, and λ and p (0 , · ) can be recoveredby λ ( i ) = (cid:88) j ∈S a ij and p (0 , j ) = a ij λ ( i ) ( ∀ i, j ) . Then one can define ( X ( a ); a ≥ as the S -valued process with rates given by the matrix Λ (with diagonal entries λ ii = − (cid:80) j (cid:54)∈S λ ij ) and initial distribution ( p (0 , i )) . Denote as inthe rest of the text p ( a, i ) = P ( X ( a ) = i ) , so that the row matrix p ( a, · ) can be computed as the product p (0 , · ) exp( a Λ) .Now let us consider the solution ( N i ( t ); t ≥ to (5) and assume there is some agedistribution g (integrable but possibly not summing to 1) such that N i (0) = (cid:90) ∞ g ( a ) p ( a, i )d a. (6)Then by uniqueness of ( N i ) and thanks to Proposition 1, for all i ∈ S , N i ( t ) = (cid:90) ∞ n ( t, a ) p ( a, i )d a, where n is the solution to the McKendrick-von Foerster PDE with initial condition g andboundary condition n ( t,
0) = (cid:90) ∞ τ ( a ) n ( t, a )d a, where τ ( a ) := (cid:80) j ∈S τ j p ( a, j ) . This shows that the solution to any linear system of ODEsof the form (5) has a simple representation in terms of the solution to the McKendrick-von Foerster PDE decorated with types, provided there is a representation of the ini-tial condition in the form (6). Note that this last property is not necessarily fulfilled.For example, if X is ergodic and started in its stationary probability distribution, then p ( a, i ) = p (0 , i ) and (6) would only hold if ( n i (0)) i ∈S were proportional to ( p (0 , i )) i ∈S .If the matrix with generic entries a ij does not have rank 1, one could derive a similarrepresentation of the solutions to (5), but using the multi-dimensional version of theMcKendrick-von Foerster equation (3). Non-Markov epidemic models have already been investigated, see e.g. [15, 16, 17, 18].The idea of representing a general branching population by its age structure has a richhistory in probability theory [19, 20, 21, 22, 23, 24, 25, 26] and the connection withthe McKendrick-von Foerster PDE has been acknowledged several times [20, 24]. In thelatter two works, the authors allow for birth and death rates that may depend not onlyon abundances of each type, but also on the whole age structure of the population. Thisimpressive level of generalization comes at the cost of assuming that the process describingthe evolution of the empirical measure on ages and types is Markovian. In particular, birthand death rates are not allowed to depend on past individual birth events. The Markovproperty then allows the use of a generator for the empirical measure and with some10xtra finite second moment assumptions on the intensity measure, this approach allowsthe authors to obtain a law of large numbers and a central limit theorem.We acknowledge that the current work is certainly not as mathematically challengingas the works alluded to above, and that some of our results are almost implicit in someof the previous works. However, we believe that our point of view (1d PDE decoratedwith types) does deserve to be highlighted in the current sanitary crisis since it providesa natural and efficient inference methodology. More than 70,000 publications related tothe COVID-19 crisis have appeared since the onset of the pandemic, with many differentmodeling approaches. One of the modest aims of the present note is to convey the ideathat individual-based stochastic models suggest a simple and tractable framework fortackling some of the complex features of the disease. Furthermore, since we ignore finitepopulation effects, our proofs are quite elementary compared to [20, 24] and should beaccessible to a much wider audience interested in such a modeling approach.Finally, we already pointed out that the connection between branching processes andMcKendrick-von Foerster PDE has been discussed in previous works. However, as faras we can tell, the duality result exposed in Section 2.4 is new and can presumably beextended to more general branching processes where birth and death rates are allowed tobe frequency-dependent. In [10], some of the authors of the present work show that thisduality result has a natural counterpart in a model with a finite but large population.
CMJ branching process with suppression.
Recall that the infection process is mod-eled by a Crump-Mode-Jagers (CMJ) branching process [19, 4] with no death, startingfrom one individual called the progenitor (or root of the tree). Each individual x ischaracterized by an independent pair ( P x , X x ) embodying respectively the processes ofsecondary infection events from x and of types carried by x . Each pair ( P x , X x ) is a copyof the pair ( P , X ) with law L , except when x is the root, where it is distributed as ( ˜ P , ˜ X ) with law ˜ L (more on that below).Also recall some infection events can be suppressed using a suppression function ( c ( t ); t ≥ . Given a realization of the CMJ tree, for each branching point occurringat time t , we trim the tree by independently pruning the subtree stemming from it withprobability c ( t ) . See Figure 1.For simplicity, we will assume that the suppression function c is a piecewise right-continuous function, and that for any t ≥ , the process ( X ( a ); a ≥ is a.s. continuousat t . Define the average infection measure (that is, the intensity measure of the pointprocess P ) as τ (d u ) := E ( P )(d u ) . We assume that τ is absolutely continuous w.r.t. the Lebesgue measure in such a waythat there exists a measurable non-negative function β such that τ (d u ) = β ( u )d u and R := (cid:90) ∞ β ( u )d u < ∞ . (7)11 =0 Figure 1: The initial individual ( ˜ P , ˜ X ) is represented by a black segment. In Section 2.1,we assume that at time t = 0 , the age of the initial individual (length of the grey segment)is distributed according to a probability density g . If a branching event is observed attime t (see e.g., black dots), the infection occurs with probability c ( t ) . In the CMJ, thisamounts to prune the corresponding subtree with probability c ( t ) (dotted red tree).We also assume that there exists a unique parameter α ∈ R (the so-called Malthusianparameter of the CMJ process) such that (cid:90) exp( − αu ) τ (d u ) = 1 . (8)The parameter α can be either positive (supercritical) or negative (subcritical). Finally,we will also enforce the Kesten and Stigum criterium [27] E (cid:0) R α log + ( R α ) (cid:1) < ∞ , (9)where R α := (cid:90) ∞ exp( − αt ) P (d t ) . Initial shifted law.
For any finite measure m on R + , we define θ t ◦ m as the measureshifted by t , i.e., (cid:90) ∞ f ( u ) θ t ◦ m (d u ) = (cid:90) ∞ t f ( u − t ) m (d u ) , where f is any measurable, bounded function f on R + . For any measurable function F : R + → S , we similarly define θ t ◦ F by θ t ◦ F ( u ) = F ( u + t ) . (We make a slight abuseof notation by using the same symbol for the shift operator on measures and functions.)Let g be a probability density on R + . We now specify the law ˜ L g of the pair ( ˜ P , ˜ X ) characterizing the root. In order to connect the CMJ process with the McKendrick-von Foerster equation, we will focus on the case where ( ˜ P , ˜ X ) is identical in law to ( θ A ◦ P , θ A ◦ X ) , where A is a r.v. independent of ( P , X ) and distributed according to g .The distribution ˜ L g will be referred to as the shifted law . In particular, we have ˜ τ (d a ) := E ( ˜ P (d a )) = (cid:18)(cid:90) ∞ β ( x + a ) g ( x )d x (cid:19) d a. otation. We assume that individuals are indexed by the standard Ulam-Harris label-ing. Namely, individuals are indexed in I = ∪ n ( N ∗ ) n . If x ∈ I , then xi (the concatenationof x and i ) is interpreted as the i -th child of x . Children are ranked according to theirbirth time: ( x, is the oldest child of x , ( x, the second oldest, etc . (Since we assumedthat τ has a density, there is no simultaneous births and the atoms of P are distinct a.s.)We denote by σ x the date of birth of x with the convention that σ x = ∞ if the individualis never born. For instance, if σ x < ∞ and x has k children, then σ xj = ∞ for j > k .Finally, ∅ will denote the root of the tree. In this section, we consider the time-inhomogeneous, linear McKendrick-von FoersterPDE (2). The first line in (2) is the transport equation with unit velocity, i.e., agesof individuals increase at rate . The second line gives the number of newly infectedindividual (age ) at time t .In order to motivate our definition of weak solutions, we start by giving a formalresolution of the PDE using the method of characteristics. Fix a > . Let A ( t ) = a − t Then dd s n ( t − s, A ( s )) = − ∂ t n ( t − s, A ( s )) − ∂ a n ( t − s, A ( s )) = 0 , so that s (cid:55)→ n ( t − s, a − s ) is conserved along the characteristics, i.e., ∀ s < a, n ( t, a ) = n ( t − s, a − s ) . It follows that n ( t, a ) = (cid:40) g ( a − t ) when a > tb ( t − a ) when a ≤ t (10)where b ( t ) = c ( t ) (cid:82) ∞ τ (d a ) n ( t, a ) . We now determine the function b . Injecting the pre-vious expression into the ‘spatial’ boundary condition of the PDE, we obtain a delayedequation for b : for every t > b ( t ) = c ( t ) (cid:90) t τ (d a ) b ( t − a ) + c ( t ) (cid:90) ∞ t τ (d a ) g ( a − t ) . (11) Lemma 2.
There exists a unique solution b to (11) which is locally integrable. Moreover,for any δ ≥ such that δ > α we have b ∈ L ,δ , where L ,δ denotes the set of all functions f : R + → R such that (cid:107) f (cid:107) L ,δ := (cid:82) ∞ e − δt | f ( t ) | d t < ∞ .Proof. Fix δ > α and denote by L ,δ the space L ,δ quotiented by the relation ∼ δ , where f ∼ δ g if (cid:107) f − g (cid:107) L ,δ = 0 . Then define the linear operator Φ : L ,δ → L ,δ by Φ f : t (cid:55)→ c ( t ) (cid:90) t f ( t − u ) β ( u )d u. (cid:107) Φ f (cid:107) L ,δ = (cid:90) ∞ e − δt Φ f ( t )d t = (cid:90) ∞ e − δt c ( t ) (cid:90) t f ( t − u ) β ( u )d u d t = (cid:90) ∞ e − δu f ( u ) (cid:90) ∞ u β ( t − u ) e − δ ( t − u ) c ( t )d t d u. Now using that (cid:90) ∞ u β ( t − u ) e − δ ( t − u ) c ( t )d t ≤ (cid:90) ∞ β ( t ) e − δt d t < we obtain that (cid:107) Φ (cid:107) < . Define Ψ := Id − Φ . Then Ψ is invertible with inverse (cid:80) k ≥ Φ k . Note that equation (11) can be written as Ψ( b ) = F, where F : t (cid:55)→ c ( t ) (cid:90) ∞ t β ( a ) g ( t − a )d a. The proof ends noting that F ∈ L ,δ as (cid:90) ∞ e − δt F ( t )d t ≤ (cid:90) ∞ (cid:90) ∞ t β ( a ) g ( t − a )d a d t < ∞ . We have thus proved existence and uniqueness of the solution b to (11) in L ,δ . Now forany two functions b and b such that b ∼ δ b and b and b both satisfy (11), we have b = b (i.e., there is a single element in the equivalence class of b ). This shows uniquenessof the solution b to (11) in L ,δ .Since all elements of L ,δ are locally integrable, this also shows the existence of alocally integrable solution to (11). Its uniqueness can be proved following the exactsame reasoning as previously, replacing integrations on [0 , ∞ ) by integration on compactintervals. Definition 3.
We say that ( n ( t, a ); a, t ≥ is the L , loc weak solution to the McKendrick-von Foerster PDE with initial condition g if it satisfies the relation (10) where ( b ( t ); t ≥ is the unique locally integrable solution to (11) displayed in the previous lemma. Define Z ( t ) := (cid:88) x ( σ x ∈ (0 , t ]) , B ( t ) := E (cid:0) Z ( t ) (cid:1) where Z ( t ) is interpreted as the number of infections between and t . Recall that R := (cid:82) ∞ β ( u )d u < ∞ guarantees that B ( t ) < ∞ for all t ≥ . Finally, B is non-decreasingand we denote by d B the Stieljes measure associated to B . Lemma 4.
There exists a locally integrable function ( b ( t ); t ≥ such that d B ( t ) = b ( t )d t. Further, b coincides with the unique locally integrable solution of the delayed equation (11) . roof. The fact that d B has a density easily follows from the fact that τ has a density.The fact that B ( t ) < ∞ ensures that b is locally integrable.Define ¯ P x the infection measure obtained from P x after random thinning by the func-tion ( c ( t ); t ≥ . Namely, conditional on σ x and the atoms a < a < · · · of P x , weremove independently each of the atoms with respective probabilities − c ( σ x + a ) , − c ( σ x + a ) . . . , whereas the other atoms remain unchanged.Fix t > . Let k ≤ n ∈ N . Define T k,n ( P x ) as the measure obtained from P x asfollows. Conditional on the atoms a < a < · · · of P x , we remove independently each ofthe atoms with respective probabilities − max z ∈ ( t kn ,t k +1 n ] c ( z + a ) , − max z ∈ ( t kn ,t k +1 n ] c ( z + a ) · · · and leave other atoms unchanged. Note that the thinning procedure is now independentof the starting time σ x . Further, if σ x ∈ ( t kn , t k +1 n ] , the point measure T k,n ( P x ) dominates ¯ P x .We decompose the births on (0 , t ] into two parts: individuals stemming from the root ∅ and a second part from subsequent births. Using the fact that for every individual x ,the (un-suppressed) random measure P x is independent of its birth time σ x (see secondequality below), and setting M ( t ) := (cid:82) t (cid:82) ∞ g ( a ) β ( a + u ) c ( u )d a d u , we get B ( t ) = n − (cid:88) k =0 (cid:88) x E (cid:18) (cid:18) σ x ∈ ( t kn , t k + 1 n ] (cid:19) (cid:90) [0 ,t − σ x ] ¯ P x (d a ) (cid:19) + M ( t ) ≤ n − (cid:88) k =0 (cid:88) x E (cid:32) (cid:18) σ x ∈ ( t kn , t k + 1 n ] (cid:19) (cid:90) [0 ,t − t kn ] T k,n ( P x )(d a ) (cid:33) + M ( t )= n − (cid:88) k =0 (cid:88) x E (cid:18) (cid:18) σ x ∈ ( t kn , t k + 1 n ] (cid:19)(cid:19) E (cid:18) (cid:90) [0 ,t − t kn ] T k,n ( P )(d a ) (cid:19) + M ( t )= n − (cid:88) k =0 E (cid:32)(cid:88) x (cid:18) σ x ∈ ( t kn , t k + 1 n ] (cid:19)(cid:33) E (cid:18) (cid:90) [0 ,t − t kn ] T k,n ( P )(d a ) (cid:19) + M ( t )= n − (cid:88) k =0 (cid:18) B ( t k + 1 n ) − B ( t kn ) (cid:19) E (cid:18) (cid:90) [0 ,t − t kn ] T k,n ( P )(d a ) (cid:19) + M ( t )= n − (cid:88) k =0 (cid:18) B ( t k + 1 n ) − B ( t kn ) (cid:19) (cid:90) [0 ,t − t kn ] c k,n ( u ) β ( u )d u + M ( t ) . with c k,n ( y ) = max v ∈ ( t kn ,t k +1 n ] c ( y + v ) . In particular, if tk/n → x , and x + y is a continuitypoint of c , we have c k,n ( y ) → c ( x + y ) . We will pass to the limit n → ∞ in the latterinequality. Recall that c is bounded (and valued in [0 , ) and that c is right-continuous.The first term on the RHS can be written under the form n − (cid:88) k =0 (cid:18) B ( t k + 1 n ) − B ( t kn ) (cid:19) (cid:90) t − [ x ] n c k,n ( u ) β ( u )d u = (cid:90) t f n ( x )d B ( x ) , where f ( n ) ( x ) = (cid:90) t − [ x ] n sup v ∈ ([ x ] n , [ x n ]+ tn ] c ( v + u ) β ( u )d u and [ x ] n = tn (cid:98) nx/t (cid:99) .
15e will now apply twice the Bounded Convergence Theorem. On the one hand, for afixed value of x , as n → ∞ [0 ,t − [ x ] n ] ( u ) sup v ∈ ([ x ] n , [ x ] n + tn ] c ( v + u ) β ( u ) → [0 ,t − x ] ( u ) c ( x + u ) β ( u ) Lebesgue a.e.Further, the latter term (i.e., the integrand in the integral defining f n ) is uniformlybounded by β and (cid:82) ∞ β ( u ) du < ∞ . A first application of the Bounded ConvergenceTheorem implies that for every x , as n → ∞ f n ( x ) → (cid:90) t − x c ( x + u ) β ( u )d u. On the other hand, the uniform bound, f n ( x ) ≤ R = (cid:82) ∞ β ( u )d u for all x, n , allows usto again apply the Bounded Convergence Theorem, so we get B ( t ) ≤ (cid:90) t b ( x )d x (cid:90) [0 ,t − x ] c ( x + u ) β ( u )d u + (cid:90) t (cid:90) ∞ g ( a ) β ( a + u ) c ( u )d a d u. By an analog argument, one can establish the same lower bound and strengthen the latterinequality into an equality. After a simple change of variable and interchanging the orderof integration, this yields B ( t ) = (cid:90) t c ( v ) (cid:90) v β ( v − x ) b ( x )d x d v + (cid:90) t (cid:90) ∞ g ( a ) β ( a + u ) c ( u )d a d u. Finally, differentiating with respect to t yields the desired result. Corollary 5 (Forward Feynman-Kac formula) . For every t ≥ , define ¯ µ t (d a × { i } ) := n ( t, a ) × P ( X ( a ) = i )d a, where n is the unique L ,loc weak solution to the McKendrick-von Foerster PDE with initialcondition g . Then ¯ µ t (d a × { i } ) = E (cid:18) (cid:88) x σ x
We note that there are as many dual processes as there are point processeswith intensity measure τ . Here are a few natural choices:1. M = P .2. Let M be a Poisson Point Process with intensity measure τ (d u ) . In this particularcase, the dual process is a Bellman-Harris branching process (i.e., the offspringlifetime durations are independent conditional on offspring number). We note that τ (d u ) appears naturally when considering the ancestral spine of a critical CMJ, seee.g. [28]. The measure τ can be obtained by size-biasing P (i.e. biasing by the totalmass of P ) and then recording the age of the individual at a uniformly chosen birthevent. Let ( Y t ; t ≥ be the stochastic process valued in ∪ n ∈ N R n + recording the residual life-times at time t listed in increasing order, i.e. if Y t = ( Y (1) t , · · · , Y ( n ) t ) there are n particlesalive at time t and Y ( k ) t is the residual life-time of the k th -individual with Y (1) t < · · · < Y ( n ) t .(We assumed that τ has a density so that the residual lifetimes are distinct a.s.). Inparticular, the particle labelled at any given time t will be the first to expire, and atdeath time t + Y (1) t a random number of children is produced. We let dim( Y T ) denote thenumber of particules alive at time t , i.e., the dimension of the vector Y t .Finally, we will say that n is a right-continuous version to the McKendrick-von Foersterequation, if n is a L ,loc weak solution and for every T, x ≥ the function s → n ( T − s, x − s ) is right-continuous on [0 , x ] . 18 roposition 6 (Backward Feynman-Kac formula) . Assume that the suppression function ( c ( t ); t ≥ is right-continuous. Then there is a unique right-continuous solution to theMcKendrick-von Foerster equation, and for every a, T ≥ n ( T, a ) = (cid:98) E a (cid:32) (cid:88) i ≤ dim( Y T ) g ( Y ( i ) T ) (cid:33) (13) where (cid:98) E a is the distribution of the process ( Y t ; t ≥ starting with an individual withresidual lifetime a .Proof. We first assume that there exists a right-continuous solution to the PDE. Let t < · · · < t k < · · · be the successive branching times of the dual branching process.Since τ has a density, there is a single branching particle at the successive branchingtimes t , · · · . Define the càdlàg process Z s := (cid:88) i ≤ dim( Y s ) n ( T − s, Y ( i ) s ) See also Figure 2.4 for a pictorial representation of the process. It is plain from thedefinition that n is preserved along the characteristics of the PDE, i.e., that for every x the function s → n ( T − s, x − s ) remains constant on [0 , x ) . As a consequence, ( Z s , s ≥ remains constant on every interval [ t n , t n +1 ) , with the convention t = 0 . Define z n := Z t n the value of the process ( Z t ; t ≥ at the n -th branching time. Let ( F n , n ∈ N ) be thefiltration induced by the process ( z n ; n ∈ N ) . For every n > E a ( z n | F n − ) = (cid:88) ≤ i ≤ dim( z n ) n ( T − t n , Y ( i ) t n ) + c ( T − t n ) (cid:90) ∞ n ( T − t n , a ) τ (d a )= (cid:88) ≤ i ≤ dim( z n ) n ( T − t n , Y ( i ) t n ) + n ( T − t n , z n − , where the second equality follows from the spatial boundary of the McKendrick-von Fo-erster equation. As already mentioned, the process ( Z s ; s ≥ is constant between twobranching times. As a consequence, ( Z s ; s ≥ is a martingale (w.r.t. its natural filtra-tion) so for every s ≥ , n ( T, a ) = ˆ E a (cid:18) (cid:88) i ≤ dim( Y s ) n ( T − s, Y ( i ) s ) (cid:19) . Relation (13) follows by taking s = T in the latter expression.The proof ends by checking that when c is right-continuous, the RHS of (13) indeedis a right-continuous solution to the PDE. This elementary step is left to the interestedreader. 19 a(a,T)T- 𝙩 𝙩 𝙩 Figure 2: Graphical representation of the process ( Z s ; s ≥ . We start with a singleindividual with residual lifetime a . In this picture, time flows downwards for the branchingprocess. The residual lifetime of the initial individual decreases linearly until reaching (this corresponds to time T − t in our representation). At this time, the particle diesand produces 2 red particules. Residual lifetimes travel along the characteristics of theMcKendrick-von Foerster PDE until reaching the spatial boundary condition { a = 0 } where a new branching occurs. Theorem 7 ( N individuals) . Start with N individuals at time with independent infectionand life-processes distributed according to the initial shifted law ˜ L g . Define the empiricalrandom measure for ages and types at time tµ Nt (d a × { i } ) := (cid:88) x σ x
Take t K = inf { t > { x ∈ ∪ n N n : σ x < t, X x ( t − σ x ) = D } ≥ K } , i.e., t K is the first time that the accumulated number of deaths reaches level K .Further take C ( t ) = 1 if t ≤ and C ( t ) = r < if t > . This corresponds to alockdown strategy where transmission is reduced by a factor r upon reaching K deaths. Theorem 8 (One individual) . Conditional on non-extinction • There exists a r.v. W ∞ such that W ∞ > a.s. and (cid:88) x (0 < σ x < t K ) exp( − αt K ) −−−→ t K →∞ W ∞ a.s. and in L Fix t ≥ . We have exp( − αt K ) µ t K t K + t === ⇒ t K →∞ W ∞ ¯ µ t in probability , where the convergence is meant in the weak topology, and ¯ µ t (d a × { i } ) = n ( t, a ) × P ( X ( a ) = i )d a with n the unique L , loc weak solution to the McKendrick-von Foerster PDE withinitial condition g ( a ) = α exp( − αa ) and suppression function ( C ( t ); t ≥ .Proof. We recall some basic facts about the method of random characteristics. We con-sider a plain CMJ with no death (no suppression function, no types). Every individualis characterized by an independent pair of random variables ( P , χ ) . As before, P is theinfection measure recording the times of secondary infections. Now χ is a general stochas-tic process indexed by the age of the individual called a random characteristic with theconvention χ ( − a ) = 0 for a ≥ . In great generality, P and χ may exhibit a non-trivialcorrelation. Define Z χ ( T ) = (cid:88) x χ x ( T − σ x ) the branching process counted by the random characteristic χ at time T . We now recallone of the main results of Jagers and Nerman [5] (see also Theorem 5, Appendix A in[6]). On top of all the assumptions above, we make the two following extra assumptions(a) There exists α (cid:48) < α such that E (cid:18) sup T ≥ exp( − α (cid:48) T ) χ ( T ) (cid:19) < ∞ . (16)(b) T → E ( χ ( T )) is continuous a.e. with respect to the Lebesgue measure.Then there exists a positive r.v. W ∞ (independent of the choice of χ ) such that conditionalon non-extinction Z χ ( T ) exp( − αT ) → W ∞ (cid:90) ∞ α exp( − αt ) E ( χ ( t ))d t a.s. and in L (d x ) as T → ∞ .(Note that the L convergence holds thanks to the x log( x ) condition (9).)To illustrate the method, we recall that if we take χ ( T ) = R + ( T ) then Z χ ( T ) coin-cides with B ( T ) , the total number of births before time T . For this particular choice of(deterministic) characteristic, the two properties above are immediately satisfied (recallthat α > ), so that conditional on non-extinction (cid:88) x (0 < σ x < u ) exp( − αu ) → W ∞ a.s. and in L (d x ) as u → ∞ . (17)The a.s. convergence ensures that the first item of our theorem is satisfied.Next, the second part of the theorem requires a choice of characteristic, called ‘generalcharacteristic’, that depends on the descendance of each extant individual at time t K .Because we need to prove an a.s. convergence result, whereas limit theorems on branchingprocesses counted with general characteristics only hold in distribution, we have to design22 T+t T+s
Figure 3: The individual x with characteristic ( P x , X x ) under consideration is representedby a black segment. We graft independent ( P ∗ , X ∗ ) -CMJ processes with suppressionfunction t (cid:55)→ C ( t − T ) to the atoms of P x occurring at time T + s , independently withprobability C ( s ) if s ≥ and 0 if s < . We ignore all other atoms and their descendances(lower dotted red tree s < , upper dotted tree s > ).by hand an individual characteristic that has the same distribution as the requestedgeneral characteristic.In order to define our next random characteristics, we start with some definition. Let ( P ∗ , X ∗ ) be a pair of infection and life-cycle processes (that may or may not be identicalto ( P , X ) in distribution). One can construct a collection ( Z (∆) ; ∆ ≥ of ( P ∗ , X ∗ ) -CMJprocesses starting at respective times ∆ and with a suppression mechanism C , i.e., • Z (∆) is a ( P ∗ , X ∗ ) -CMJ process starting at time ∆ with one progenitor. • The suppression function applied to infection events is C . • The previous rule applies at time t = ∆ , that is, Z (∆) is identically empty withprobability − C (∆) .For any t ≥ ∆ , we let ν (∆) ( t ) denote the empirical measure for ages and types of Z (∆) ( t ) .Finally, we define { ( ν (∆) i ; ∆ ≥ } i ∈ N ∗ as the collection made of independent copies of thecollection ( ν (∆) ; ∆ ≥ .We are now ready to construct our random characteristics by enlarging the initialCMJ process in the following way. Fix t ∈ R + and f a bounded non-negative continuousfunction on R × S with compact support in R + × S . Consider a typical individual ∅ , withinfection and life processes ( P , X ) . Denote by ( r i ) the atoms of P listed in increasingorder. For any T ≥ , define the individual random characteristic χ ( t,f ) ( T ) , by χ ( t,f ) ( T ) := f ( T + t, X ( T + t )) + (cid:88) i : r i ∈P∩ [ T,T + t ] (cid:68) ν ( r i − T ) i ( t ) , f (cid:69) , See Figure 3 for a intuitive constructing of the random characteristics.From now on, we assume that ( P ∗ , X ∗ ) is identical in law to ( P , X ) , which impliesthe following two crucial facts. 23i) E ( (cid:82) ∞ χ ( t,f ) ( a ) g ( a ) da ) = (cid:104) ¯ µ t , f (cid:105) where ¯ µ t is defined as in Corollary 5 with initialcondition g and suppression function C .(ii) Let us now count our branching process by its random characteristic Z χ ( t.f ) ( T ) = (cid:88) x χ ( t,f ) x ( T − σ x ) . Since t K is a stopping time with respect to the filtration ( F t ; t ≥ , the branchingproperty implies that Z χ ( t,f ) ( t K ) is identical in distribution to (cid:10) µ t K t K + t , f (cid:11) where weremember that µ t K t is the empirical measure w.r.t. the CMJ process with suppressionfunction c K : t (cid:55)→ C ( t − t K ) .In order to apply the aforementioned result of Jagers and Nerman, we need to checkthat condition (a), (b) above are satisfied. Condition (b) easily follows from the factthat τ has a density. We now check the first condition. Consider a ( P , X ) -CMJ process,assuming that the initial individual is un-shifted. For every s , let v s (d a ×{ i } ) the empiricalmeasure of ages and types at time s . We assumed f to be non-negative and thus E (cid:18) sup s ∈ [0 ,t ] (cid:104) v s , f (cid:105) (cid:19) ≤ (cid:107) f (cid:107) ∞ B ( t ) . Let a < a < · · · be the atoms of P between T and T + t . The random characteristicunder consideration is obtained by attaching independent (suppressed) CMJ processesbetween T and T + t to the a i ’s and by summing up the respective empirical measures attime T + t − a i . By construction, T ≤ T + t − a i ≤ T + t so that E (cid:18) sup T ≥ exp( − α (cid:48) T ) χ ( T ) (cid:19) ≤ E (cid:32) sup s ∈ [0 ,t ] (cid:104) v s , f (cid:105) (cid:33) E (cid:18) sup T ≥ exp( − α (cid:48) T ) (cid:90) T + tT P ( du ) (cid:19) = (cid:107) f (cid:107) ∞ B ( t ) E (cid:18) sup T ≥ exp( − α (cid:48) T ) (cid:90) T + tT P ( du ) (cid:19) ≤ R B ( t ) (cid:107) f (cid:107) ∞ . This shows that property (a) is satisfied for any ≤ α (cid:48) < α . This shows that as v → ∞ ,conditional on non-extinction Z χ ( t,f ) ( v ) exp( − αv ) → W ∞ E (cid:18)(cid:90) ∞ α exp( − αu ) χ ( t,f ) ( u ) (cid:19) a.s.As already pointed out in (i), E (cid:18)(cid:90) ∞ α exp( − αu ) χ ( t,f ) ( u )d u (cid:19) = (cid:104) ¯ µ t , f (cid:105) where ¯ µ t is defined as in Corollary 5 with initial condition α exp( − αt ) and suppressionfunction ( C ( t ); t ≥ . Since the latter convergence is a.s. and Z χ ( t,f ) ( t K ) is identical inlaw with (cid:10) µ t K t K + t , f (cid:11) (see (ii) above), the result follows.24 Inference
In this section, we illustrate how to use our framework to make inferences from macro-scopic observables of the epidemic, e.g., incidence of positively tested patients, hospital orICU (intensive care unit) admissions, deaths, etc . We show how to use those observablesto extract the underlying age structure of the population, estimate model parameters andforecast the future of the epidemic.We focused on a longitudinal case study in France. From March 18 2020, the Frenchgovernment has provided daily reports of the numbers of ICU and hospital admissions, ofdeaths, of discharged patients, and of occupied ICU and hospital beds. Moreover, severaltheoretical studies have already been conducted on the same dataset. This allowed usto fix the values of some crucial biological parameters that had already been estimatedand to carry out a comparison with our method. We want to emphasize that the aimof this section is to provide a mathematical framework in which convergence results canbe rigorously proved while remaining flexible enough for other applications. Our goal isnot to provide new estimates of epidemiological parameters for France, as many robustestimates are already available. For instance we do not provide confidence intervals forour estimates, and neither do we conduct a sensibility analysis.The remainder of the section is laid out as follows. In Section 4.2 we identify themathematical quantities that impact the dynamics of the epidemic for large populationsizes. Section 4.3 then presents the choice of distribution we made for these quantitiesand the parameters that need to be estimated. Finally, estimation of these parametersfrom the French incidence data is performed in Section 4.4. We start by fitting a simplemodel and then show how this model can be made more complex to account for morecomplex incidence data.
As mentioned previously, the age structure of the population cannot be directly accessed.What is observed is a subset of the population with some characteristic of interest, forinstance individuals that have been tested positively, deceased individuals or dischargedpatients. Recall that under the assumptions stated in Section 3, the number of individualsthat are in a given state i at time t converges to (cid:90) ∞ n ( t, a ) P ( X ( a ) = i )d a, where ( n ( t, a )) is the solution to the McKendrick-von Foerster equation. Note that theassumptions of Theorem 8 are in essence that the epidemic has been ongoing for a longenough time at the lockdown onset for the infected population to be large, which weassume to hold true for France as the number of infected individuals on March 16 2020was on the order of tens of thousands of individuals.The McKendrick-von Foerster equation is determined by two quantities: (i) the averageinfection measure τ defined in 2.1 and (ii) an initial condition, which is of the form ∀ a ≥ , n (0 , a ) = W αe − αa W and a parameter α which correspondsto the exponential growth rate of the epidemic before the lockdown onset.Therefore, using Theorem 8 to obtain a theoretical prediction for some observablesunder consideration requires the knowledge of:1. The intensity measure τ of secondary infections;2. The initial number of infected W and the parameter α ;3. For each state i of interest the probability P ( X ( a ) = i ) .The next section exposes how we have parametrized these quantities. Average infection measure.
Recall the definition of τ from Section 2.1. Let us furtherdefine R = τ ([0 , ∞ )) , ˆ τ (d a ) = τ (d a ) τ ([0 , ∞ )) , so that τ can be expressed as τ (d a ) = R ˆ τ (d a ) . The total mass of τ , R , is the mean number of secondary infections induced by a singleinfected individual. Thus R corresponds to the basic reproduction number of the epi-demic, and we leave it as a parameter to infer. The epidemiological interpretation of theprobability measure ˆ τ is the following. Consider a large population of infected individuals.Then, as the size of that population goes to infinity, the distribution of the time betweenthe infection of a uniformly sampled individual and the infection of its “parent” convergesto ˆ τ . Therefore, ˆ τ is the distribution of the so-called generation time of the epidemic,which has already been estimated in several previous studies. We used the estimationof [29], and assumed that ˆ τ is a Weibull distribution, that is ˆ τ (d a ) = kλ (cid:16) aλ (cid:17) k − e − ( a/λ ) k d a, (18)where the values of the shape parameter k and scale parameter λ are recalled in Table 1. Initial condition.
The growth rate α is defined implicitly through equation (8). More-over, we know that α corresponds to the exponential growth rate of any of the observablesof the epidemic before the lockdown onset. We chose to estimate α from the cumulativenumber of deaths, which appeared to be more reliable than the number of positive testsas the number of tests that have been conducted in the early phase of the epidemic inFrance varied greatly. We simply estimated α using a linear regression on the logarithmof the number of deaths from March 7 to March 20, 2020. The estimated α as well as thecorresponding R pre-lockdown are shown in Table 1. The number of infected individualsat the time of lockdown was left as a parameter to infer.26 ife-cycle. The last quantities that need to be defined are the one-dimensional marginalsof the life-cycle process ( X ( a ); a ≥ . From now on, we assume that the sequence of statesvisited by ( X ( a ); a ≥ is a Markov chain and that the sojourn times in each state isGamma distributed.More specifically, we suppose that we are given a global dispersion parameter γ > and that for each state i ∈ S , the time D i spent in state i follows a Gamma distributionwith expectation t i > and variance t i /γ , • E ( D i ) = t i • Var( D i ) = t i /γ This assumption has the following consequence. Let i , . . . , i k be a possible sequence ofstates of ( X ( a ); a ≥ , and denote by T i , . . . , T i k the respective entrance times in thesestates. Then conditional on ( X ( a ); a ≥ successively visiting the states i , . . . , i k in thisorder, we have ( T i , . . . , T i k ) ∼ (0 , Y t , . . . , Y t + ··· + t k − ) where ( Y t ) t ≥ is a Gamma process such that Y t ∼ γ γt Γ( γt ) a γt − e − γa d a. Another advantage of this parametrization is that the 1d marginals of ( X ( a ); a ≥ canbe computed efficiently, while only requiring one parameter for each state of interest, anda global dispersion parameter. In this section, we fit six time series of the French epidemic: the number of ICU andhospital admissions, the number of deaths, the number of occupied ICU and hospitalbeds and the number of discharged patients. We provide two examples of Markov chainsthat we use to fit these data. The first Markov chain is only used to fit the first threecurves, that we will call incidence curves , i.e., the daily number of ICU admissions, ofhospital admissions and of deaths. With the parametrization of the previous section, thepredictions of our model for these time series only involves the delay between infectionand death, ICU, or hospital admissions. We do not need to estimate the time betweenhospitalization and discharge, which makes the model and the inference procedure easier.Then we show how this model can be made more complex to fit the last three curves,that we call prevalence curves . All parameter estimations were realized using the datafrom March 18 to May 11, 2020. This time frame corresponds to the lockdown period inFrance.
Fitting incidence data.
In order to fit the incidence curves, we considered the simplestmodel that can account for these three time series. The model is illustrated in Figure 4.Upon infection, with probability − p hosp , an individual develops a mild form ofCOVID-19 and is placed in state I , which encompasses all cases that do not require ahospitalization. With probability p hosp the individual has a severe infection and is placedin state C . Individuals in state C are eventually hospitalized and moved to state H .27 RR R RDUHIC i n f e c t i o n h o s p i t a l a d m i ss i o n I C U a d m i ss i o n d e a t h T H T U T D p hosp p hosp p ICU p ICU p death p death age ofinfection Figure 4: Illustration of the admission model.Then, with probability p ICU individuals in state H are admitted in ICU and move to state U . Otherwise they eventually recover and are discharged. Finally, individuals in state U die with probability p death , or recover with probability − p death . In this model, onlyindividuals in ICU may die.As we are fitting the number of individuals that enter a state, and not the number ofindividuals that are currently in that state, we only need to track the times T H , T U , and T D elapsed between infection and hospital admission, ICU admission and death, respectively.We will refer to this model as the admission model .Estimations for p hosp , p ICU and of death probability conditional on hospitalization(equal in our setting to p ICU × p death ) in France have already been conducted in [11].We used these estimates and considered the values of p hosp , p ICU and p death to be fixed.All other parameters were estimated using a maximum likelihood procedure which isdescribed in Appendix A. The parameter estimations are provided in Table 2, and thecorresponding predicted values for the time series under consideration are displayed inFigure 5. Overall, our simple model seems to match the observed data. Note howeverthat the model overestimates the number of ICU admissions in the second part of thelockdown. This is likely due to a temporal reduction in the ICU admission probabilitywhich has been reported in [11].Our estimation of the basic reproduction number during the lockdown period is R =0 . . This suggests that lockdown has reduced the basic reproduction number by a factor . compared to the beginning of the epidemic. Moreover, we estimated that . × infections have occurred in France before March 17th. Both these values are in line withprevious estimates for France [30, 11].We did not impose that T H < T U in the inference procedure. Interestingly we foundthat the data are best explained by assuming that the mean of T H is . days, whereasthe mean of T U is . days. This indicates that the delay between infection and hospital28igure 5: Best fit of the admission model. Solid line correspond to the number of hospitaladmissions, ICU admissions and deaths predicted by the admission model. The dots arethe corresponding observed values. Notation Description Value Source α Pre-lockdown exponential growth rate 0.315 E R pre Basic reproduction number before lock-down 3.25 E k Shape parameter of the generation time 2.83 [29] λ Scale parameter of the generation time 5.67 [29]Table 1: Parameter values common to both models. In the “Source” column, “E” indicatesthat the parameter has been estimated in the present work.admission is shorter for individuals that end up in ICU, compared to the average timebetween infection and hospitalization. Therefore it would be more appropriate to allowindividuals to have an admission to hospital delay that is different depending on whetherthey will end up in ICU or not, modeling the fact that they have a more severe form ofthe disease. We estimated the mean of T D , the time between infection and death, to be . days. This estimate is lower than but consistent with previous estimates based onthe study of individual-case data [3, 31, 32]. Fitting prevalence data.
A first attempt to fit the prevalence curves could be to keepthe admission model of Figure 4 and to estimate the time between hospital admissionand discharge using the observed number of occupied ICU, hospital beds, and dischargedpatients. However this only yields a poor fit of the data (see Appendix B). We identifiedtwo main reasons for this discrepancy. First, we assumed that all individuals are admittedto ICU prior to death. Using the probability estimated in [11] then yields that theprobability of dying conditional on being in ICU is . . This value is unrealisticallyhigh, and we need to assume that a fraction of hospital deaths occur without going throughthe ICU. Second, under the admission model, the delay between hospital admission anddischarge is almost unimodal. However, the observed number of occupied hospital bedsrises fast but falls slowly. Such a shape cannot be easily accounted for by a unimodaldistribution for the time spent in hospital.Taking into account the previous two points required us to make the model more com-29 otation Description Value Source R Basic reproduction number duringlockdown . E W Total number of infections beforeMarch 17 2020 . × E p hosp Probability of being hospitalized . [11] p ICU
Probability of entering ICU conditionalon being at the hospital . [11] p ICU · p death Death probability conditional on beinghospitalized . [11] T H Delay between infection and hospitaladmission . days E T U Delay between infection and ICU ad-mission . days E T D Delay between infection and death . days E γ Scale parameter common to all Gammadistributions 0.463 ETable 2: Inferred parameter set for the admission model. The values indicated for thedurations correspond to the means of the Gamma distributions. In the “Source” column,“E” indicates that the parameter has been estimated in the current work.plex. The resulting model, referred to as the occupancy model , is illustrated in Figure 6.We now consider that upon infection, individuals go to one of three states depending onthe severity of their infection: • The state C u which gathers critical infections that lead to death or ICU admission.The probability of having a critical infection is denoted by p crit . • The state C h which corresponds to severe infections that require a hospitalizationbut are not critical. Such infections occur with probability p sev . • The I state which consists of all mild infections that do not lead to a hospitaladmission, and occur with probability − p crit − p sev .Individuals in state C h are admitted to hospital after a duration D C h . Then, withprobability p short they are discharged after a duration D short , while with probability − p short they are discharged after a duration D long .Critically infected individuals are admitted to hospital after a duration D C u . Uponarrival at hospital, they die immediately with probability d hosp , or go to ICU after aduration D H u . Individuals in ICU die with probability d ICU after a delay D D . Otherwisethey are discharged after a stay of length D U .30 R RDUHI i n f e c t i o n h o s p i t a l a d m i ss i o n d e a t h d hosp p sev -p crit age ofinfection HD RR I C U a d m i ss i o n h o s p i t a l a d m i ss i o n d i s c h a r g e d i s c h a r g e D short D Ch D U D U D long D Cu D D d hosp d icu d icu p short p short p crit p sev C u C h Figure 6: Illustration of the occupancy modelIn our model, the probability of hospital admission is p crit + p sev , the probability ofICU admission is p crit (1 − d hosp ) and that of death is p crit ( d hosp + (1 − d hosp ) d ICU ) . Wehave fixed these three values to those estimated in [11], and we only had one remainingparameter out of 4 ( p crit , p sev , d short , d ICU ) to estimate from the data. We have fixed thetime D U to . days as estimated in [11]. All other parameters were estimated using thesame likelihood method as previously, which is described in Section A. The estimatedparameter set is shown in Table 3, while Figure 7 shows the best-fitting model.The estimated parameters provide a good fit of the six observed time series. Again,the model has a tendency to overestimate the ICU admissions in the second part of thelockdown, which has the same interpretation as before.Under the occupancy model, we estimated that R = 0 . , and W = 9 . × . Theseestimates are extremely close to those made with the admission model. The estimatedmean time between infection and death or hospital, ICU admission are respectively . days, . days and . days. Again we see that these estimates in the more complexmodel are consistent with those of the simple model. The mean recovery time fromhospital is . days for severe infections, and . days for critical infections. Thisyields an overall mean recovery time of . days. Finally, we estimated that the deathprobability conditional on being in ICU is . . This yields that in our model a fraction . of all deaths occur shortly after hospital admission. This result is consistent with [11]that estimated that a fraction . of all deaths occurred within the first day after hospitaladmission. However, it has been reported in [33] that the death probability of ICU patientsis . . Our estimated value is thus unrealistically high. This indicates that there is afraction of hospital deaths that occur without any ICU admission, and not quickly afterhospital admission, that our model is not accounting for.31igure 7: Best fit of the admission model. The solid lines correspond to the numberof deaths, discharges, occupied ICU and hospital beds and ICU and hospital admissionspredicted by the occupancy model. The dots are the corresponding observed values.Our estimates, though they are not the key message of the present paper, can neverthe-less draw attention to potential heterogeneities in the infected population. We estimatedthat the mean time between infection and ICU admission is shorter than that betweeninfection and hospital admission. This suggests that the time between infection and se-vere symptom onset is shorter for critical infection, that lead to ICU admission, than formilder ones. Moreover, fitting the prevalence time series required to divide the hospitaland death compartments in two subcompartments, indicating that the data are not wellexplained by a simple homogeneous model, as seen in Figure 8. Such heterogeneity couldoriginate from underlying structuring variables, such as comorbidity or (real) age, that weare not accounting for. Many estimates of clinical features, such as the incubation period,are obtained from a pooled dataset that does not take heterogeneity in the populationinto account [34, 31, 35, 36, 37, 38, 14]. When estimating the total number of infectedindividuals using only a fraction of the detected cases, e.g., using the hospital admissionsor deaths, it is interesting to keep in mind that the time periods estimated from pooledstudies could be inaccurate for the fraction of infected individuals under consideration. Acknowledgements.
The authors thank the
Center for Interdisciplinary Research inBiology (CIRB) for funding and the taskforce MODCOV19 of the INSMI (CNRS) fortheir technical/scientific support during the epidemic. P.C. has received funding fromthe European Union’s Horizon 2020 research and innovation program under the MarieSkłodowska-Curie grant agreement PolyPath 844369.32 otation Description Value Source R Basic reproduction number duringlockdown . E W Total number of infections beforeMarch 17 2020 . × E p crit + p sev Probability of being hospitalized . [11] p crit (1 − d hosp ) p crit + p sev Probability of entering ICU conditionalon being at the hospital . [11] d hosp +(1 − d hosp ) d ICU p sev /p crit Death probability conditional on beinghospitalized . [11] d ICU
Probability of death conditional on be-ing in ICU . E p short Probability of a short stay at hospital . E D C h Delay between severe infection and hos-pital admission . days E D short Delay between hospital admission andquick discharge . days E D long Delay between hospital admission andslow discharge . days E D C u Delay between critical infection andhospital admission . days E D H Delay between hospital admission andICU admission . days [11] D U Delay between ICU admission and dis-charge . days E D D Delay between ICU admission anddeath . days E γ Scale parameter common to all Gammadistributions . ETable 3: Inferred parameter set for the occupancy model. The values indicated for thedurations correspond to the means of the Gamma distributions. In the “Source” column,“E” indicates that the parameter has been estimated in the current work.33 eferences [1] H. Krauss.
Zoonoses : infectious diseases transmissible from animals to humans .ASM Press, 2003.[2] Dorothy H. Crawford.
Deadly companions: how microbes shaped our history . OxfordUniversity Press, second edition, 2018.[3] Joseph T. Wu, Kathy Leung, Mary Bushman, Nishant Kishore, Rene Niehus,Pablo M. de Salazar, Benjamin J. Cowling, Marc Lipsitch, and Gabriel M. Leung.Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan,China.
Nature Medicine , 26(4):506–510, 2020.[4] Olle Nerman. On the convergence of supercritical general (C-M-J) branching pro-cesses.
Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete , 57(3):365–395, 1981.[5] Peter Jagers and Olle Nerman. The growth and composition of branching popula-tions.
Advances in Applied Probability , 16(2):221–259, 1984.[6] Ziad Taïb.
Branching Processes and Neutral Evolution , volume 93 of
Lecture Notesin Biomathematics . Springer Berlin Heidelberg, Berlin, Heidelberg, 1992.[7] Tom Britton, Frank Ball, and Pieter Trapman. The disease-induced herd immunitylevel for COVID-19 is substantially lower than the classical herd immunity level. arXiv preprint arXiv:2005.03085 , 2020.[8] Tom Britton. Epidemics in heterogeneous communities: estimation of R and securevaccination coverage. Journal of the Royal Statistical Society: Series B (StatisticalMethodology) , 63(4):705–715, 2001.[9] Harry Cohn. Multitype finite mean supercritical age-dependent branching processes.
Journal of applied probability , 26(2):398–403, 1989.[10] Jean-Jil Duchamps, Felix Foutel-Rodier, and Emmanuel Schertzer.
From individual-based epidemic models to McKendrick-von Foerster PDE’s: the non-linear case . Inprogress, 2020.[11] Henrik Salje, Cécile Tran Kiem, Noémie Lefrancq, Noémie Courtejoie, Paolo Bosetti,Juliette Paireau, Alessio Andronico, Nathanaël Hozé, Jehanne Richet, Claire-LiseDubost, Yann Le Strat, Justin Lessler, Daniel Levy-Bruhl, Arnaud Fontanet, LullaOpatowski, Pierre-Yves Boelle, and Simon Cauchemez. Estimating the burden ofSARS-CoV-2 in France.
Science , 369(6500):208–211, 2020.[12] Lionel Roques, Etienne K Klein, Julien Papaix, Antoine Sar, and Samuel Soubeyrand.Using early data to estimate the actual infection fatality ratio from COVID-19 inFrance.
Biology , 9(5):97, 2020.[13] Theodoros Evgeniou, Mathilde Fekom, Anton Ovchinnikov, Raphael Porcher,Camille Pouchol, and Nicolas Vayatis. Epidemic models for personalised COVID-19 isolation and exit policies using clinical risk predictions. medRxiv , 2020.3414] Ramses Djidjou-Demasse, Yannis Michalakis, Marc Choisy, Micea T Sofonea, andSamuel Alizon. Optimal COVID-19 epidemic control until vaccine deployment. medRxiv , 2020.[15] Thomas Sellke. On the asymptotic distribution of the size of a stochastic epidemic.
Journal of Applied Probability , 20(2):390–394, 1983.[16] Guodong Pang and Etienne Pardoux. Functional limit theorems for non-Markovianepidemic models. arXiv preprint arXiv:2003.03249 , 2020.[17] Frank Ball. A unified approach to the distribution of total size and total area underthe trajectory of infectives in epidemic models.
Advances in Applied Probability ,18(2):289–310, 1986.[18] Andrew D Barbour. The duration of the closed stochastic epidemic.
Biometrika ,62(2):477–482, 1975.[19] Peter Jagers.
Branching processes with biological applications . Wiley, 1975.[20] Jie Yen Fan, Kais Hamza, Peter Jagers, and Fima Klebaner. Convergence of the agestructure of general schemes of population processes.
Bernoulli , 26(2):893–926, 2020.[21] Peter Jagers and Olle Nerman. Limit theorems for sums determined by branching andother exponentially growing processes.
Stochastic Processes and their Applications ,17(1):47–71, 1984.[22] Peter Jagers and Fima C Klebaner. Population-size-dependent, age-structuredbranching processes linger around their carrying capacity.
Journal of Applied Prob-ability , 48(A):249–260, 2011.[23] Peter Jagers and Fima C Klebaner. Population-size-dependent and age-dependentbranching processes.
Stochastic Processes and their Applications , 87(2):235–254,2000.[24] Kais Hamza, Peter Jagers, and Fima C Klebaner. The age structure of population-dependent general branching processes in environments with a high carrying capacity.
Proceedings of the Steklov Institute of Mathematics , 282(1):90–105, 2013.[25] Viet Chi Tran. Large population limit and time behaviour of a stochastic particlemodel describing an age-structured population.
ESAIM: Probability and Statistics ,12:345–386, 2008.[26] Regis Ferriere and Viet Chi Tran. Stochastic and deterministic models for age-structured populations with genetically variable traits. In
ESAIM: Proceedings , vol-ume 27, pages 289–310. EDP Sciences, 2009.[27] Peter Olofsson. The x log x condition for general branching processes. Journal ofapplied probability , 35(3):537–544, 1998.[28] Emmanuel Schertzer and Florian Simatos. Height and contour processes of Crump-Mode-Jagers forests (I): general distribution and scaling limits in the case of shortedges.
Electronic Journal of Probability , 23:43pp., 2018.3529] Luca Ferretti, Chris Wymant, Michelle Kendall, Lele Zhao, Anel Nurtay, LucieAbeler-Dörner, Michael Parker, David Bonsall, and Christophe Fraser. Quantify-ing SARS-CoV-2 transmission suggests epidemic control with digital contact tracing.
Science , 368(6491), 2020.[30] Mircea T. Sofonea, Bastien Reyné, Baptiste Elie, Ramsès Djidjou-Demasse, Chris-tian Selinger, Yannis Michalakis, and Samuel Alizon. Epidemiological monitoringand control perspectives: application of a parsimonious modelling framework to theCOVID-19 dynamics in France. medRxiv , 2020.[31] Natalie M. Linton, Tetsuro Kobayashi, Yichi Yang, Katsuma Hayashi, Andrei R.Akhmetzhanov, Sung-mok Jung, Baoyin Yuan, Ryo Kinoshita, and Hiroshi Nishiura.Incubation period and other epidemiological characteristics of 2019 novel coronavirusinfections with right truncation: A statistical analysis of publicly available case data.
Journal of Clinical Medicine , 9(2):538, 2020.[32] Robert Verity, Lucy C. Okell, Ilaria Dorigatti, Peter Winskill, Charles Whittaker,Natsuko Imai, Gina Cuomo-Dannenburg, Hayley Thompson, Patrick G. T. Walker,Han Fu, Amy Dighe, Jamie T. Griffin, Marc Baguelin, Sangeeta Bhatia, AdhirathaBoonyasiri, Anne Cori, Zulma Cucunubá, Rich FitzJohn, Katy Gaythorpe, WillGreen, Arran Hamlet, Wes Hinsley, Daniel Laydon, Gemma Nedjati-Gilani, StevenRiley, Sabine van Elsland, Erik Volz, Haowei Wang, Yuanrong Wang, Xiaoyue Xi,Christl A. Donnelly, Azra C. Ghani, and Neil M. Ferguson. Estimates of the severityof coronavirus disease 2019: a model-based analysis.
The Lancet Infectious Diseases ,20(6):669–677, 2020.[33] Santé Publique France. COVID-19 : point épidémiologique du 4 juin 2020, 2020.[34] Jantien A. Backer, Don Klinkenberg, and Jacco Wallinga. Incubation period of 2019novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China, 20-28January 2020.
Eurosurveillance , 25(5), 2020.[35] Stephen A Lauer, Kyra H Grantz, Qifang Bi, Forrest K Jones, Qulu Zheng, Hannah RMeredith, Andrew S Azman, Nicholas G Reich, and Justin Lessler. The incubationperiod of coronavirus disease 2019 (COVID-19) from publicly reported confirmedcases: Estimation and application.
Annals of internal medicine , 172(9):577–582,2020.[36] Lauren Tindale, Michelle Coombe, Jessica E Stockdale, Emma Garlock, WingYin Venus Lau, Manu Saraswat, Yen-Hsiang Brian Lee, Louxin Zhang, DongxuanChen, Jacco Wallinga, and Caroline Colijn. Transmission interval estimates suggestpre-symptomatic spread of COVID-19. medRxiv , page 2020.03.03.20029983, 2020.[37] Qifang Bi, Yongsheng Wu, Shujiang Mei, Chenfei Ye, Xuan Zou, Zhen Zhang, Xiao-jian Liu, Lan Wei, Shaun A Truelove, Tong Zhang, Wei Gao, Cong Cheng, XiujuanTang, Xiaoliang Wu, Yu Wu, Binbin Sun, Suli Huang, Yu Sun, Juncen Zhang, TingMa, Justin Lessler, and Teijian Feng. Epidemiology and transmission of COVID-19in Shenzhen China: Analysis of 391 cases and 1,286 of their close contacts. medRxiv ,2020. 3638] Clément Massonnaud, Jonathan Roux, and Pascal Crépey. COVID-19: Forecastingshort term hospital needs in France. medRxiv , 2020.[39] Santé Publique France. Données hospitalières relatives à l’épidémie de COVID-19,2020.
A Maximum likelihood method