Branching Processes: Their Role in Epidemiology

Branching processes are stochastic individual-based processes leading consequently to a bottom-up approach. In addition, since the state variables are random integer variables (representing population sizes), the extinction occurs at random finite time on the extinction set, thus leading to fine and realistic predictions. Starting from the simplest and well-known single-type Bienaymé-Galton-Watson branching process that was used by several authors for approximating the beginning of an epidemic, we then present a general branching model with age and population dependent individual transitions. However contrary to the classical Bienaymé-Galton-Watson or asymptotically Bienaymé-Galton-Watson setting, where the asymptotic behavior of the process, as time tends to infinity, is well understood, the asymptotic behavior of this general process is a new question. Here we give some solutions for dealing with this problem depending on whether the initial population size is large or small, and whether the disease is rare or non-rare when the initial population size is large.


Introduction
Mathematical models of propagation of a disease in given populations play a central role for understanding this propagation, for predicting the future extension of the outbreak, its extinction time, and for evaluating the efficiency of control measures. Of course the validity and the richness of results of a model strongly depend on the reliability and the accuracy of the model. A fine predictive model should be built as far as possible in a rigorous mechanistic way starting from the mechanism of exposure/infection of each individual and taking into account their variability. The population dynamic described by births, deaths, migrations should also be taken into account, especially when the incubation time is relatively long in comparison with this dynamic. Of course, the time unit of the model should be chosen in keeping with the respective durations of each health state and the data.
Nevertheless, since the asymptotic behavior, as time tends to infinity, of deterministic models in continuous time are more easily studied than that of discrete time models or stochastic models, a large literature in applied mathematical journals is devoted to such theoretical studies [19], while for statistical purposes, either descriptive non dynamical models or very simple stochastic models are used, such as a chain binomial model for the propagation of SIS or SIR diseases in closed populations [5,6], or BGW (Bienaymé-Galton-Watson) branching processes on the clinical (or diagnosed) cases as an approximate model of the beginning of an outbreak [2,3,4,6,10,11,20]. These two approaches are both individualbased (each individual transition is given) and stochastic (individual variability is taken into account) but the first one takes into account the individual health state evolution S → I → S or S → I → R, where S means susceptible, I means infective or clinical cases (according to authors and purposes) and R means removed from the susceptible population, either by immunization or death; while the second approach directly models the process of I individuals, without analyzing this mechanism of exposure/infection. The behavior of these two types of processes are well-known. The chain binomial models are just Markov chains in a finite state space, and since there is no immigration and no incubation period, the state "0 infected individual" is an absorbing state. Concerning BGW processes, a huge literature exists for describing their properties. In fact these two different types of models may be written as particular cases of a large class of branching processes with individual transitions that may be population and age dependent, and which can take into account both the population dynamics and the disease dynamic.
Branching processes were initiated in the nineteen century by Bienaymé and then by Galton and Watson, for studying the extinction of some family names. Since this time, the complexity of these processes continues to increase allowing to describe more and more realistic population dynamics. These processes are based on the simple property that the population size of each considered type (such as clinical cases here) at each time is calculated as the sum of all the new individuals ("offspring") of this type who are generated by the individuals of the population at the previous times. Since the modelled variables are integers, then the extinction time of the population of each type is finite on the set of trajectories which extinct. This is a finer and more realistic property than the asymptotic extinction time given by a deterministic model. Since the population dynamic may influence the disease dynamic, and conversely, these two dynamics should be explicitly taken into account in the model, thus leading to rigorous, but not simple, multitype models where each type should represent an health state crossed with influence factors levels. Typical examples of such influence factors are age, geographical locations. However, some authors model directly the time evolution of the incidence of clinical cases, without explicating all the intermediary steps.
In the following subsections, we present such direct models starting from the simplest one, the singletype Bienaymé-Galton-Watson process, and finishing by a general and rigorous approach that takes into account the intermediary steps and the population dynamic, and is based on age-dependent and population-dependent individual transitions. We focus on models in discrete time since they have the double advantage to be easily written as recursive models and to easily correspond to the time unit of observation, which offers a pedagogic framework and moreover facilitates the model validation and the estimation of unknown parameters. We study here the behavior of these models. From now on, all results are given conditionally to the initial value F 0 of the process and the notation "|F 0 " will be therefore omitted for the sake of simplicity of formulas. The proofs are given in detail in [23] and will be omitted here.

