Epidemic transmission on SEIR stochastic models with
Transcript Of Epidemic transmission on SEIR stochastic models with
Epidemic transmission on SEIR stochastic models with nonlinear incidence rate
M.J. Lopez-Herrero
[email protected] Faculty of Statistical Studies, Complutense University of Madrid
28040 Madrid
Our interest is to quantify the spread of an infective process with latency period and generic incidence rate that takes place in a …nite and homogeneous population.
Within a stochastic framework, two random variables are de…ned to describe the variations of the number of secondary cases produced by an index case inside of a closed population. Computational algorithms are presented in order to characterize both random variables. Finally, theoretical and algorithmic results are illustrated by several numerical examples.
Keywords: Nonlinear incidence rate; Stochastic SEIR epidemic model; Index case; Secondary cases.
1 Introduction
A very large number of epidemiological models involves a compartmental division of individuals. It is also commonly assumed that individuals make contacts at random showing no preferences. One of the fundamental compartmental models is the SIR proposed by Kermack & McKendrick [1], where individuals are classi…ed as susceptible to the infection, labeled S, infected by the pathogen, I, or recovered from the infection, R. There are many hypothesis behind this model, which refer to population’s size (a large one, a closed one, with or without natural births and deaths during the outbreaks), the lack of a latency period (individuals become infectious as soon as they become infected), lifetime immunity after recovery, and homogeneous mixing.
This paper regards the study of the spread of a communicable disease involving latency period and a general mode of transmission so, we generalize the classic SIR model in two ways. First we deal with an SEIR model, that is, we split the infected population group into two ones: latent or in exposed period individuals and infectious individuals able to spread the infection. Second, we consider a transmission of nonlinear type.
The mode of transmission determines the possible response of the disease upon study. Its mathematical description is based in a function called force
1
of infection or incidence rate. In most models transmission is assumed to occur via so called mass action [2], involving a bilinear function of S and I: This particular choice makes that theoretical models show an exponential growth in the early part of the epidemics ([1], [2]).
Bilinear incidence rate is a natural assumption in accordance with the homogeneous mixing hypothesis but it is dubious upon heterogeneous. In that sense, recent studies ([3], [4]), based on empirical infection diseases datasets, revealed that a slower than exponential initial growth is observed in populations with highly constrained contact structures that are su¤ering from infectious diseases transmitted via close contacts (e.g. Ebola, SARS, STDs, etc.). These works provide an evidence that considering nonlinear transmission rates seems a more realistic assumption than linearity, specially with regard to recover sub-exponential growth pro…les. Some authors ([5], [6]) generalize the bilinear rate to a power function, IpSq, or introduce a deviation from the bilinear choice by using a rate of the form IS(1+ I) (see [7] and the references therein). Korobeinikov & Maini [8] consider nonlinear transmission to incorporate the population structure and heterogeneous mixing into the model. More precisely, they consider an increasing function of S that is also increasing and concave respect to the number of infective individuals. Zhang et al. [9] and Lahrouz & Omari [10] consider a transmission rate of the form IS= (I), with (I) positive and non decreasing; which covers well-known incidence functions such as the proposed by Capasso & Serio [11] arising from saturation e¤ects ( (I) = 1 + I), or the one re‡ecting psychological e¤ects ( (I) = 1 + I2) proposed Xiao & Ruan [12]. More recently, Li et al. [13] choose a general incidence rate of the form SG(I), with G(I) not necessarily increasing, which includes the mass action, saturated incidence and incidence with media coverage ([14], [15]).
More examples of nonlinear infection rates and their e¤ect on the population can be found in the works by Hethcote & Levin [16], Alexander & Moghadas [17] and Turnes & Monteiro [18].
Measuring the transmission potential in a population is the role of R0, the basic reproduction number (see [19] or [20]). This is the most widely used measure of disease spread in epidemiology and population dynamics. R0 is a threshold parameter which tries to measure the initial spread of an epidemic process by estimating the average of secondary cases of infection arising from an index case in a virgin population, during its infectious period. Artalejo & Lopez-Herrrero include in [21] a summary of remarkable features and drawbacks of R0, which is supported by an important set of references. It is not the aim of the author to repeat same arguments here so, we kindly address the interested reader to [21] and the references therein.
2
In the stochastic framework it is frequently assumed a value for R0 inherited from the deterministic counterpart. In that sense, for our stochastic SEIR model with non-linear incidence rate, we call basic reproduction number to the quantity appearing at expression (3.7) in the paper by Korobeinikov & Maini [8], determined by van den Driessche & Watmough [22] using the next generation method, whenever the non-linear incidence rate satis…es conditions stated on [22]. We remark that, for a deterministic model without vital dynamics, the expression for R0, is the same in SIR and SEIR models; showing that the average number of individuals infected by an index case does not depend on the exposed population. And one wonders if this coincidence is not a consequence of estimating an average number in a large population.
When dealing with a small population, where the stochastic approach is preferable [23], the validity of the linearization hypothesis assumed in a deterministic framework is questionable. To neglect the impact of depletion of susceptible individuals due to the infective process produces an overestimation on the number of secondary cases which correspond to contacts between the index case and individuals previously infected.
In the paper we present two measures alternative to R0 focusing on a stochastic perspective. For an epidemic process described by a continuous time Markov chain (CTMC ) we quantify the steady-state behavior of the epidemic spread by means of two random variables. Namely, Re0: the exact reproduction number and Rp or the population transmission number. The aim of this paper is two fold:
To extent the existing methodology to a more involved epidemic model and to …nd stable computational schemes able to handle the high dimensionality of the state spaces arising when an additional component (i.e., latency period) is introduced in the Markovian model representing the epidemic evolution within the population.
To study the in‡uence of the latency parameter to the epidemic spread in an SEIR model, when transmission depends on a general nonlinear incidence function. In particular we investigate whether Re0 and Rp distributions depend or not on latency rate for a speci…c group of incidence functions including the usual bilinear mass action.
The organization of the rest of the paper is as follows. In Section 2, we introduce the mathematical description of the stochastic SEIR model as a continuous time Markov chain. In Sections 3 and 4, starting from a …xed initial state, we de…ne above variables and provide algorithmic recursive
3
schemes for the computation of the generating functions, mass probability functions and factorial moments of Re0 and Rp, respectively. Section 5 illustrates numerically the behavior of both random variables. Last, in Section 6, we discuss our results.
2 Model description
Consider a closed population of N individuals which is a¤ected by a contagious disease conferring immunity and having a latency period between being infected and becoming infective. At any time, individuals in the population are classi…ed into four homogeneous compartments according to his disease state. More speci…cally, at a given point time t we denote by S(t) the number of susceptible individuals, by E(t) the number of exposed or individuals in latent period, by I(t) the number of infectious individuals and by R(t) the number of recovered or immune individuals.
Latency periods are assumed to be independent and identical distributed with length described by an exponential law having a progression rate
. Each individual recovers independently of the rest of infectious individuals according to an exponential law with rate . Disease transmission is produced through direct contact between an infective individual and a susceptible; these contact occur at the time points of a time homogeneous Poisson process; infection rate will depend on the number of infective and susceptible through a …nite incidence function (s; i) si. After recovery, individuals play no role in the epidemic spread; consequently the epidemic stops as soon as the population contains no infected (i.e., exposed or infectious) individuals. In addition, all the involved Poisson processes are assumed to be independent.
As we are dealing with a closed population, there is no need of recording all compartment’s occupancy level. So, the number of recovered individuals is not reported and the stochastic SEIR model is a three dimensional continuous time Markov chain X = f(S(t); E(t); I(t)) : t > 0g that provides the number of susceptible, exposed and infectious individuals at any time t and whose state space is S = S1 S2 with S1 = f(s; e; i) : 0 s; e; i N; s + e + i N g and S2 = f(s; e; 0) : 0 s; e N; s + e = N g. Notice that S is a …nite set having jSj = N (N + 1)(N + 5)=6 states in. We assume that any epidemic outbreak starts with a single infectious individual in a totally susceptible population, hence the initial state of X is (S(0); E(0); I(0)) = (N 1; 0; 1): If we order the states in S according to natural order, that is, S = f(s; e; i) : 0 s N 1; 0 e N s 1; 0
4
i N s eg, we notice that a state (s; e; i) is accessible from another (s0; e0; i0) when (s; e; i) (s0; e0; i0) respect to the natural order above. More-
over, we observe that once the CTMC X departs from a state never returns
to it. In long terms, the process enters into the set of absorbing states
SA = f(s; 0; 0) : 0 s N 1g which correspond to the end of the infec-
tious disease within the closed population. Due to the reducible character
of X, the stationary distribution becomes degenerate in the sense that it
only assingns probability mass to states in SA, which guarantees the certain
extinction of the epidemic process.
Soujourn times at each state (s; e; i) are independent exponential random
variables with rates determined by means of the in…nitesimal generator Q;
whose entries are described as follows
8
>>>><
si; e;
if (s0; e0; i0) = (s 1; e + 1; i); if (s0; e0; i0) = (s; e 1; i + 1);
q(s;e;i)(s0;e0;i0) = >>>
i; q
;
if if
(s0; e0; i0) = (s; e; i 1); (s0; e0; i0) = (s; e; i);
>: sei
0;
otherwise.
where qsei = si + e + i, for (s; e; i) 2 S: Matrix Q can be expressed in a block lower triangular form by parti-
tioning the state space in level states: S = [Ns=01ls, where ls is the s-th level that contains states with exactly s susceptible individuals. Each level state
is partitioned, in turns, in sub-levels according to the number of exposed individuals, that is ls = [Ne=0s 1lse. Details can be found in [24]. This partition will be useful when searching for numerical solutions of systems of
equations trough recursive computationally stable schemes.
3 The exact reproduction number Re0
The distribution of the o¤spring distribution of secondary infections at a particular stage of an epidemic assuming Markov chain compartmental models was focused …rst by Ross in [25] and, in an independent way, by Artalejo and Lopez-Herrero in [21]. Both papers evaluate the probability mass function of secondary infections by iteratively solving a set of linear equations. Methodology in [25] is applied to phase type infectious periods, producing e¢ cient results for small populations but the size of the matrices involved are extremely large and computational costs prevent the extension of this methodology to more complicated models. On the other hand, methodology given in [21] involves …rst-step arguments and solves systems of linear
5
equations using a stable iterative procedure that avoids the use and storage of large matrices, which permits to handle populations of up to 10,000 individuals, when dealing with classical SIS and SIR models.
Our aim is to extend the underlying methodological ideas to characterize the distribution of the exact reproduction number, Re0, for the Markov chain SEIR compartmental model and observe the in‡uence of long latent periods on the epidemic spread.
In this Section we present stable algorithmic schemes for determining the distribution of Re0, via generating function, mass probability function and moments. Numerical experiments in Section 5, obtained by implementing their corresponding computational algorithms, will show that the distribution of the exact reproduction number depends upon the latency time even when the mass action transmission function is considered, in contradistinction to what happens in deterministic SIR and SEIR models without demographic parameters.
Let us begin by introducing the random variable Re0 associated to the infection expansion which is de…ned as follows:
Re0 is the exact number of secondary cases arising from a single infective individual during its entire infectious period.
We notice that this random variable is de…ned to correct the e¤ect of the linearization assumption because vain contacts (i.e., those contacts occurring between the selected infective and any contacted individual previously infected) are ignored. In particular, the basic reproduction number, R0, is de…ned as the expected value of Re0 conditioned to the initial state (N 1; 0; 1) for the process X, that is a single infective in a virgin population. As Re0 is a random variable we can get more information about it than an expected value under a …xed circumstance, so in what follows we proceed by studying the whole distribution (via generating and mass probability functions) and also its moments.
We start by introducing an appropriate notation for the generating and probability mass functions of Re0, and factorial moments conditioned to the current situation.
Xs
'sei(z) =
znP fRe0 = n j(S(0); E(0); I(0)) = (s; e; i) g; for jzj 1;
n=0
xnsei = P fRe0 = n j(S(0); E(0); I(0)) = (s; e; i) g; for 0 n s;
mlsei =
0; if l = 0; E[Re0(Re0
1) : : : (Re0
l + 1)]; if l > 0; ;
with (s; e; i) 2 ST = S SA.
6
In order to develop algorithmic schemes for determining 'sei(z); xnsei and mlsei we begin by marking the infective individual starting the infective process, next we split the transmission process by distinguishing whether a contact involves or not the initial marked infective individual. In consequence, transmission rate is decomposed as
si = si + esi; for (s; e; i) 2 ST and i 1;
where si = si=i and esi = (i 1) si=i. By appealing to a …rst step argument, we …nd that for a …xed state
(s; e; i) 2 ST the next transition from that state of the Markov chain is associated to one of the following events: (i) an e¤ective contact between the tagged infective and any susceptible individual, (ii) an e¤ective contact between a non- tagged infective and any susceptible individual, (iii) the end latency period for any exposed individual, (iv) the recovery of a nontagged infective and (v) the recovery of the marked infective individual. The following stable algorithm provides the generating function of the number of e¤ective contacts produced by the tagged infective during its infectious time, from the current state (s; e; i): Notation xy stands for Kronecker’s delta, de…ned as one for x = y and zero otherwise.
Algorithm 1 The conditional generating function 'sei(z) of the number of secondary cases produced by a typical infective given an a current state (s; e; i) 2 ST is computed with the help of the following scheme: Step 1a: Set s = 0, Step 1b: Set e = 0, Step 1c: Set i = 1 Step 1d: Compute
'sei(z) = (1 s0)
si
z's 1;e+1;i(z)
(1)
si + e + i
+ (1 + (1
e
)(1 i1)
si
's 1;e+1;i(z)
s0
si + e + i
e
e0) si + e + i 's;e 1;i+1(z)
+ (1
(i 1) i1) si + e + i 's;e;i 1(z) + si + e + i ;
Step 2a: Set i = i + 1. If i N s e, then go to Step 1d. Step 2b: Set e = e + 1. If e N s 1, then go to Step 1c. Step 2c: Set s = s + 1. If s N 1, then go to Step 1b.
7
By di¤erentiating l 1 times the equation (1) with regard to z and
setting z = 1, we obtain a new equation that involves conditional moments of Re0, that is, mlsei, for (s; e; i) 2 ST are determined from
mlsei = (1 + (1 + (1
s0) si + sei + i lmls 11;e+1;i
s0)
si
mls 1;e+1;i + (1
si + e + i
i1)
(i 1) ml ; i 1:
s;e;i 1
si + e + i
(2)
e0)
e
ml
s;e 1;i+1
si + e + i
We notice that for l = 0 moments satisfy m0sei = 1, whenever (s; e; i) 2 ST and i 1: So, there is a recursive stable algorithm, similar to Algorithm 1, that provides conditional moments up to any order l 1 based on moments of one order less. In more details, we start from order zero moments and substitute expression (1) in Algorithm 1 by equation (2). In particular, we can get numerical values for Re0 = E[Re0 j(S(0); E(0); I(0)) = (N 1; 0; 1) ] which provides the expected number of secondary cases produced by the single infective individual inserted in a completely susceptible population that is aimed to be registered by the basic reproduction number R0:
Algorithm 1 is also the starting point for getting the mass function of Re0, just by using a numerical inversion method (see for instance [26]) applied to generating functions provided by (1). Instead of that, we present a direct method for computing the conditional probability distribution of Re0 given the current state for the CTMC X.
A new appeal to the …rst step analysis methodology, based on the …rst transition from a state (s; e; i) 2 ST gives
x00ei = 1; 0 e N 1; 1 i N e;
(3)
xnsei = (1 s0)(1 n0) si + sei + i xns 11;e+1;i
e
+(1 s0)(1 i1)(1 ns)
si
xns 1;e+1;i
si + e + i
+(1 e0)
e
xns;e 1;i+1 + (1
i1)
(i 1) xn
s;e;i 1
si + e + i
si + e + i
+ n0 si + e + i ; (4)
which can be recursively solved directly after the procedure described in the next algorithm.
8
Algorithm 2 The set of probabilities xnsei describing the mass function of Re0 conditioned to the current situation (s; e; i) 2 ST , are determined as follows: Step 1: Set n = 0. Step 2a: Set s = 0. Step 2b: Set e = 0. Step 2c: Set i = 1. Step 2d: Compute xnsei via equations (3) or (4). Step 3a: Set i = i + 1. If i N s e, go to Step 2d. Step 3b: Set e = e + 1. If e N s 1, go to Step 2c. Step 3c: Set s = s + 1. If s N 1, go to Step 2b. Step 4: Set n = n + 1. If n N 1 go to Step 2a.
4 The population transmission number Rp
This section focuses on the random variable Rp, or the population transmission number, which is de…ned as the exact number of secondary cases produced by all currently infective individuals prior to the …rst recovery. This random variable was …rst studied in [21] and, for SIS and SIR epidemics models with classical rates, its distribution characteristics are determined through simple explicit formulas.
Assuming that we are able to identify the infection when …rst removal occurs, maybe because it is a death case, then we notice that if we interpret a recovery as a removal, the random variable Rp is a measure of the epidemic spread before the …rst detection for the outbreak of the epidemic process.
As we did in the previous section, the objective is to characterize the distribution of the random population transmission number in terms of its generating and mass probability function and also providing its moments. First we introduce some notation.
Let us condition on the current situation, described in terms of the CTMC X, for (s; e; i) 2 ST and jzj 1 function
Xs sei(z) = znP fRp = n j(S(0); E(0); I(0)) = (s; e; i) g
n=0
is the generating function of Rp. Given an integer l 0, the factorial moment of order l conditioned to the current situation is
Mslei = l0+(1 l0)E[Rp(Rp 1) (Rp l+1) j(S(0); E(0); I(0)) = (s; e; i) ]:
9
Finally, we denote by ysnei the probability of having n secondary cases prior to …rst recovery, starting from the situation (s; e; i) 2 ST , that is
ysnei = P fRp = n j(S(0); E(0); I(0)) = (s; e; i) g:
As in the precedent section, it is possible to derive generating, mass probability functions and factorial moments in terms of an algorithmic scheme; where recursive equations are obtained using a …rst step argument. More speci…cally, conditioning on the …rst possible event starting from the state (s; e; i) 2 ST we achieve the following set of equations involving the generating functions at jzj 1 :
0ei(z) = 1; 0 e N 1; 1 i N e;
(5)
sei(z) = (1 s0)(1 i0)
si
z s 1;e+1;i(z)
(6)
si + e + i
e
i
+(1 e0) si + e + i s;e 1;i+1(z) + (1 i0) si + e + i :
We can compute numerical values of sei(z) in the natural order, that is for 0 s N 1; 0 e N s 1 and e0 i N s e.
Taking l 1 derivatives on (5) and (6) regarding z, and evaluating at
z = 1, we get a new set of linear equations involving conditional factorial
moments. That is,
M0lei = 0; 0 e N 1; 1 i N e;
(7)
Mslei = (1 s0)(1 i0) si + sei + i Msl 11;e+1;i (8)
+(1 s0)(1 i0)
si
Msl 1;e+1;i
si + e + i
+(1 e0)
e
Ml
:
s;e 1;i+1
si + e + i
We observe that system of equations (7) - (8) involve conditional moments of one order less, which is an advantage for minimizing computational storage.
A new …rst step argument drives us to the following set of equations governing probabilities, ysnei for 0 n s and (s; e; i) 2 ST :
y00ei = 1; 0 e N 1; 1 i N e;
ysnei = (1 s0)(1 i0)(1 n0) si + sei + i ysn 11;e+1;i
+(1 e0)
e
ysn;e 1;i+1 + (1 i0) n0
i:
si + e + i
si + e + i
10
M.J. Lopez-Herrero
[email protected] Faculty of Statistical Studies, Complutense University of Madrid
28040 Madrid
Our interest is to quantify the spread of an infective process with latency period and generic incidence rate that takes place in a …nite and homogeneous population.
Within a stochastic framework, two random variables are de…ned to describe the variations of the number of secondary cases produced by an index case inside of a closed population. Computational algorithms are presented in order to characterize both random variables. Finally, theoretical and algorithmic results are illustrated by several numerical examples.
Keywords: Nonlinear incidence rate; Stochastic SEIR epidemic model; Index case; Secondary cases.
1 Introduction
A very large number of epidemiological models involves a compartmental division of individuals. It is also commonly assumed that individuals make contacts at random showing no preferences. One of the fundamental compartmental models is the SIR proposed by Kermack & McKendrick [1], where individuals are classi…ed as susceptible to the infection, labeled S, infected by the pathogen, I, or recovered from the infection, R. There are many hypothesis behind this model, which refer to population’s size (a large one, a closed one, with or without natural births and deaths during the outbreaks), the lack of a latency period (individuals become infectious as soon as they become infected), lifetime immunity after recovery, and homogeneous mixing.
This paper regards the study of the spread of a communicable disease involving latency period and a general mode of transmission so, we generalize the classic SIR model in two ways. First we deal with an SEIR model, that is, we split the infected population group into two ones: latent or in exposed period individuals and infectious individuals able to spread the infection. Second, we consider a transmission of nonlinear type.
The mode of transmission determines the possible response of the disease upon study. Its mathematical description is based in a function called force
1
of infection or incidence rate. In most models transmission is assumed to occur via so called mass action [2], involving a bilinear function of S and I: This particular choice makes that theoretical models show an exponential growth in the early part of the epidemics ([1], [2]).
Bilinear incidence rate is a natural assumption in accordance with the homogeneous mixing hypothesis but it is dubious upon heterogeneous. In that sense, recent studies ([3], [4]), based on empirical infection diseases datasets, revealed that a slower than exponential initial growth is observed in populations with highly constrained contact structures that are su¤ering from infectious diseases transmitted via close contacts (e.g. Ebola, SARS, STDs, etc.). These works provide an evidence that considering nonlinear transmission rates seems a more realistic assumption than linearity, specially with regard to recover sub-exponential growth pro…les. Some authors ([5], [6]) generalize the bilinear rate to a power function, IpSq, or introduce a deviation from the bilinear choice by using a rate of the form IS(1+ I) (see [7] and the references therein). Korobeinikov & Maini [8] consider nonlinear transmission to incorporate the population structure and heterogeneous mixing into the model. More precisely, they consider an increasing function of S that is also increasing and concave respect to the number of infective individuals. Zhang et al. [9] and Lahrouz & Omari [10] consider a transmission rate of the form IS= (I), with (I) positive and non decreasing; which covers well-known incidence functions such as the proposed by Capasso & Serio [11] arising from saturation e¤ects ( (I) = 1 + I), or the one re‡ecting psychological e¤ects ( (I) = 1 + I2) proposed Xiao & Ruan [12]. More recently, Li et al. [13] choose a general incidence rate of the form SG(I), with G(I) not necessarily increasing, which includes the mass action, saturated incidence and incidence with media coverage ([14], [15]).
More examples of nonlinear infection rates and their e¤ect on the population can be found in the works by Hethcote & Levin [16], Alexander & Moghadas [17] and Turnes & Monteiro [18].
Measuring the transmission potential in a population is the role of R0, the basic reproduction number (see [19] or [20]). This is the most widely used measure of disease spread in epidemiology and population dynamics. R0 is a threshold parameter which tries to measure the initial spread of an epidemic process by estimating the average of secondary cases of infection arising from an index case in a virgin population, during its infectious period. Artalejo & Lopez-Herrrero include in [21] a summary of remarkable features and drawbacks of R0, which is supported by an important set of references. It is not the aim of the author to repeat same arguments here so, we kindly address the interested reader to [21] and the references therein.
2
In the stochastic framework it is frequently assumed a value for R0 inherited from the deterministic counterpart. In that sense, for our stochastic SEIR model with non-linear incidence rate, we call basic reproduction number to the quantity appearing at expression (3.7) in the paper by Korobeinikov & Maini [8], determined by van den Driessche & Watmough [22] using the next generation method, whenever the non-linear incidence rate satis…es conditions stated on [22]. We remark that, for a deterministic model without vital dynamics, the expression for R0, is the same in SIR and SEIR models; showing that the average number of individuals infected by an index case does not depend on the exposed population. And one wonders if this coincidence is not a consequence of estimating an average number in a large population.
When dealing with a small population, where the stochastic approach is preferable [23], the validity of the linearization hypothesis assumed in a deterministic framework is questionable. To neglect the impact of depletion of susceptible individuals due to the infective process produces an overestimation on the number of secondary cases which correspond to contacts between the index case and individuals previously infected.
In the paper we present two measures alternative to R0 focusing on a stochastic perspective. For an epidemic process described by a continuous time Markov chain (CTMC ) we quantify the steady-state behavior of the epidemic spread by means of two random variables. Namely, Re0: the exact reproduction number and Rp or the population transmission number. The aim of this paper is two fold:
To extent the existing methodology to a more involved epidemic model and to …nd stable computational schemes able to handle the high dimensionality of the state spaces arising when an additional component (i.e., latency period) is introduced in the Markovian model representing the epidemic evolution within the population.
To study the in‡uence of the latency parameter to the epidemic spread in an SEIR model, when transmission depends on a general nonlinear incidence function. In particular we investigate whether Re0 and Rp distributions depend or not on latency rate for a speci…c group of incidence functions including the usual bilinear mass action.
The organization of the rest of the paper is as follows. In Section 2, we introduce the mathematical description of the stochastic SEIR model as a continuous time Markov chain. In Sections 3 and 4, starting from a …xed initial state, we de…ne above variables and provide algorithmic recursive
3
schemes for the computation of the generating functions, mass probability functions and factorial moments of Re0 and Rp, respectively. Section 5 illustrates numerically the behavior of both random variables. Last, in Section 6, we discuss our results.
2 Model description
Consider a closed population of N individuals which is a¤ected by a contagious disease conferring immunity and having a latency period between being infected and becoming infective. At any time, individuals in the population are classi…ed into four homogeneous compartments according to his disease state. More speci…cally, at a given point time t we denote by S(t) the number of susceptible individuals, by E(t) the number of exposed or individuals in latent period, by I(t) the number of infectious individuals and by R(t) the number of recovered or immune individuals.
Latency periods are assumed to be independent and identical distributed with length described by an exponential law having a progression rate
. Each individual recovers independently of the rest of infectious individuals according to an exponential law with rate . Disease transmission is produced through direct contact between an infective individual and a susceptible; these contact occur at the time points of a time homogeneous Poisson process; infection rate will depend on the number of infective and susceptible through a …nite incidence function (s; i) si. After recovery, individuals play no role in the epidemic spread; consequently the epidemic stops as soon as the population contains no infected (i.e., exposed or infectious) individuals. In addition, all the involved Poisson processes are assumed to be independent.
As we are dealing with a closed population, there is no need of recording all compartment’s occupancy level. So, the number of recovered individuals is not reported and the stochastic SEIR model is a three dimensional continuous time Markov chain X = f(S(t); E(t); I(t)) : t > 0g that provides the number of susceptible, exposed and infectious individuals at any time t and whose state space is S = S1 S2 with S1 = f(s; e; i) : 0 s; e; i N; s + e + i N g and S2 = f(s; e; 0) : 0 s; e N; s + e = N g. Notice that S is a …nite set having jSj = N (N + 1)(N + 5)=6 states in. We assume that any epidemic outbreak starts with a single infectious individual in a totally susceptible population, hence the initial state of X is (S(0); E(0); I(0)) = (N 1; 0; 1): If we order the states in S according to natural order, that is, S = f(s; e; i) : 0 s N 1; 0 e N s 1; 0
4
i N s eg, we notice that a state (s; e; i) is accessible from another (s0; e0; i0) when (s; e; i) (s0; e0; i0) respect to the natural order above. More-
over, we observe that once the CTMC X departs from a state never returns
to it. In long terms, the process enters into the set of absorbing states
SA = f(s; 0; 0) : 0 s N 1g which correspond to the end of the infec-
tious disease within the closed population. Due to the reducible character
of X, the stationary distribution becomes degenerate in the sense that it
only assingns probability mass to states in SA, which guarantees the certain
extinction of the epidemic process.
Soujourn times at each state (s; e; i) are independent exponential random
variables with rates determined by means of the in…nitesimal generator Q;
whose entries are described as follows
8
>>>><
si; e;
if (s0; e0; i0) = (s 1; e + 1; i); if (s0; e0; i0) = (s; e 1; i + 1);
q(s;e;i)(s0;e0;i0) = >>>
i; q
;
if if
(s0; e0; i0) = (s; e; i 1); (s0; e0; i0) = (s; e; i);
>: sei
0;
otherwise.
where qsei = si + e + i, for (s; e; i) 2 S: Matrix Q can be expressed in a block lower triangular form by parti-
tioning the state space in level states: S = [Ns=01ls, where ls is the s-th level that contains states with exactly s susceptible individuals. Each level state
is partitioned, in turns, in sub-levels according to the number of exposed individuals, that is ls = [Ne=0s 1lse. Details can be found in [24]. This partition will be useful when searching for numerical solutions of systems of
equations trough recursive computationally stable schemes.
3 The exact reproduction number Re0
The distribution of the o¤spring distribution of secondary infections at a particular stage of an epidemic assuming Markov chain compartmental models was focused …rst by Ross in [25] and, in an independent way, by Artalejo and Lopez-Herrero in [21]. Both papers evaluate the probability mass function of secondary infections by iteratively solving a set of linear equations. Methodology in [25] is applied to phase type infectious periods, producing e¢ cient results for small populations but the size of the matrices involved are extremely large and computational costs prevent the extension of this methodology to more complicated models. On the other hand, methodology given in [21] involves …rst-step arguments and solves systems of linear
5
equations using a stable iterative procedure that avoids the use and storage of large matrices, which permits to handle populations of up to 10,000 individuals, when dealing with classical SIS and SIR models.
Our aim is to extend the underlying methodological ideas to characterize the distribution of the exact reproduction number, Re0, for the Markov chain SEIR compartmental model and observe the in‡uence of long latent periods on the epidemic spread.
In this Section we present stable algorithmic schemes for determining the distribution of Re0, via generating function, mass probability function and moments. Numerical experiments in Section 5, obtained by implementing their corresponding computational algorithms, will show that the distribution of the exact reproduction number depends upon the latency time even when the mass action transmission function is considered, in contradistinction to what happens in deterministic SIR and SEIR models without demographic parameters.
Let us begin by introducing the random variable Re0 associated to the infection expansion which is de…ned as follows:
Re0 is the exact number of secondary cases arising from a single infective individual during its entire infectious period.
We notice that this random variable is de…ned to correct the e¤ect of the linearization assumption because vain contacts (i.e., those contacts occurring between the selected infective and any contacted individual previously infected) are ignored. In particular, the basic reproduction number, R0, is de…ned as the expected value of Re0 conditioned to the initial state (N 1; 0; 1) for the process X, that is a single infective in a virgin population. As Re0 is a random variable we can get more information about it than an expected value under a …xed circumstance, so in what follows we proceed by studying the whole distribution (via generating and mass probability functions) and also its moments.
We start by introducing an appropriate notation for the generating and probability mass functions of Re0, and factorial moments conditioned to the current situation.
Xs
'sei(z) =
znP fRe0 = n j(S(0); E(0); I(0)) = (s; e; i) g; for jzj 1;
n=0
xnsei = P fRe0 = n j(S(0); E(0); I(0)) = (s; e; i) g; for 0 n s;
mlsei =
0; if l = 0; E[Re0(Re0
1) : : : (Re0
l + 1)]; if l > 0; ;
with (s; e; i) 2 ST = S SA.
6
In order to develop algorithmic schemes for determining 'sei(z); xnsei and mlsei we begin by marking the infective individual starting the infective process, next we split the transmission process by distinguishing whether a contact involves or not the initial marked infective individual. In consequence, transmission rate is decomposed as
si = si + esi; for (s; e; i) 2 ST and i 1;
where si = si=i and esi = (i 1) si=i. By appealing to a …rst step argument, we …nd that for a …xed state
(s; e; i) 2 ST the next transition from that state of the Markov chain is associated to one of the following events: (i) an e¤ective contact between the tagged infective and any susceptible individual, (ii) an e¤ective contact between a non- tagged infective and any susceptible individual, (iii) the end latency period for any exposed individual, (iv) the recovery of a nontagged infective and (v) the recovery of the marked infective individual. The following stable algorithm provides the generating function of the number of e¤ective contacts produced by the tagged infective during its infectious time, from the current state (s; e; i): Notation xy stands for Kronecker’s delta, de…ned as one for x = y and zero otherwise.
Algorithm 1 The conditional generating function 'sei(z) of the number of secondary cases produced by a typical infective given an a current state (s; e; i) 2 ST is computed with the help of the following scheme: Step 1a: Set s = 0, Step 1b: Set e = 0, Step 1c: Set i = 1 Step 1d: Compute
'sei(z) = (1 s0)
si
z's 1;e+1;i(z)
(1)
si + e + i
+ (1 + (1
e
)(1 i1)
si
's 1;e+1;i(z)
s0
si + e + i
e
e0) si + e + i 's;e 1;i+1(z)
+ (1
(i 1) i1) si + e + i 's;e;i 1(z) + si + e + i ;
Step 2a: Set i = i + 1. If i N s e, then go to Step 1d. Step 2b: Set e = e + 1. If e N s 1, then go to Step 1c. Step 2c: Set s = s + 1. If s N 1, then go to Step 1b.
7
By di¤erentiating l 1 times the equation (1) with regard to z and
setting z = 1, we obtain a new equation that involves conditional moments of Re0, that is, mlsei, for (s; e; i) 2 ST are determined from
mlsei = (1 + (1 + (1
s0) si + sei + i lmls 11;e+1;i
s0)
si
mls 1;e+1;i + (1
si + e + i
i1)
(i 1) ml ; i 1:
s;e;i 1
si + e + i
(2)
e0)
e
ml
s;e 1;i+1
si + e + i
We notice that for l = 0 moments satisfy m0sei = 1, whenever (s; e; i) 2 ST and i 1: So, there is a recursive stable algorithm, similar to Algorithm 1, that provides conditional moments up to any order l 1 based on moments of one order less. In more details, we start from order zero moments and substitute expression (1) in Algorithm 1 by equation (2). In particular, we can get numerical values for Re0 = E[Re0 j(S(0); E(0); I(0)) = (N 1; 0; 1) ] which provides the expected number of secondary cases produced by the single infective individual inserted in a completely susceptible population that is aimed to be registered by the basic reproduction number R0:
Algorithm 1 is also the starting point for getting the mass function of Re0, just by using a numerical inversion method (see for instance [26]) applied to generating functions provided by (1). Instead of that, we present a direct method for computing the conditional probability distribution of Re0 given the current state for the CTMC X.
A new appeal to the …rst step analysis methodology, based on the …rst transition from a state (s; e; i) 2 ST gives
x00ei = 1; 0 e N 1; 1 i N e;
(3)
xnsei = (1 s0)(1 n0) si + sei + i xns 11;e+1;i
e
+(1 s0)(1 i1)(1 ns)
si
xns 1;e+1;i
si + e + i
+(1 e0)
e
xns;e 1;i+1 + (1
i1)
(i 1) xn
s;e;i 1
si + e + i
si + e + i
+ n0 si + e + i ; (4)
which can be recursively solved directly after the procedure described in the next algorithm.
8
Algorithm 2 The set of probabilities xnsei describing the mass function of Re0 conditioned to the current situation (s; e; i) 2 ST , are determined as follows: Step 1: Set n = 0. Step 2a: Set s = 0. Step 2b: Set e = 0. Step 2c: Set i = 1. Step 2d: Compute xnsei via equations (3) or (4). Step 3a: Set i = i + 1. If i N s e, go to Step 2d. Step 3b: Set e = e + 1. If e N s 1, go to Step 2c. Step 3c: Set s = s + 1. If s N 1, go to Step 2b. Step 4: Set n = n + 1. If n N 1 go to Step 2a.
4 The population transmission number Rp
This section focuses on the random variable Rp, or the population transmission number, which is de…ned as the exact number of secondary cases produced by all currently infective individuals prior to the …rst recovery. This random variable was …rst studied in [21] and, for SIS and SIR epidemics models with classical rates, its distribution characteristics are determined through simple explicit formulas.
Assuming that we are able to identify the infection when …rst removal occurs, maybe because it is a death case, then we notice that if we interpret a recovery as a removal, the random variable Rp is a measure of the epidemic spread before the …rst detection for the outbreak of the epidemic process.
As we did in the previous section, the objective is to characterize the distribution of the random population transmission number in terms of its generating and mass probability function and also providing its moments. First we introduce some notation.
Let us condition on the current situation, described in terms of the CTMC X, for (s; e; i) 2 ST and jzj 1 function
Xs sei(z) = znP fRp = n j(S(0); E(0); I(0)) = (s; e; i) g
n=0
is the generating function of Rp. Given an integer l 0, the factorial moment of order l conditioned to the current situation is
Mslei = l0+(1 l0)E[Rp(Rp 1) (Rp l+1) j(S(0); E(0); I(0)) = (s; e; i) ]:
9
Finally, we denote by ysnei the probability of having n secondary cases prior to …rst recovery, starting from the situation (s; e; i) 2 ST , that is
ysnei = P fRp = n j(S(0); E(0); I(0)) = (s; e; i) g:
As in the precedent section, it is possible to derive generating, mass probability functions and factorial moments in terms of an algorithmic scheme; where recursive equations are obtained using a …rst step argument. More speci…cally, conditioning on the …rst possible event starting from the state (s; e; i) 2 ST we achieve the following set of equations involving the generating functions at jzj 1 :
0ei(z) = 1; 0 e N 1; 1 i N e;
(5)
sei(z) = (1 s0)(1 i0)
si
z s 1;e+1;i(z)
(6)
si + e + i
e
i
+(1 e0) si + e + i s;e 1;i+1(z) + (1 i0) si + e + i :
We can compute numerical values of sei(z) in the natural order, that is for 0 s N 1; 0 e N s 1 and e0 i N s e.
Taking l 1 derivatives on (5) and (6) regarding z, and evaluating at
z = 1, we get a new set of linear equations involving conditional factorial
moments. That is,
M0lei = 0; 0 e N 1; 1 i N e;
(7)
Mslei = (1 s0)(1 i0) si + sei + i Msl 11;e+1;i (8)
+(1 s0)(1 i0)
si
Msl 1;e+1;i
si + e + i
+(1 e0)
e
Ml
:
s;e 1;i+1
si + e + i
We observe that system of equations (7) - (8) involve conditional moments of one order less, which is an advantage for minimizing computational storage.
A new …rst step argument drives us to the following set of equations governing probabilities, ysnei for 0 n s and (s; e; i) 2 ST :
y00ei = 1; 0 e N 1; 1 i N e;
ysnei = (1 s0)(1 i0)(1 n0) si + sei + i ysn 11;e+1;i
+(1 e0)
e
ysn;e 1;i+1 + (1 i0) n0
i:
si + e + i
si + e + i
10