Search tips
Search criteria 


Logo of ijerphMDPI Open Access JournalsMDPI Open Access JournalsThis articleThis JournalInstructions for authorsAdd your e-mail address to receive forthcoming issues of this journal
Int J Environ Res Public Health. 2010 March; 7(3): 1186–1204.
Published online 2010 March 19. doi:  10.3390/ijerph7031204
PMCID: PMC2872325

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.

Keywords: branching process, age-dependence, population-dependence, extinction time, epidemic size

1. 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 SIS or SIR, 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.

2. Single-type Bienaymé-Galton-Watson Branching Processes

Let In be the incidence of clinical (or diagnosed) cases at time n and Yn,i represent the numbers of secondary cases generated at time n by the previous case i. Then In is modelled from the past incidences by


where the variables {Yn,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 Fn−1 := {Ik}k≤n−1. We denote m, σ2 the mean and variance of Yn,1 given Fn−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): In=mIn1+In1ηn, where ηn:=[i=1In1(Yn,im)][In1]1 is asymptotically distributed according to N(0,σ2) on {limn In = ∞}, 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({limn In = } [union or logical sum] {limn In = 0}) = 1 and there exists an integrable random variable W such that limn Inmn = 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 L(m,σ2). This limit behavior comes from the martingale property of Inmn, that is E(Inmn|Fn−1) = In−1m−(n−1) which implies E(Inmn) = 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(Yn,1|Fn−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 Xn = mXn−1 with X0 = I0, derived from E(In|Fn−1) = In−1m, has the same bifurcation parameter, but when m = 1, the deterministic model persists with Xn = I0 while the stochastic one dies out a.s.. Moreover before dying out, the stochastic process In|In ≠ 0 presents a linear increasing behavior: limn P(In[0.5σ2n]−1x|In ≠ 0) = 1 − exp(−x).

In addition, the probability of extinction, q, is the solution of q = f(q) := E(qY1,1) and in the subcritical/critical cases m ≤ 1, the “epidemic size” N:=k=0Text.Ik (total number of cases generated until the extinction time Text.), follows a power series distribution when Yn,1 follows itself a power series distribution, that is, P(Yn,i = k) [proportional, variant] akλk [7, 11, 12]. In particular, when P(Yn,i = k) = exp(−m)[k!]−1mk (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.

3. 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, In=i=1In1Yn,i, where the {Yn,i}i are assumed i.i.d. according to a distribution L(Fn1)=:L(m(Fn1),σ2(Fn1)) given Fn−1 : {Ik}kn−1. The depletion effect of the S population at time n is replaced by the fact that m(Fn−1) (resp. P(Yn,i = 0|Fn−1)) is a decreasing (resp. increasing) function of the previous incidences of cases {Ik}kn−1.

A simple example is when L(Fn1)=L(l=1dnμlInl), μ [set membership] [0, 1], where m(l=1dnμlInl)(resp.P(Yn,1=0|Fn1)=p(l=1dnμlInl)) is a decreasing (resp. increasing) function of l=1dnμlInl. For example, E(Yn,1|Yn,1 ≠ 0) = α and P(Yn,1=0|Fn1)=1K[K+l=1dnμlInl]1, 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 dn = 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 dn = 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(Yn,i = 0|In−1 = N, Fn−1) ≥ f(N) > 0, for all N [set membership] N+ and any value of Fn−1 consistent with In−1 = N.

This assumption is checked as soon as P(Yn,i = 0|Fn−1) is a non-decreasing function of In−1, In−2,. . ., which is strictly increasing in In−1.

Proposition 1 Let us assume A1. Then P(limn In = 0 [union or logical sum] limn In = ∞) = 1.

Let A2: there exists m* < ∞ such that m(F) ≤ m*, for all F, and lim|F|→∞m(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 {In}.

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

Remark 1 When dn = d, for all n, {In} may be represented as a multitype branching process Ĩn = (In, In−1, ..., In−(d−1)) =: (In,1, In,2, ..., In,d), which is asymptotically BGW if L(Fn1) tends to a distribution L independent of the process as time tends to ∞ on {limn In = ∞}. But the asymptotic BGW process is not positively regular at the opposite of classical BGW branching processes, since (Fn−1) being defined by E(Ĩn|Fn−1) =: Ĩn−1(Fn−1), then


which implies that, for all nd − 1, n[1, j] = mn−(j−1), and n[i, j] = 0, i > 1, where m := lim|F|→∞ m(F) and := lim|F|→∞ (F).

Remark 2 If m((1, 0, ..., 0)) > 1 (supercritical assumption at the beginning of the disease spread), since m(Fn−1) < 1 under A2, for |Fn−1| and n sufficiently large, then we may observe oscillations until the a.s. extinction, the process increasing when |Ĩn−1|is small enough, and decreasing when it is too large. It is the case of the logistic model m(Fn1)=αK[K+l=1dnμlInl]1, when dn = d with d > 1 and m((1, 0, ..., 0)) = α(1 + μK−1)−1 > 1. Moreover since R0 = α > 1, then {Xn} tends to a stable limit cycle, as n → ∞ [9, 27].

4. 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 {Yn,i(j)}i are i.i.d. according to L(j)(Fn1)=:L(m(j)(Fn1),σ2(j)(Fn1)) given Fn1:{{Ik(j)}j=1,J}kn1. We assume that m(Fn−1) is non increasing in Ik(j), for each j, each kn − 1.

A simple example is when as previously L(j)(Fn1)=L(j)(l=1dnμljInl(j)), μ [set membership] [0, 1], and m(j)(l=1dnμljInl(j))(resp.P(Yn,1(j)=0|Fn1)=p(j)(l=1dnμljInl(j))) is a decreasing (resp. increasing) function of l=1dnμljInl(j). 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):=lim¯|F|m(j)(F)<1, then, for all j, P(limnIn(j)=0)=1, 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 In(j)=i=1In1(j)Yn,i(j), Yn,i(j)|Fn1Poisson(αK[K+ldμl(Inl(1)+ ...

Let the deterministic associate trajectory {Xn(j)} be defined by Xn(j):=m(j)(Fn1X)Xn1(j), n≥ 1, X0(j):=I0(j).

Proposition 3 The reproductive number R(j):=lim¯|F|m(j)(F) is the bifurcation parameter for {In(j)}, that is, if R(j)<1, then limnIn(j)=a.s.0, and the reproductive number R0(j):=lim_|FX|0m(j)(FX) is the bifurcation parameter of {Xn(j)}, that is if R0(j)<1, then limnXn(j)=0. Moreover if R0(j)>1, then lim¯nXn(j)>0.

A possible extension of the previous models consists in adding an immigration. Then:


where Jn 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 ({Yn,i}i, Jn) has a distribution independent of n and Fn−1, with either an immigration allowed only in the periods of extinction of the process {In} (regenerative branching process), or when the {δn} are i.i.d. independent of Fn−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 ({Yn,i}i, Jn, δn) is not time-dependent and when 1{In−1=0}Jnδn = 0, for all n, then (3) may be written as a model of Section 3..

5. Multitype Branching Processes with Age and Population-Dependent Offsprings

Let Nn=(Nn1,NnD) satisfying


where h, k represent types and l is the maturation time for getting the “offsprings” Ynl,n,i(h),k 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 {Ynl,n,i(h),k}i,l are mutually independent given Fn−1 = {Nn−1, ..., N0} and that the {Ynl,n,i(h),k}i are i.i.d. given Fn−1.

The models of the previous sections (Sections 2. and 4.) belong to this class when dn = 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” {Ynl,n,i(h),k} 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 {δ2,nl+1,n,iE,I|S}i are i.i.d. Bernoulli variables given Fn−1, with δ2,nl+1,n,iE,I|S equal to 1 if the individual i undergoes transitions SE, at time n − l + 1, and EI, at time n, which means that his incubation period is l, the {δ1,nl+1,n,i,jE,I|(h)}i,j are similar i.i.d. Bernoulli variables given Fn−1, relative to the newborn j of i, and the {Ynl+1,i(h)}i are i.i.d variables given Fn−1, Ynl+1,i(h) being the number of newborns of i at time n − l + 1. So Ynl,n,i(h),I 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 Nn := (Nn, Nn−1, ..., Nn−(d−1)). Then E(Nn|Fn−1) = Nn−1 M(Fn−1), where


I is the D × D identity matrix, and Mnl,n(Fn1)[h,k]:=E(Ynl,n,i(h),k|Fn1), for h = 1, ..., D, k = 1, ..., D. Then, assuming that M(F) is invertible, for any value F of Fn−1, Nn[ M(F0)... M(Fn−1)]−1 is a martingale with E(Nn [ M(F0)... M(Fn−1)]−1|F0) = N0, implying by the convergence theorem for martingale [15] that limn Nn [ M(F0)... M(Fn−1)]−1 Ut = WU, for any vector U, where WU is an integrable random variable. But, except in the simple cases M(Fn−1) = M (BGW process) or M(Fn1)=M(kNn1k) with limN M(N) = M (asymptotically BGW process) [30, 31, 32, 33], we cannot deduce from this martingale convergence any information on the asymptotic behavior of {Nn} 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.

5.1. The Population Size Remains Bounded

Let Nn:=kNnk 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 Yn,i at time n of the animal i is bounded whatever i, n, and we use the following control: if Nn−1 > NM, for some chosen threshold NM, then Yn,i = 0 (all the newborn animals are sold), and when Nn−1 < Nm, for some chosen threshold Nm < NM, then new animals are bought. Therefore assuming that P(Nn = N2|Nn−1 = N1) =: Q(N1, N2) is independent of n, {Nn} 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 limnNn[N0ρn]1=a.s.W determined by martingale theory in the frame of BGW processes.

5.2. 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: Nn/N0 =: Dn (empirical densities) or Nn/Nn =: Pn (empirical probabilities) or Nn[N0ρn]−1 =: Wn (empirical normalized densities), when assuming that limN0 D0 = D0 or limN0 P0 = 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]: In=i=1In1Yn,i, where the {Yn,i}i are i.i.d. given Fn−1 = {In−1, ..., I0} and depend on the previous cases incidences only through the condition: Yn,i = 0 when In−1 > cI0, that is, a massive vaccination is done as soon as the density In1I01 crosses the threshold c. This implies that the process dies out a.s.. The author studied Dn := In/I0 and proved that limI0 limn Dn/Dn ≠ 0 belongs to the set of stable fixed points of Dn, where Dn is the deterministic limit model Dn = Dn−1m(Dn−1), n ≥ 1, D0 = 1, with m(Dn−1) := E(Yn,1|Dn−1 = Dn−1). 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 limnlimN0P^n|Nn0=PlimN0limnP^n|Nn0 (and similarly for {Wn}), which allows to approximate limn Pn or limn Wn.

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 k=1DYnl,n,1(h),k only depends on l and is denoted k=1DYnl,n,i(h),k=Ynl,n,i. Let us denote Φl:=E(Ynl,n,1|Fn1),σl2:=Var(Ynl,n,1|Fn1). Then the total size of the alive population at time n, Nn=k=1DNnk, may be expressed as a multitype branching process Ñn := (Nn,1, ...Nn,d) := (Nn, Nn−1, ..., Nn−(d−1)), where Nn=l=1di=1NnlYnl,n,i and


with, for j = 1, Yn,i(l,1)=Ynl,n,i, for j > 1, Yn,i(j1,j)=1 and Yn,i(l,j)=0, for lj − 1, and the {Yn,i(l,j)}i are independent given Fn−1.

Moreover, assuming 0 < Φl < ∞, l = 1, …, d, 0<l=1dσl2<, and P(Yn−l,n,i ≥ 2|Fn−1) > 0 for some l, then {Nn}n satisfies the following properties:

  1. {Nn} is positively regular, nonsingular, and checks the xlogx property;
  2. P(limn Nn = 0 [union or logical sum] limn Nn = ∞) = 1
  3. Let ρ be the first eigenvalue of defined by E(Ñn|Ñn−1) =: Ñn−1, that is,
    Then E(Ñnξt) = Ñ0ρnξt, where ξt = ρξt. Moreover ρ ≤ 1 (subcritical and critical cases) implies that P(limn Nn = 0) = 1 (a.s. extinction), and ρ > 1 (supercritical case) implies the existence of an integrable random variable W such that limnWn=a.s.W , where Wn := Nn[N0ρn]−1, and P(limn Nn = ∞) = P(W > 0). Finally let us assume that Nl = ρ−lN0, l = 1, ..., d − 1. Then if ρ [set membership] Q, limN0Wn=a.s.limN0limnWn=a.s.1.
  4. ρ is solution of l=1dΦlρl=1, and ρ ≤ 1 [left and right double arrow ]R ≤ 1, where R=l=1dΦl (total mean number of offspring generated by an individual).

Let us notice that the process {Nn} 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, dn := dn, nd = [left floor]n/d[right floor] (integer part of n/d) and let us define:

A3: limnρn/2k=nd+1nl1,,lk[αl1ρl1/2αlkρlk/2]1{j=1kljn}<.

For d = 1, or d > 1 with [var phi]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 kn, α1kl1,,lk[αl1αlk]1{j=1kljn}1 and l1,,lk[αl1αlk]1{j=1kljn} is decreasing in k.

Proposition 5 Let us assume, as in Proposition 4, that the distribution of k=1DYnl,n,1(h),k depends only on l. Let us moreover assume ρ > 1, A3, Mn−l,n(Fn−1) = Mn−l,n(Pn−1) (density-dependence), limN0→∞P0 = P0, and Nl = clN0, l = 1, ..., d − 1, where cl is independent of N0. Then, for all η > 0, limN0 P(|PnPn| > η|Nn ≠ 0)=0 and


where Pn is the vector of probabilities defined by:


Proposition 5 allows to use the attractors of the dynamical model {Pn} 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: limN0limnNn[N0ρn]1=a.s.1 deduced from limnNn[N0ρn]1=a.s.W (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 {Pn} (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.

5.3. The Initial Population Size is Large and the Disease is Rare

Another approach for studying the asymptotic behavior of the process {Nn} 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): Nnk=l=1dh=1Di=1NnlhYnl,n,i(h),k, where the {Ynl,n,i(h),k}i,l are mutually independent given Fn−1 = {Nn−1, ..., N0} and the {Ynl,n,i(h),k}i are i.i.d. given Fn−1. 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 {δ2,nl+1,n,iE,I|h}i are i.i.d. Bernoulli variables given Fn−1, the {δ1,nl+1,n,i,jE,I|(h)}i,j are i.i.d. Bernoulli variables given Fn−1, the {Ynl+1,i(h)}i are i.i.d. given Fn−1, and the {δ2,nl+1,n,iE,I|h}i and the {δ1,nl+1,n,i,jE,I|(h)}i,j are mutually independent given Fn−1. Therefore the incidence of cases aged a at time n is given by:


where NY,k(h):=i=1Nk1hYk,i(h).

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 = Itot × a, where Itot represents the clinical state (new or not), then


where δnl,n,iI,Itot|I×al 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 {NnItot×a}a,n is easily calculated from the distribution of {NnI×a}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 NnI×a given Fn1*:={{NkI}kn1,{Nk}kn1,{NY,k}kn}, where NY,k=hNY,k(h). Of course, expressions similar to (5) and (6) can be written for NnE×a and NnEtot×a.

Let δn−l,h(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 limN0NnH×a exist, for H [set membership] {I, E, Itot, Etot} and all a, and let us assume the following conditions, for l ≥ 2,


where Ψ1,a,l|n−l and Ψ2,a,l|n−l are independent of Fnl*. Then the limit process Ia,n:=DlimN0NnI×a of incidence of cases aged a at time n is a single-type Markovian process of order d with a Poissonian transition law:



where In := ∑a Ia,n. Moreover we may write In=l=1di=1InlYnl,n,i, where the {Yn−l,n,i}i are i.i.d. given Fn−1 = {In−1, In−2, ...} with Yn−l,n,1|Fn−1 ~ Poissonl|n−l), Ψl|n−l = ∑a Ψa,l|n−l, and the {Yn−l,n,i}i,l 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, F˜nl*:={{NkI}knl+1,{Nk}knl,{NY,k}knl+1} instead of Fnl*. Then the assumptions of this proposition 6 were satisfied under the following assumptions [28]:

  1. The S and E individuals have the same time-homogeneous survival law {Sa}a;
  2. There is no over-contamination during the incubation period or the clinical state;
  3. The number of newborn animals Yl,i(h) 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);
  4. The population is roughly stable: limN0N1N01=a.s.1;
  5. The disease is rare at the initial time: limN0(N0E+N0I)2N01=a.s.0;
  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
    where NkI*tot 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 Ψl|nl=a=ld[(θal,Rc+θal,Rφnl)+1{a=l}pmat]P(a)Pinc.(l), where θlal,Rc and θlal,R(expλ1)1φnl 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), [var phi]n−l [set membership] [0, 1] represents the efficiency at time n − l of the main current control regulations: [var phi]k = 1, for k ≤ 1988, [var phi]k = [var phi], for k [set membership] (1989, 1996), and [var phi]k = 0, for k ≥ 1997, where [var phi] represents the efficiency of the Meat and Bone Meal ban of July 1988, pmat. 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 Pinc.(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 {In} 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 Ĩn =: (In,1, In,2, ..., In,d) := (In, In−1, ..., In−(d−1)) and Fn−1 = {In−1, In−2, ..., I0}.

Let be defined by E(Ĩn|Fn−1) =: Ĩn−1, where


Proposition 7 {Ĩn} is a multitype BGW branching process:


where, for j = 1, Yn,i(l,1):=Ynl,n,i, for j > 1, Yn,i(j1,j):=1 and Yn,i(l,j):=0, for lj − 1, and the {Yn,i(l,j)}i,l,j are independent given Fn–1 and Yn,i(l,1)|I˜n1Poisson(Ψl).

Moreover E(Ĩnξt) = Ĩ0ρnξt, where ρ and ξ are the first eigenvalue and corresponding eigenvector with ξ1 = 1, of , that is, ξt = ρξt, implying that ξl=ρl1l=ldρlΨl, l = 1, …, d.

Let s := (s1, .., sd), sd+1 := 1, f(l)(s):=E(s1Yn,1(l,1)sdYn,1(l,d)), f(s) := (f(1)(s), ..., f(d)(s)).

Proposition 8 The generating function of {Ĩn}, Fn(s):=E(s1In,1sdIn,d) is equal to


where fn(s) := f(fn−1(s)), and f(h)(s) = sh+1 exp(−Ψh(1 − s1)).

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

Proposition 9 ρ is solution of l=1dΨlρ1=1, and ρ ≤ 1 is equivalent to R ≤ 1, where R=l=1dΨl 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,..., Ψd > 0.

Proposition 10

  1. P(limn In = 0 [union or logical sum] limn In = ∞) = 1;
  2. ρ > 1 is equivalent to R > 1 (supercritical case) which itself implies the existence of a positive integrable random variable W such that limnIn[I0ρn]1=a.s.W, where P(limn In = ∞) = P(W > 0). Moreover let us assume that Il = clI0, l = 1, ..., d − 1, where cl is independent of I0. Then limI0In[I0ρn]1=a.s.limI0limnIn[I0ρn]1=a.s.1.
  3. ρ ≤ 1 is equivalent to R ≤ 1 (subcritical and critical cases) which implies P(limn In = 0) = 1 (a.s. extinction).

Let us point out that the deterministic model derived from E(Ĩn|Ĩn−1) = Ĩn−1 is defined by: n := n−1, 0 := Ĩ0. Therefore n = 0n and the bifurcation parameter of {n} is ρ, but for ρ = 1, nξt = 0ξt (persistence), while limnI˜n=a.s.0 (extinction).

Let Text. be the (random) extinction time of {In}. Then the extinction probability is q := PĨ0(limn Ĩn = 0) = PĨ0(Text. < ∞).

Proposition 11 We have q=q1I0qdI(d1), where qh=exp(l=hdΨl(1q1)), 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 q=q1I0, 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 qh = 1, for all h = 1, ..., d.

Let Gn := P(Text.n) = P(Ĩn = 0) (distribution of the extinction time Text.). Then according to [1], p.187, if R < 1, the distribution of the extinction time has an exponential form as n → ∞: limn ρn(1 − Gn) = Q(0)Ĩ0.ξt > 0 where Q(0) := limn ρnζ.(1fn(0))t, Mξt = ρξt, ζM = ρζ, ξ.ζt = 1, ξ.1t = 1, and for h = 1, ..., d, ζh [proportional, variant] ρ−(h−1), ξhρh1l=hdρ1Ψl. So for n large enough, P(Text. > n) [proportional, variant]ρn, but the coefficient of proportionality is not easy to calculate in practice. So we investigate other ways to calculate Gn.

Proposition 12 For any value of R, the extinction time distribution satisfies, for nd,


Let us notice that Gn 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 {Gn}:

Proposition 13 For any value of R, the distribution {Gn} of Text. satisfies:


where s0 = (0, ..., 0) and fn(s) = (fn,1(s), ..., fn,d(s)) is the offspring generating function (see Proposition 8). In the particular single-type case d = 1, Gn=[f(Gn11/I0)]I0=exp[ΨI0(1Gn11/I0)].

Let N:=l=(d1)Text.1Il 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λ, sl) GLPD (Generalized Lagrange Probability Distribution), where fλ (s) = eλ(1−s) is the generating function of the Poisson distribution with parameter λ [7, 14] (see (2)).

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

Proposition 14 Let us assume that l=1dΨl<1 (subcritical case).

  1. Let Ĩ0 := (I0, 0, ..., 0). Then g(s)=sI0fI0l=1dΨl(g(s)), that is NBorelTanner(l=1dΨl,I0):
    and E(N)=I0(1l=1dΨl)1, Var(N)=I0(l=1dΨl)(1l=1dΨl)3.
  2. Let Ĩ0 := (I0, I−1, , ..., I−(d−1)). Then N=Dh=1di=1I(h1)(j=1Y0,h,iNi,j+1), where Y0,h,iPoisson(l=hdΨl), the {Ni,j} are i.i.d. with Ni,jBorelTanner(l=1dΨl,1) and [plus sign in circle] means the mutual independence, that is

6. 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 limnlimN0P^n|Nn0=PlimN0limnP^n|Nn0, we proved, for a large N0, that the proportion of infected individuals in the whole population, P^nE+P^nI, behaves, as n → ∞, as the probability PnE+PnI 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 {Nn} 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 limnlimN0NnI=limN0limnNnI. 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 [similar, equals] 1. This is the generalization of the difference observed between the BGW process and its deterministic counterpart, when R0(= ρ)= 1.


1. Athreya KB, Ney PE. Branching Processes. Springer-Verlag; Berlin, Germany: 1972.
2. Becker NG. On parametric estimation for mortal branching processes. Biometrika. 1974;61:393–399.
3. Becker NG. Estimation for an epidemic model. Biometrics. 1976;32:769–777. [PubMed]
4. Becker NG. Estimation for discrete time branching processes with applications to epidemics. Biometrics. 1977;33:515–522. [PubMed]
5. Becker NG. A general chain binomial model for infectious diseases. Biometrics. 1981;37:251–258. [PubMed]
6. Becker NG. Analysis of Infectious Disease Data. Chapman and Hall/CRC Press; Boca Raton, FL, USA: 1989.
7. Devroye L. The branching process method in Lagrange random variate generation. Commun. Statistics-Simulat. Comput. 1992;21:1–14.
8. Diekmann O, Heesterbeek JAP. Mathematical Epidemiology of Infectious Diseases Model Building, Analysis and Interpretation Wiley Series in Mathematical and Computational Biology. John Wiley Ltd; Chichester, UK: 2000
9. Elaydi S. An Introduction to Difference Equations 3rd ed Undergraduate Texts in Mathematics. Springer; New York, NY, USA: 2005
10. Farrington CP, Grant AD. The distribution of time to extinction in subcritical branching processes: Applications to outbreaks of infectious disease. J. Appl. Probab. 1999;36:771–779.
11. Farrington CP, Kanaan MN, Gay NJ. Branching process models for surveillance of infectious diseases controlled by mass vaccination. Biostatistics. 2003;4:279–295. [PubMed]
12. Good IJ. The Lagrange distributions and branching processes. SIAM J. Appl. Math. 1975;28:270–275.
13. Haccou P, Jagers P, Vatutin VA. Branching Processes Variation, Growth, and Extinction of Populations. Cambridge Studies in Adaptative Dynamics; Cambridge, UK: 2005.
14. Haight FA, Breuer AB. The Borel-Tanner distribution. Biometrika. 1960;47:143–150.
15. Hall P, Heyde CC. Martingale Limit Theory and its Application. Probability and Mathematical Statistics; New York, NY, USA: 1980.
16. Harris TE. The Theory of Branching Processes Die Grundlehren der Mathematischen Wissenschaften, Bd 119 Springer-Verlag; Berlin, Germany: Prentice-Hall, Inc; Englewood Cliffs, N.J., USA: 1963
17. Heesterbeek JAP, Dietz K. The concept of Ro in epidemic theory. Stat. Neerl. 1996;50:89–110.
18. Hefferman JM, Smith RJ, Wahl LM. Perspectives on the basic reproductive ratio. J. R. Soc. Interface. 2005;2:281–293. [PMC free article] [PubMed]
19. Hethcote HW. The Mathematics of Infectious Diseases. SIAM Rev. 2000;42:599–653.
20. Heyde CC. On assessing the potential severity of an outbreak of a rare infectious disease: a Bayesian approach. Austral. J. Statist. 1979;21:282–292.
21. Jacob C, Viet AF. Epidemiological modeling in a branching population. Particular case of a general SIS model with two age classes. Math. Biosci. 2003;182:93–111. [PubMed]
22. Jacob C, Magal P. Influence of the routine slaughtering on the evolution of BSE. Example of the British and French slaughterings. Risk Anal. 2007;27:1151–1167. [PubMed]
23. Jacob C. Branching Processes: Their rôle in Epidemiology Technical report, UR341, MIA department, INRA; Jouy-en-Josas, France: 2009
24. Jagers P. Branching Processes with Biological Applications. Wiley Series in Probability and Mathematical Statistics, John Wiley; London, UK: 1975.
25. Jacob C, Lalam N, Yanev N. Statistical inference for processes depending on environments and application in regenerative processes. Pliska Stud. Math. Bulgar. 2005;17:109–136.
26. Jacob C, Viet AF. A new class of processes for formalizing and generalizing individual-based models: the semi-semi-Markov processes. Pliska Stud Math Bulgar. 2007:121–144.
27. Jacob C. Saturation Effects in Population Dynamics: Use Branching Processes or Dynamical Systems? In: Deutsch A, Bravo de la Parra R, de Boer R, Dieckmann O, Jagers P, Kisdi E, Kretzschmar M, Lansky P, Metz H, editors. Mathematical Modeling of Biological Systems. Birkhäuser; Boston, MA, USA: 2008. pp. 339–352.
28. Jacob C, Maillard-Teyssier L, Denis JB, Bidot C. A branching process approach for the propagation of the Bovine Spongiform Encephalopathy in Great-Britain Accepted in Branching processes and their Applications, Lecture Notes in Statistics, Springer-Verlag; Berlin, Germany: 2010. 227–242.242
29. Jagers P. A general stochastic model for population development. Skan. Aktuarietidskr. 1969;52:84–103.
30. Klebaner FC. Population-size-dependent branching process with linear rate of growth. J. Appl. Probab. 1983;20:242–250.
31. Klebaner FC. On population-size-dependent branching processes. Adv. Appl. Probab. 1984;16:30–55.
32. Klebaner FC. Linear growth in near-critical population-size-dependent multitype Galton-Watson processes. J. Appl. Probab. 1989;26:431–445.
33. Klebaner FC. Correction: “Linear growth in near-critical population-size-dependent multitype Galton-Watson processes” [J. Appl. Probab. 26 (1989), no. 3, 431–445] and “Asymptotic behavior of near-critical multitype branching processes” [ibid. 28 (1991), no. 3, 512–519] J. Appl. Probab. 1992;29:246.
34. Klebaner FC. Population-dependent branching processes with a threshold. Stochastic Process. Appl. 1993;46:115–127.
35. Klebaner FC, Nerman O. Autoregressive approximation in branching processes with a threshold. Stochastic Process. Appl. 1994;51:1–7.
36. Klebaner FC, Zeitouni O. The exit problem for a class of density-dependent branching systems. Ann. Appl. Probab. 1994;4:1188–1205.
37. Klebaner FC. Population and density dependent branching processes Classical and Modern Branching Processes IMA Vol Math Appl, 84, Springer; New York, NY, USA: 1997. 165–169.169
38. Klebaner FC, Lazar J, Zeitouni O. On the quasi-stationary distribution for some randomly perturbed transformations of an interval. Ann. Appl. Probab. 1998;8:300–315.
39. Yanev NM, Mitov KV. Regenerative Branching Processes Records and Branching processes Ahsanullah M, Yanev GP, editors. Nova Science Publishers, Inc; New York, NY, USA: 2008. Chapter 3; 37–62.62
40. Viet AF, Fourichon C, Seegers H, Jacob C, Guihenneuc-Jouyaux C. A model of the spread of the bovine viral-diarrhoea virus within a dairy herd. Prev. Vet. Med. 2004;63:211–236. [PubMed]

Articles from International Journal of Environmental Research and Public Health are provided here courtesy of Multidisciplinary Digital Publishing Institute (MDPI)