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 individual-based (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 single-type 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 F0 of the process and the notation “|F0” 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 be the incidence of clinical (or diagnosed) cases at time n and Y represent the numbers of secondary cases generated at time n by the previous case i. Then I is modelled from the past incidences by where the variables {Y} are assumed i.i.d. (independently and identically distributed) according to a distribution independent of the time and of the past process given F1 := {I}1. We denote m, σ2 the mean and variance of Y1 given F1. Since these first two moments are the moments influencing the behavior of the process on the non-extinction set, we will also write rather than . This comes from the following writing of (1): , where is asymptotically distributed according to on {lim I = ∞}, 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 I = ∞} ∪ {lim I = 0}) = 1 and there exists an integrable random variable W such that lim I− = W, a.s. (almost surely), where P(W > 0) > 0 if and only if m > 1 (supercritical case). The random variable W depends on I0 and of . This limit behavior comes from the martingale property of I−, that is E(I−|F) = Im−( which implies E(I−) = I0. 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|F) 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 = mX with X0 = I0, derived from E(I) = Im, has the same bifurcation parameter, but when m = 1, the deterministic model persists with X = I0 while the stochastic one dies out a.s.. Moreover before dying out, the stochastic process I ≠ 0 presents a linear increasing behavior: lim P(I[0.5σ2n]−1 ≤ x|I ≠ 0) = 1 − exp(−x).

In addition, the probability of extinction, q, is the solution of q = f(q) := E(q) and in the subcritical/critical cases m ≤ 1, the “epidemic size” (total number of cases generated until the extinction time T), follows a power series distribution when Y follows itself a power series distribution, that is, P(Y = k) ∝ a [7, 11, 12]. In particular, when P(Y = k) = exp(−m)[k!]−1m (Poisson(m) distribution), then N follows the Borel − Tanner(m, I0) distribution, for m < 1, and E(N) = I0(1 − m)−1, Var(N) = I0m(1 − m)−3 are increasing functions of m, I0.

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, , where the {Y} are assumed i.i.d. according to a distribution given F : {I}. The depletion effect of the S population at time n is replaced by the fact that m(F) (resp. P(Y = 0|F)) is a decreasing (resp. increasing) function of the previous incidences of cases {I}.

A simple example is when , μ ∈ [0, 1], where is a decreasing (resp. increasing) function of . For example, E(Y|Y ≠ 0) = α and , 0 < μ ≤ 1, K > 0. 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 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 = 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 = 0|I = N, F) ≥ f(N) > 0, for all N ∈ ℕ+ and any value of F consistent with I = N.

This assumption is checked as soon as P(Y = 0|F) is a non-decreasing function of I, I,. . ., which is strictly increasing in I.

Proposition 1 Let us assume A1. Then P(lim I = 0 ∪ lim I = ∞) = 1.

Let A2: there exists m < ∞ such that m(F) ≤ m, for all F, and limm(F) < 1, where |F| is the L1-norm of F.

Let us define the bifurcation parameter R by where R is any real quantity depending on the parameters of the distribution of {I}.

Proposition 2 The bifurcation parameter R∞ is equal to limm(F). Let us assume A1 and A2. Then R∞ < 1, i.e., the process dies out a.s.. Moreover for the deterministic trajectory {X} defined by , where , then the bifurcation parameter is the reproductive number R0 = limm(F), that is, if R0 < 1, then lim X = 0, while if R0 > 1, then lim ≠ 0.

Remark 1 When d = d, for all n, {I} may be represented as a multitype branching process Ĩ = (I1, ..., I(1)) =: (I1, I2, ..., I), which is asymptotically BGW if tends to a distribution independent of the process as time tends to ∞ on {lim I = ∞}. But the asymptotic BGW process is not positively regular at the opposite of classical BGW branching processes, since M͂(F) being defined by E(Ĩ) =: ĨM͂(F), then which implies that, for all n ≥ d − 1, M͂[1, j] = m, and M͂[i, j] = 0, i > 1, where m := lim m(F) and M͂ := lim M͂(F).

Remark 2 If m((1, 0, ..., 0)) > 1 (supercritical assumption at the beginning of the disease spread), since m(F) < 1 under A2, for |F| and n sufficiently large, then we may observe oscillations until the a.s. extinction, the process increasing when |Ĩ|is small enough, and decreasing when it is too large. It is the case of the logistic model , when d = d with d > 1 and m((1, 0, ..., 0)) = α(1 + μK−1)−1 > 1. Moreover since R0 = α > 1, then {X} 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: where, for each j, the are i.i.d. according to given . We assume that m(F) is non increasing in , for each j, each k ≤ n − 1.

