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”
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
Fn−1, 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
Fn−1, relative to the newborn
j of
i, and the
are i.i.d variables given
Fn−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.
5.1. 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
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(ℕ
n =
N2|ℕ
n−1 =
N1) =:
Q(
N1, N2) is independent of
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(
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.
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 =: D̂n (empirical densities) or Nn/Nn =: P̂n (empirical probabilities) or Nn[N0ρn]−1 =: Wn (empirical normalized densities), when assuming that limN0 D̂0 = D0 or limN0 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 {
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
crosses the threshold
c. This implies that the process dies out a.s.. The author studied
D̂n :=
In/
I0 and proved that lim
I0 lim
n D̂n/
D̂n ≠ 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|D̂n−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
(and similarly for {Wn}), which allows to approximate limn P̂n 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 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 Ñn := (
Nn,1, ...Nn,d) := (
Nn, Nn−1, ..., Nn−(d−1))
, where andwith, for j = 1,
,
for j > 1,
and,
for l ≠
j − 1,
and the are independent given Fn−1.
Moreover, assuming 0 < Φ
l < ∞,
l = 1, …,
d,
, and P(
Yn−l,n,i ≥ 2
|Fn−1) > 0
for some l, then {
Nn}
n satisfies the following properties:{Nn} is positively regular, nonsingular, and checks the xlogx property;
P(limn Nn = 0 ∪ limn Nn = ∞) = 1
Let ρ be the first eigenvalue of M͂ defined by E(
Ñn|
Ñn−1) =:
Ñn−1M͂,
that is,Then E(
Ñnξt) =
Ñ0ρnξt,
where M͂ξt =
ρξt. Moreover ρ ≤ 1
(subcritical and critical cases) implies that P(lim
n Nn = 0) = 1
(a.s. extinction), and ρ > 1
(supercritical case) implies the existence of an integrable random variable W such that , where Wn :=
Nn[
N0ρn]
−1, and P(lim
n Nn = ∞) =
P(
W > 0)
. Finally let us assume that N−l =
ρ−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 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 := d ∧ n, nd = ⌊n/d⌋ (integer part of n/d) and let us define:
A3:
.
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,
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,
Mn−l,n(
Fn−1) =
Mn−l,n(
P̂n−1)
(density-dependence), lim
N0→∞P̂0 =
P0, and N−l =
clN0, l = 1
, ..., d − 1
, where cl is independent of N0. Then, for all η > 0
, lim
N0 P(
|P̂n −
Pn| >
η|Nn ≠ 0)=0
andwhere 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:
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 {
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):
, where the
are mutually independent given
Fn−1 = {
Nn−1, ..., N0} and the
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
are i.i.d. Bernoulli variables given
Fn−1, the
are i.i.d. Bernoulli variables given
Fn−1, the
are i.i.d. given
Fn−1, and the
and the
are mutually independent given
Fn−1. 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 =
Itot ×
a, where
Itot 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 δ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 exist, for H ∈ {
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. 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 In := ∑
a Ia,n. Moreover we may write, 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 ∼
Poisson(Ψ
l|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,
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 {Sa}a;
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), φ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, 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
M͂ be defined by
E(
Ĩn|Fn−1) =:
Ĩn−1M͂, where
Proposition 7 {
Ĩn} is a multitype BGW branching process:where, for j = 1,
,
for j > 1,
and,
for l ≠
j − 1
, and the are independent given Fn–1 and.
Moreover E(Ĩnξt) = Ĩ0ρnξt, where ρ and ξ are the first eigenvalue and corresponding eigenvector with ξ1 = 1, of M͂, that is, M͂ξt = ρξt, implying that, l = 1, …, d.
Let s := (s1, .., sd), sd+1 := 1,
, f(s) := (f(1)(s), ..., f(d)(s)).
Proposition 8 The generating function of {
Ĩn},
is equal towhere 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 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,..., Ψd > 0.
Proposition 10P(limn In = 0 ∪ limn In = ∞) = 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(limn In = ∞) = P(W > 0). Moreover let us assume that I−l = clI0, l = 1, ..., d − 1, where cl is independent of I0. Then
.
ρ ≤ 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 M͂ is defined by: X͂n := X͂n−1M͂, X͂0 := Ĩ0. Therefore X͂n = X͂0M͂n and the bifurcation parameter of {X͂n} is ρ, but for ρ = 1, X͂nξt = X͂0ξt (persistence), while
(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, 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 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 → ∞: lim
n ρ−n(1
− Gn) =
Q(
0)
Ĩ0.
ξt > 0 where
Q(
0) := lim
n ρ−nζ.(
1 −
fn(
0))
t,
Mξt =
ρξt,
ζM =
ρζ,
ξ.ζt = 1,
ξ.
1t = 1, and for
h = 1
, ..., d,
ζh ∝
ρ−(h−1),
. So for
n large enough,
P(
Text. >
n) ∝
ρ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 n ≥
d,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,
.
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λ, 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 (subcritical case).
Let Ĩ0 := (
I0, 0, ..., 0)
. Then, that is:
and,
.
Let Ĩ0 := (
I0,
I−1, , ...,
I−(d−1))
. Then,
where,
the {
Ni,j}
are i.i.d. with and ⊕
means the mutual independence, that isand