Single-type Bienaymé-Galton-Watson Branching Processes
Let I n be the incidence of clinical (or diagnosed) cases at time n and Y n,i represent the numbers of secondary cases generated at time n by the previous case i. Then I n is modelled from the past incidences by where the variables {Y n,i } i are assumed i.i.d. (independently and identically distributed) according to a distribution L independent of the time and of the past process given F n−1 := {I k } k≤n−1 . We denote m, σ 2 the mean and variance of Y n,1 given F n−1 . Since these first two moments are the moments influencing the behavior of the process on the non-extinction set, we will also write L(m, σ 2 ) rather than L. This comes from the following writing of (1): is asymptotically distributed according to N (0, σ 2 ) on {lim n I n = ∞}, which is the non-extinction set, as n → ∞.
The behavior of this process has been deeply analysed for a long time (see for example [1,13,16,24]). We recall here the main properties. The trajectories of this process either die out or explode in an exponential way: P ({lim n I n = ∞} ∪ {lim n I n = 0}) = 1 and there exists an integrable random variable W such that lim n I n m −n = W , a.s. (almost surely), where P (W > 0) > 0 if and only if m > 1 (supercritical case). The random variable W depends on I 0 and of L(m, σ 2 ). This limit behavior comes from the martingale property of I n m −n , that is E(I n m −n |F n−1 ) = I n−1 m −(n−1) which implies E(I n m −n ) = I 0 . So this process reproduces the initial phase of exponential growth of an epidemic and can be used for describing this phase when the incubation period is negligible compared to the time unit. The quantity m := E(Y n,1 |F n−1 ) is the current reproductive number of the process (mean number of secondary cases produced by one case during a time unit) and is known to be the bifurcation parameter of the process, that is, m ≤ 1 implies the a.s. extinction of the process. The deterministic model X n = mX n−1 with X 0 = I 0 , derived from E(I n |F n−1 ) = I n−1 m, has the same bifurcation parameter, but when m = 1, the deterministic model persists with X n = I 0 while the stochastic one dies out a.s.. Moreover before dying out, the stochastic process I n |I n ̸ = 0 presents a linear increasing behavior: lim n P (I n [0.5σ 2 n] −1 ≤ x|I n ̸ = 0) = 1 − exp(−x).
In addition, the probability of extinction, q, is the solution of q = f (q) := E(q Y 1,1 ) and in the subcritical/critical cases m ≤ 1, the "epidemic size" N := ∑ Text. k=0 I k (total number of cases generated until the extinction time T ext. ), follows a power series distribution when Y n,1 follows itself a power series distribution, that is, P (Y n,i = k) ∝ a k λ k [7,11,12]. In particular, when P (Y n,i = k) = exp(−m)[k!] −1 m k (P oisson(m) distribution), then N follows the Borel − T anner(m, I 0 ) distribution, for m < 1, and