A simple example is when as previously , μ ∈ [0, 1], and is a decreasing (resp. increasing) function of . A typical example is given by different influenza viruses.

As in the single-type case, assuming m((F) < m < ∞, for all F, with , then, for all j, , that is each disease dies out. 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).

Figure 1.

Two populations of infectives from similar diseases in competition following the same logistic Poisson model , , j = 1, 2, with , K = 105, μ = 1, d = 20. Each line of graphics concerns a trajectory of the process , 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).

Let the deterministic associate trajectory be defined by , n≥ 1, .

Proposition 3 The reproductive number is the bifurcation parameter for , that is, if , then , and the reproductive number is the bifurcation parameter of , that is if , then . Moreover if , then .

A possible extension of the previous models consists in adding an immigration. Then: where J is an immigration at time n when δ = 1, and δ is a Bernoulli variable allowing or not an immigration at time n. When {δ} 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}) has a distribution independent of n and F, with either an immigration allowed only in the periods of extinction of the process {I} (regenerative branching process), or when the {δ} are i.i.d. independent of F (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}) is not time-dependent and when 1{J = 0, for all n, then (3) may be written as a model of Section 3..

Multitype Branching Processes with Age and Population-Dependent Offsprings

Let satisfying where h, k represent types and l is the maturation time for getting the “offsprings” of the individual i (number of k individuals at time n generated by i belonging to the type h at time n − l). We assume that the are mutually independent given F = {N, ..., N0} and that the are i.i.d. given F.

The models of the previous sections (Sections 2. and 4.) belong to this class when d = d.

When rigorously modeling the propagation of a disease in a given population, taking explicitly into account the population dynamic and the disease dynamic, then a type should be a health state (S, E (incubation) or I, for example) crossed with an influence factors level, and the “offsprings” 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, where the are i.i.d. Bernoulli variables given F, with equal to 1 if the individual i undergoes transitions S → E, at time n − l + 1, and E → I, at time n, which means that his incubation period is l, the are similar i.i.d. Bernoulli variables given F, relative to the newborn j of i, and the are i.i.d variables given F−1, being the number of newborns of i at time n − l + 1. So is the number of I individuals generated at time n with an incubation time l from i who is in the h type at n − l.

Let us write ℕ := (N, N, ..., N). Then E(ℕ) = ℕ (F), where I is the D × D identity matrix, and , for h = 1, ..., D, k = 1, ..., D. Then, assuming that (F) is invertible, for any value F of F, ℕ[ (F0)... (F)]−1 is a martingale with E(ℕ [ (F0)... (F1)]−1|F0) = ℕ0, implying by the convergence theorem for martingale [15] that lim ℕ [ (F0)... (F1)]−1 = WU, for any vector , where W is an integrable random variable. But, except in the simple cases (F1) = (BGW process) or with lim (N) = M (asymptotically BGW process) [30, 31, 32, 33], we cannot deduce from this martingale convergence any information on the asymptotic behavior of {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 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 at time n of the animal i is bounded whatever i, n, and we use the following control: if N1 > N, for some chosen threshold N, then Y = 0 (all the newborn animals are sold), and when N1 < N, for some chosen threshold N < N, then new animals are bought. Therefore assuming that P(ℕ = N2|ℕ = N1) =: Q(N1, N2) is independent of 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(N1, N2)} 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 determined by martingale theory in the frame of BGW processes.

The Initial Population Size N0 is Large and the Disease is not Rare

We assume here that the initial population size N0 is large which allows to study the limit, as N0 → ∞, of the following quantities: N0 =: D̂ (empirical densities) or N =: P̂ (empirical probabilities) or N[N0ρ]−1 =: W (empirical normalized densities), when assuming that lim D̂0 = D0 or lim P̂0 = P0.

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]: , where the {Y} are i.i.d. given F = {I, ..., I0} and depend on the previous cases incidences only through the condition: Y = 0 when I > cI0, that is, a massive vaccination is done as soon as the density crosses the threshold c. This implies that the process dies out a.s.. The author studied D̂ := I/I0 and proved that lim lim D̂/D̂ ≠ 0 belongs to the set of stable fixed points of D, where D is the deterministic limit model D = D1m(D1), n ≥ 1, D0 = 1, with m(D1) := E(Y1|D̂1 = D1). Thus, for I0 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 density-dependence assumption. We prove in Proposition 5 that (and similarly for {W}), which allows to approximate lim P̂ or lim W.

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.

Proposition 4 Let us assume that the distribution of only depends on l and is denoted . Let us denote . Then the total size of the alive population at time n, , may be expressed as a multitype branching process Ñ := (N1, ...N) := (N1, ..., N(1)), where and with, for j = 1, , for j > 1, and , for l ≠ j − 1, and the are independent given F.