Single-type Branching Processes with Population-Dependent Offsprings
When the population is small or when the individual migrations are slow compared to the infection process, the depletion of susceptible individuals due to the infection should be taken into account, which is not the case in the BGW process. The typical bell curve form of outbreaks is due to this depletion. The simplest such model is just an extension of the BGW process, that is, The depletion effect of the S population at time n is replaced by the fact that m(F n−1 ) (resp. P (Y n,i = 0|F n−1 )) is a decreasing (resp. increasing) function of the previous incidences of cases We will see in this section that this kind of models may exhibit a random cycle and therefore may be used for modeling recurrent epidemics. The extremal case µ = 0 means no depletion in S, which is got by an immediate healing with no immunization, after one time unit in the state I (SIS disease) and this model is reduced to the previous BGW process.
The other extremal case µ = 1 with d n = n corresponds either to a not necessarily immediate healing with persistent immunization or to a fatal issue (SIR disease). This process is a branching process with a long memory and has not be studied in the literature excepted by simulation [27].
The intermediate case 0 < µ ≤ 1 with d n = d represent the possibility of recovering but with a transient immunization only, and leads to a Markovian chain of order d. The behavior of this process has been studied in [27] and is described in the following paragraphs.
Let A1: there exists a positive function f (.) such that P (Y n,i = 0|I n−1 = N, F n−1 ) ≥ f (N ) > 0, for all N ∈ N + and any value of F n−1 consistent with This assumption is checked as soon as P (Y n,i = 0|F n−1 ) is a non-decreasing function of I n−1 , I n−2 ,. . . , which is strictly increasing in I n−1 .
Let us define the bifurcation parameter R ∞ by where R is any real quantity depending on the parameters of the distribution of {I n }.
Remark 1 When d n = d, for all n, {I n } may be represented as a multitype branching processĨ n = (I n , I n−1 , ..., I n−(d−1) ) =: (I n,1 , I n,2 , ..., I n,d ), which is asymptotically BGW if L(F n−1 ) tends to a distribution L independent of the process as time tends to ∞ on {lim n I n = ∞}. But the asymptotic BGW process is not positively regular at the opposite of classical BGW branching processes, sinceM Moreover since R 0 = α > 1, then {X n } tends to a stable limit cycle, as n → ∞ [9,27].

Set of Single-type Branching Processes with Population-Dependent Offsprings
We generalize here the model of the previous section to a set of J similar diseases in competition: A simple example is when as previously n−l . A typical example is given by different influenza viruses. As in the single-type case, assuming m (j) (F ) < m * < ∞, for all F , with R (j) Therefore as soon as the number of cases due to any pathogenic agent decreases or dies out, then the epidemics due to the other pathogenic agents may increase (see Figure 1).
Let the deterministic associate trajectory {X A possible extension of the previous models consists in adding an immigration. Then: where J n is an immigration at time n when δ n = 1, and δ n is a Bernoulli variable allowing or not an immigration at time n. When {δ n } follows some seasonality, then the model is suitable for recurrent seasonal epidemics with seasonal immigrant, such as influenza. This type of processes has been deeply studied in the subcritical case m < 1, when ({Y n,i } i , J n ) has a distribution independent of n and F n−1 , with either an immigration allowed only in the periods of extinction of the process {I n } (regenerative branching process), or when the {δ n } are i.i.d. independent of F n−1 (see the review article [39] and the estimation point of view in [25]).
We may also write (3) in the following way: Consequently, when ({Y n,i } i , J n , δ n ) is not time-dependent and when 1 {I n−1 =0} J n δ n = 0, for all n, then (3) may be written as a model of Section 3..

Multitype Branching Processes with Age and Population-Dependent Offsprings
where h, k represent types and l is the maturation time for getting the "offsprings" Y Each line of graphics concerns a trajectory of the process {I (1) n , I (2) n }, and on each line, the graphic on the left concerns population 1 and the graphic on the right, population 2. On each graphic, the red line represents the deterministic limit cycle (reached very quickly) and the blue one, the stochastic cycle. We see that when one population is small during a long enough period, then the other population may be large, but both populations may also die out very quickly (second line). : bl ue, X n ( 2) : r ed should be expressed according to the number of newborn individuals of i, who are of the k type, and to the proper transition of i from the state h at n − l to the state k at n. For example assuming that there is no influence factors and that the set of individual health states is {S, E, I} and assuming that the life duration in the state I is at most one time unit, then, for k = I,   [30,31,32,33], we cannot deduce from this martingale convergence any information on the asymptotic behavior of {N n } itself.
In Section 5.1., we give some approaches for studying the behavior of model (4) when the total population size remains bounded. In Sections 5.2. and 5.3., we study the asymptotic behavior of limit models derived from (4) as the initial population size increases to infinity.

The Population Size Remains Bounded
Let N n := ∑ k N k n be the total population size at time n. Let us first assume that this population size remains bounded by some control. For example the population is a herd of farm animals. We assume that the number of newborns Y n,i at time n of the animal i is bounded whatever i, n, and we use the following control: if N n−1 > N M , for some chosen threshold N M , then Y n,i = 0 (all the newborn animals are sold), and when N n−1 < N m , for some chosen threshold N m < N M , then new animals are bought. Therefore assuming that P (N n = N 2 |N n−1 = N 1 ) =: Q (N 1 , N 2 ) is independent of n, {N n } is a homogeneous Markov chain on a finite space. Since the population is open to immigration, then the healthy state "0 case" is not an absorbing state and there may exist some endemic behavior: if there exists an asymptotic distribution π, it satisfies ∑ But due to combinatorial aspects, the transition probabilities {Q (N 1 , N 2 )} may be difficult to compute. A solution is then to work in continuous time: in [26], we determine the distribution of the population process from the individual transitions that may be age and population-dependent and that may be all different, and in [40], we use this type of process for studying the propagation of the BVD (Bovine Diarrhoea Virus) within a dairy herd. For that, a renewal approach is used, which generalizes to a population the semi-Markov process theory relative to one individual. This approach also allows to build a rigorous and general simulation algorithm for this individual based model, but when used for non bounded branching populations, it cannot lead to fine stochastic behaviors such as lim n N n [N 0 ρ n ] −1 a.s. = W determined by martingale theory in the frame of BGW processes.

The Initial Population Size N 0 is Large and the Disease is not Rare
We assume here that the initial population size N 0 is large which allows to study the limit, as N 0 → ∞, of the following quantities: N n /N 0 =: D n (empirical densities) or N n /N n =: P n (empirical probabilities) or N n [N 0 ρ n ] −1 =: W n (empirical normalized densities), when assuming that lim N 0 D 0 = D 0 or lim N 0 P 0 = P 0 .
A simple theoretical example that we (artificially) apply here on epidemics is the following model on densities with D = 1, d = 1 [34,35,36,37,38]: given F n−1 = {I n−1 , ..., I 0 } and depend on the previous cases incidences only through the condition: Y n,i = 0 when I n−1 > cI 0 , that is, a massive vaccination is done as soon as the density I n−1 I −1 0 crosses the threshold c. This implies that the process dies out a.s.. The author studied D n := I n /I 0 and proved that lim I 0 lim n D n | D n ̸ = 0 belongs to the set of stable fixed points of D n , where D n is the deterministic limit model D n = D n−1 m(D n−1 ), n ≥ 1, D 0 = 1, with m(D n−1 ) := E(Y n,1 | D n−1 = D n−1 ). Thus, for I 0 very large (which occurs when the initial time is chosen when the epidemic is large enough), the random empirical densities has the same asymptotic behaviour, as n → ∞, as the limit densities.
Let us study now the empirical probabilities in the general case D ≥ 1, d ≥ 1 under a densitydependence assumption. We prove in Proposition 5 that lim n lim N 0 P n |N n ̸ = 0 P = lim N 0 lim n P n |N n ̸ = 0 (and similarly for {W n }), which allows to approximate lim n P n or lim n W n .
Let us first study the behavior of the total population size, under the assumption that the number of newborns is independent of the mother health state. n,i } i are independent given F n−1 .

total mean number of offspring generated by an individual).
Let us notice that the process {N n } is the counterpart in discrete time of a single-type age-dependent Crump-Mode-Jagers branching process in continuous-time [29].
Let us write α l := Φ l ρ −l , l = 1, ..., d, d n := d ∧ n, n d = ⌊n/d⌋ (integer part of n/d) and let us define: For d = 1, or d > 1 with ϕ l = 0, for all l < d, then the sum in A3 is reduced to 0 implying that A3 is satisfied. In the general case A3 is satisfied under the stronger assumption: where, we notice that, for all k ≤ n, α k 1 ≤ is decreasing in k.