Moreover, assuming 0 < Φ < ∞, l = 1, …, d, , and P(Y ≥ 2|F1) > 0 for some l, then {N} satisfies the following properties: Let us notice that the process {N} is the counterpart in discrete time of a single-type age-dependent Crump-Mode-Jagers branching process in continuous-time [29].

{N} is positively regular, nonsingular, and checks the xlogx property;

P(lim N = 0 ∪ lim N = ∞) = 1

Let ρ be the first eigenvalue of M͂ defined by E(Ñ|Ñ−1) =: Ñ−1M͂, that is, Then E(Ñ) = Ñ0ρ, where M͂ = ρ ≤ 1 (subcritical and critical cases) implies that P(lim N = 0) = 1 (a.s. extinction), and ρ > 1 (supercritical case) implies the existence of an integrable random variable W such that , where W := N[N0ρ]−1, and P(lim N = ∞) = P(W > 0). Finally let us assume that N− = ρ−lN0, l = 1, ..., d − 1. Then if ρ ∈ , .

ρ is solution of , and ρ ≤ 1 ⇔R∞ ≤ 1, where (total mean number of offspring generated by an individual).

Let us write α := Φ, l = 1, ..., d, d := d ∧ n, n = ⌊n/d⌋ (integer part of n/d) and let us define:

A3: .

For d = 1, or d > 1 with φ = 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, and is decreasing in k.

Proposition 5 Let us assume, as in Proposition 4, that the distribution of depends only on l. Let us moreover assume ρ > 1, A3, M(F) = M(P̂1) (density-dependence), limP̂0 = P0, and N− = c0, l = 1, ..., d − 1, where c is independent of N0. Then, for all η > 0, lim P(|P̂ − P > η|N ≠ 0)=0 and where P is the vector of probabilities defined by:

Proposition 5 allows to use the attractors of the dynamical model {P} 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: deduced from (Section 2.). But the main drawbacks of this approach consist in the loss of the population variability, since the quantities are normalized by N0 with N0 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} (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} 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): , where the are mutually independent given F1 = {N1, ..., N0} and the are i.i.d. given F1. We assume here that h and k are health states H and K crossed with age a of the individual. For k = I × a, where I represents here the first time unit in the clinical state, then the number of new clinical cases aged a at time n generated by an individual i with a delay of l time units, is thus defined: where we assume that the are i.i.d. Bernoulli variables given F, the are i.i.d. Bernoulli variables given F, the are i.i.d. given F1, and the and the are mutually independent given F. Therefore the incidence of cases aged a at time n is given by: where .

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 × a, where I represents the clinical state (new or not), then where 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 is easily calculated from the distribution of 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 given , where . Of course, expressions similar to (5) and (6) can be written for and .

Let δ(i) be the Bernoulli variable equal to 1 if i belongs to the h population at time n − l.

Proposition 6 Let us assume that exist, for H ∈ {I, E, I} and all a, and let us assume the following conditions, for l ≥ 2, where Ψ1, and Ψ2, are independent of . Then the limit process of incidence of cases aged a at time n is a single-type Markovian process of order d with a Poissonian transition law: where I := ∑ I , where the {Y} are i.i.d. given F = {I1, I2, ...} with Y1|F1 ∼ Poisson(Ψ), Ψ = ∑ Ψ {Y} are independent.

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, instead of . Then the assumptions of this proposition 6 were satisfied under the following assumptions [28]:

The S and E individuals have the same time-homogeneous survival law {S};

There is no over-contamination during the incubation period or the clinical state;

The number of newborn animals at each birth per individual at time l, is independent of l, i, and of the health state h of i (but the health state of each newborn and his survival during the first time unit may depend on h);

The population is roughly stable: ;

The disease is rare at the initial time: ;

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 where 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 and 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), φ ∈ [0, 1] represents the efficiency at time n − l of the main current control regulations: φ = 1, for k ≤ 1988, φ = φ, for k ∈ (1989, 1996), and φ = 0, for k ≥ 1997, where φ represents the efficiency of the Meat and Bone Meal ban of July 1988, p 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(l) is the probability for an infected animal to have an intrinsic incubation time equal to l (conditioned on survival). So Ψ 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 Φ depends only on l. This is achieved by modeling the process on a period with a constant control. Then {I} may be written as a multitype BGW branching process and therefore results of Proposition 4 are valid for this process, replacing Φ by Ψ := Ψ. More precisely, let Ĩ =: (I1, I2, ..., I) := (I1, ..., I(1)) and F1 = {I1, I2, ..., I0}.

Let M͂ be defined by E(Ĩ1) =: Ĩ1M͂, where