Proposition 5 Let us assume, as in Proposition 4, that the distribution of
where P n is the vector of probabilities defined by: Proposition 5 allows to use the attractors of the dynamical model {P n } as an approximation of the asymptotic behavior of the empirical frequencies if the initial population size is very large, and moreover generalizes the result of the BGW process: = W (Section 2.). But the main drawbacks of this approach consist in the loss of the population variability, since the quantities are normalized by N 0 with N 0 tending to infinity, and the loss of the possibility for the extinction time of any type to be finite. Moreover the global asymptotic stability of the healthy state for {P n } (extinction of the disease propagation starting from any initial infected population size) may be difficult to prove, the difficulty increasing with the number d × D of dimensions. This global asymptotic stability was studied for the propagation of a SIS disease in a branching population (D = 2) [21] and for the one of a SEI disease (D = 3) [22]. Both studies assumed d = 1 and we showed that the bifurcation parameter for the disease process was based on the comparison of the capacity of infection and the capacity for the population to renew its susceptible population.

The Initial Population Size is Large and the Disease is Rare
Another approach for studying the asymptotic behavior of the process {N n } is to study the asymptotic behavior of the limit model assuming now that the disease is rare at the initial time, which forbids the use of densities or probabilities as in Section 5.2. We present here such an approach generalizing the case of the BSE propagation at the scale of a country [28]. Let us recall model (4): 1,n−l+1,n,i,j } i,j are mutually independent given F n−1 . Therefore the incidence of cases aged a at time n is given by: where N (h) k,i . Moreover, if a I individual may survive in this state a longer time than a time unit, the infection process will depend on the total number of infectives including all the I individuals. In this case, let k = I tot × a, where I tot represents the clinical state (new or not), then (6) where δ I,I tot |I×a−l n−l,n,i is a Bernoulli variable, equal to 1 if i survives in the state I from age a − l at n − l to age a at n at least. Since the distribution of {N I tot ×a n } a,n is easily calculated from the distribution of {N I×a n } a,n and since clinical cases are generally rapidly isolated and then removed from the infection process, and since moreover only the new cases are counted by surveillance systems, we will study here the limit of N I×a n given F * n−1 : Y,k . Of course, expressions similar to (5) and (6) where Ψ 1,a,l|n−l and Ψ 2,a,l|n−l are independent of F * n−l . Then the limit process I a,n : D = lim N 0 N I×a n of incidence of cases aged a at time n is a single-type Markovian process of order d with a Poissonian transition law: Ψ a,l|n−l := 1 {a>l} Ψ 2,a,l|n−l + 1 {a=l} Ψ 1,a,l|n−l , where I n := ∑ a I a,n . Moreover we may write I n = As an example of such an approach, we studied the propagation of the BSE in Great-Britain by a general model of type (4). Since the time unit is large (one year), we used Proposition 6 with, for l ≥ 2,F * n−l : Then the assumptions of this proposition 6 were satisfied under the following assumptions [28] 5. The disease is rare at the initial time: 6. The probability for a given S to be infected at time k + 1 via the horizontal route of excretion, follows a Reed-Frost's type model, that is P (i aged a at k , is inf ected at k + 1 by excretion|F * k , i survives at k + 1)= is the number of infectious animals at time k (including those in the latest stage of their incubation period). We assume a similar expression for the infection via the horizontal route concerning contaminated meat and bone meal produced from dead infectious animals.
Under these assumptions, then where θ a−l,R c l and θ a−l,R l (exp λ − 1) −1 ϕ n−l represent the mean numbers of new infected animals aged a − l + 1 per infective at time n − l + 1 respectively via a horizontal route (excretion of alive animals and slaughtered animals respectively), ϕ n−l ∈ [0, 1] represents the efficiency at time n − l of the main current control regulations: ϕ k = 1, for k ≤ 1988, ϕ k = ϕ, for k ∈ (1989,1996), and ϕ k = 0, for k ≥ 1997, where ϕ represents the efficiency of the Meat and Bone Meal ban of July 1988, p mat. is the probability for a newborn animal to be infected by his mother (maternal route), P (a) is the probability at each time for an animal to have age a which may be expressed as a function of the survival distribution, and P inc. (l) is the probability for an infected animal to have an intrinsic incubation time equal to l (conditioned on survival). So Ψ l|n−l represents the mean number of new clinical cases produced with a delay l by a clinical case of time n − l.
Let us assume that Ψ l|n−l depends only on l. This is achieved by modeling the process on a period with a constant control. Then {I n } may be written as a multitype BGW branching process and therefore results of Proposition 4 are valid for this process, replacing Φ l by Ψ l := Ψ l|n−l . More precisely, let I n =: (I n,1 , I n,2 , ..., I n,d ) := (I n , I n−1 , ..., I n−(d−1) ) and F n−1 = {I n−1 , I n−2 , ..., I 0 }.