Proposition 7 {Ĩ} is a multitype BGW branching process: where, for j = 1, , for j > 1, and , for l ≠ j − 1, and the are independent given F–1 and .

Moreover E(Ĩ) = Ĩ0ρ are the first eigenvalue and corresponding eigenvector with ξ1 = 1, of M͂, that is, M͂ = ρ , l = 1, …, d.

Let s := (s1, .., s), s+1 := 1, , f(s) := (f(1)(s), ..., f()(s)).

Proposition 8 The generating function of {Ĩ}, is equal to where f(s) := f(f1(s)), and f()(s) = s+1 exp(−Ψ(1 − s1)).

As in Proposition 4, the bifurcation parameter of the process is given by the first eigenvalue ρ of M͂.

Proposition 9 ρ is solution of , and ρ ≤ 1 is equivalent to R ≤ 1, where is the total mean number of new cases I who begin to be generated by a I during his first time unit.

Let us assume from now on that Ψ1 > 0,..., Ψ > 0.

Proposition 10

P(lim I = 0 ∪ lim I = ∞) = 1;

ρ > 1 is equivalent to R∞ > 1 (supercritical case) which itself implies the existence of a positive integrable random variable W such that , where P(lim I = ∞) = P(W > 0). Moreover let us assume that I− = c, l = 1, ..., d − 1, where c is independent of I0. Then .

ρ ≤ 1 is equivalent to R∞ ≤ 1 (subcritical and critical cases) which implies P(lim I = 0) = 1 (a.s. extinction).

Let us point out that the deterministic model derived from E(ĨĨ1) = Ĩ1 M͂ is defined by: X͂ := X͂1M͂, X͂0 := Ĩ0. Therefore X͂ = X͂0M͂ and the bifurcation parameter of {X͂} is ρ, but for ρ = 1, X͂ = X͂0 (persistence), while (extinction).

Let T be the (random) extinction time of {I}. Then the extinction probability is q := P(lim Ĩ = 0) = P(T < ∞).

Proposition 11 We have , where , is the extinction probability starting from Ĩ0 = (0, .., 0, 1, 0, ..., 0) (one individual of the h type), h = 1, ..., d.

As a consequence, in the particular case Ĩ0 = (I0, 0, ..., 0), then , where q1 is solution of the equation: 1 − q1 = 1 − exp(−R(1 − q1)). Thus, there exists a solution q1 ≠ 1 if R > 1, implying that q1 decreases as R increases. Otherwise, if R ≤ 1, then q1 = 1, implying q = 1, for all h = 1, ..., d.

Let G := P(T ≤ n) = P(Ĩ = 0) (distribution of the extinction time T.). Then according to [1], p.187, if R < 1, the distribution of the extinction time has an exponential form as n → ∞: lim ρ−(1 − G) = Q(0)Ĩ0. > 0 where Q(0) := lim ρ−.(1 − f(0)), M = ρ, = ρ, = 1, .1 = 1, and for h = 1, ..., d, ζ ∝ ρ−(1), . So for n large enough, P(T > n) ∝ρ, but the coefficient of proportionality is not easy to calculate in practice. So we investigate other ways to calculate G.

Proposition 12 For any value of R ≥ d, Let us notice that G decreases (and R increases) as soon as at least one Ψ increases.

In addition, using Proposition 8, we get the following result, allowing an exact iterative computation of {G}:

Proposition 13 For any value of R {G} of T satisfies: where s0 = (0, ..., 0) and f(s) = (f1(s), ..., f(s)) is the offspring generating function (see Proposition 8). In the particular single-type case d = 1, .

Let be the size of the tree generated by Ĩ0 (epidemic size). Let us recall that the Borel-Tanner distribution with parameter (λ, l) is a (f) GLPD (Generalized Lagrange Probability Distribution), where fλ (s) = e−(1) is the generating function of the Poisson distribution with parameter λ [7, 14] (see (2)).

Let g(s) := E(s) be the generating function of N.

Proposition 14 Let us assume that (subcritical case).

Let Ĩ0 := (I0, 0, ..., 0). Then , that is : and , .

Let Ĩ0 := (I0, I−1, , ..., I−(−1)). Then , where , the {N} are i.i.d. with and ⊕ means the mutual independence, that is 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 , we proved, for a large N0, that the proportion of infected individuals in the whole population, , behaves, as n → ∞, as the probability 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} 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 . 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 N0, 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 N0 is too small, the limit in N0 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 (R0), 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 R0. The extinction time roughly increases as R0 increases. So the greatest difference between the behavior of this population-dependent process and its deterministic counterpart is obtained when R0 > 1 with R0 ≃ 1. This is the generalization of the difference observed between the BGW process and its deterministic counterpart, when R0(= ρ)= 1.