Proposition 12
For any value of R ∞ , the extinction time distribution satisfies, for n ≥ d, ) .
Let us notice that G n decreases (and R ∞ increases) as soon as at least one Ψ l increases. In addition, using Proposition 8, we get the following result, allowing an exact iterative computation of {G n }: Proposition 13 For any value of R ∞ , the distribution {G n } of T ext. satisfies: 2. LetĨ 0 := (I 0 , I −1 , , ..., and

Discussion
We presented a general class of branching processes in discrete time for modeling in a stochastic way some diseases propagation when the infected period is long respectively to the time frequency of births. However when the transitions are population dependent, the long-term prediction of these processes is an open problem in the general case. We indirectly solved this problem by studying the behavior of the limit models, as the total initial population size increases to infinity, assuming at the initial time, either a non-rare disease with a density-dependence assumption, or a rare disease.
Under the first assumption, since lim n lim N 0 P n |N n ̸ = 0 P = lim N 0 lim n P n |N n ̸ = 0, we proved, for a large N 0 , that the proportion of infected individuals in the whole population, P E n + P I n , behaves, as n → ∞, as the probability P E n + P I n for an individual in an infinite population to be infected. Under the second assumption, we got a branching process on the incidence of cases which has the advantage to correspond to the usual observations, to be rigorously built, starting from a detailed multitype branching process {N n } taking into account the different disease steps together with the population dynamic, to keep the population variability, and to belong to the simple class of multitype BGW processes for which many analytical results exist and that we used or generalized to this model. We calculated the distribution of the extinction time and the distribution of the epidemic size. Moreover we proved that the bifurcation parameter of this process was the total mean number of secondary cases who began to be generated by one case during his first time unit.
This result would validate and generalize the current use in epidemiology of the reproductive number as a bifurcation parameter [8,17,18], if we could establish that lim n lim N 0 N I n = lim N 0 lim n N I n . We studied the left quantity, while the valid realistic quantity is in reality the right one. In fact this equality is in general not true because the bell form of an epidemic can be modelled only by a population dependent process. But for large N 0 , the stochastic limit model is a good approximation of the epidemic growth (in the supercritical case), or of the epidemic decay (in the subcritical case). In any case, if N 0 is too small, the limit in N 0 cannot be used and therefore the limit models developed here cannot be used.
Let us finally notice that we proved that in a size-dependent model on the clinical cases, then the quantity determining the extinction of the process was the total mean number of secondary cases that will be produced in the future by a case as the whole current population is infected (R ∞ ), while the corresponding quantity in the associate deterministic model derived from the conditional expectation of the process is the total mean number of secondary cases that will be produced by a case as the whole current population is susceptible (R 0 ), which easily leads to the extinction of the process on one hand and the persistence of the deterministic trajectory on the other hand. However, simulations of the single-type process with population-dependent offsprings described in Section 3., showed that until its extinction, the process roughly behaved as its deterministic counterpart and the extinction time strongly depends on the parameters of R 0 . The extinction time roughly increases as R 0 increases. So the greatest difference between the behavior of this population-dependent process and its deterministic counterpart is obtained when R 0 > 1 with R 0 ≃ 1. This is the generalization of the difference observed between the BGW process and its deterministic counterpart, when R 0 (= ρ)= 1.