SEIRS epidemics in growing populations

TL;DR

This paper models SEIRS epidemics in growing populations using stochastic and deterministic methods, analyzing various outbreak scenarios and their impact on population dynamics and disease prevalence.

Contribution

It introduces a combined stochastic-deterministic framework for SEIRS epidemics in super-critical growing populations, exploring different outbreak outcomes.

Findings

Disease can die out quickly or lead to endemic equilibrium.

Epidemic growth can outpace population growth or cause population decline.

Various parameter regimes determine epidemic and population trajectories.

Abstract

An SEIRS epidemic with disease fatalities is introduced in a growing population (modelled as a super-critical linear birth and death process). The study of the initial phase of the epidemic is stochastic, while the analysis of the major outbreaks is deterministic. Depending on the values of the parameters, the following scenarios are possible. i) The disease dies out quickly, only infecting few; ii) the epidemic takes off, the \textit{number} of infected individuals grows exponentially, but the \textit{fraction} of infected individuals remains negligible; iii) the epidemic takes off, the \textit{number} of infected grows initially quicker than the population, the disease fatalities diminish the growth rate of the population, but it remains super critical, and the \emph{fraction} of infected go to an endemic equilibrium; iv) the epidemic takes off, the \textit{number} of infected individuals grows initially quicker than the population, the diseases fatalities turn the exponential growth of the population to an exponential decay.

Tables

Table 1. Table 1: The uniform Markovian dynamic SEIRS epidemic: type of events, their state change l 𝑙 l (The old state z = ( u , v , x , y ) 𝑧 𝑢 𝑣 𝑥 𝑦 z=(u,v,x,y) is hence changed to z + l 𝑧 𝑙 z+l ) and their rates.

Event	State change $l$	Rate
Birth	$(1,0,0,0)$	$\lambda N$
Death of susceptible	$(-1,0,0,0)$	$\mu S$
Death of exposed	$(0,-1,0,0)$	$\mu E$
Death of infective	$(0,0,-1,0)$	$\left(\mu +\sigma \right)I$
Death of recovered	$(0,0,0,-1)$	$\mu R$
Infection	$(-1,1,0,0)$	$\kappa SI/N$
End of latency	$(0,-1,1,0)$	$\nu E$
Recovery	$(0,0,-1,1)$	$\delta I$
Loss of immunity	$(1,0,0,-1)$	$\rho R$

Table 2. Table 2: Thresholds for sub-models of the SEIRS model

Model	$R_0$	$\alpha$	Parameters change
$SEIRS$	$\frac{\kappa \nu }{\left(\mu +\nu \right)\left(\mu +\delta +\sigma \right)}$	$-\left(\mu +\frac{\nu +\delta +\sigma }{2}\right)+\sqrt{\frac{{\left(\nu -\delta -\sigma \right)}^2}{4}+\kappa \nu }$
$SEIR$	$\frac{\kappa \nu }{\left(\mu +\nu \right)\left(\mu +\delta +\sigma \right)}$	$-\left(\mu +\frac{\nu +\delta +\sigma }{2}\right)+\sqrt{\frac{{\left(\nu -\delta -\sigma \right)}^2}{4}+\kappa \nu }$	$\rho =0$
$SEIS$	$\frac{\kappa \nu }{\left(\mu +\nu \right)\left(\mu +\delta +\sigma \right)}$	$-\left(\mu +\frac{\nu +\delta +\sigma }{2}\right)+\sqrt{\frac{{\left(\nu -\delta -\sigma \right)}^2}{4}+\kappa \nu }$	$\rho \to \infty$
$SEI$	$\frac{\kappa \nu }{\left(\mu +\nu \right)\left(\mu +\sigma \right)}$	$-\left(\mu +\frac{\nu +\sigma }{2}\right)+\sqrt{\frac{{\left(\nu -\sigma \right)}^2}{4}+\kappa \nu }$	$\delta =0$, $\rho =0$
$SIRS$	$\frac{\kappa }{\mu +\delta +\sigma }$	$\kappa -\left(\mu +\delta +\sigma \right)$	$\nu \to \infty$
$SIR$	$\frac{\kappa }{\mu +\delta +\sigma }$	$\kappa -\left(\mu +\delta +\sigma \right)$	$\nu \to \infty$, $\rho =0$
$SIS$	$\frac{\kappa }{\mu +\delta +\sigma }$	$\kappa -\left(\mu +\delta +\sigma \right)$	$\nu \to \infty$, $\rho \to \infty$
$SI$	$\frac{\kappa }{\mu +\sigma }$	$\kappa -\left(\mu +\sigma \right)$	$\nu \to \infty$, $\delta =0$, $\rho =0$

Table 3. Table 3: Estimation of the probability of minor outbreak π n subscript 𝜋 𝑛 \pi_{n} of the epidemic starting with one latent individual and n − 1 𝑛 1 n-1 susceptible individuals. The theoretical result is π 𝜋 \pi and the simulated result is π ^ n subscript ^ 𝜋 𝑛 \hat{\pi}_{n} from 1000 1000 1000 simulations.

$n$	$1000$	$2000$	$1000$	$2000$	$1000$	$2000$	1000	$2000$	$1000$	$2000$
$\kappa$	$10$	$10$	$20$	$20$	$30$	$30$	$50$	$50$	$100$	$100$
$R_0$	$0.654$	$0.654$	$1.307$	$1.307$	$1.961$	$1.961$	$3.268$	$3.268$	$6.536$	$6.536$
$\pi$	1.000	$1.000$	$0.770$	$0.770$	$0.520$	$0.520$	$0.320$	$0.320$	$0.170$	$0.170$
${\hat{\pi }}_n$	$1.000$	$1.000$	$0.788$	$0.782$	$0.513$	$0.520$	$0.303$	$0.329$	$0.165$	$0.172$

Equations

$${\int }_0^{\infty }{e}^{-\alpha t}c(t)dt=1,$$

$$\alpha =-\left(\mu +\frac{\nu +\delta +\sigma }{2}\right)+\sqrt{\frac{(\nu -\delta -\sigma {)}^2}{4}+\kappa \nu }.$$

$$c(t)=\kappa e^{-\mu t}{\int }_0^t\nu e^{-\nu s}{e}^{-(\delta +\sigma )(t-s)}ds=\kappa \nu e^{-(\mu +\delta +\sigma )t}{\int }_0^t{e}^{-(\nu -(\delta +\sigma ))s}ds.$$

$$\displaystyle\displaystyle c(t)=\left\{\begin{array}[]{ll}\displaystyle\frac{\kappa\nu}{\nu-(\delta+\sigma)}\left(e^{-(\mu+\delta+\sigma)t}-e^{-(\nu+\mu)t}\right)\textrm{, if }\nu\neq\delta+\sigma,\\ \\ \displaystyle\kappa\nu te^{-(\mu+\delta+\sigma)t}\textrm{, if }\nu=\delta+\sigma.\end{array}\right.$$

$$\displaystyle\displaystyle\alpha=\left\{\begin{array}[]{ll}\displaystyle-\left(\mu+\frac{\nu+\delta+\sigma}{2}\right)+\sqrt{\frac{(\nu-\delta-\sigma)^{2}}{4}+\kappa\nu}\textrm{, if }\nu\neq\delta+\sigma,\\ \\ \displaystyle\sqrt{\kappa\nu}-(\mu+\delta+\sigma)\textrm{, if }\nu=\delta+\sigma.\end{array}\right.$$

$$\alpha =-\left(\mu +\frac{\nu +\delta +\sigma }{2}\right)+\sqrt{\frac{(\nu -\delta -\sigma {)}^2}{4}+\kappa \nu }.$$

$$R_0=\frac{\nu \kappa }{(\mu +\nu )(\delta +\mu +\sigma )}.$$

$$P(Y=0)=\frac{\mu }{\mu +\nu }+\frac{\nu }{\nu +\mu }\times \frac{\delta +\mu +\sigma }{\kappa +\delta +\mu +\sigma }.$$

$$P(Y=k)=\frac{\nu }{\mu +\nu }\times {\left(\frac{\kappa }{\kappa +\delta +\mu +\sigma }\right)}^k\times \frac{\delta +\mu +\sigma }{\kappa +\delta +\mu +\sigma }.$$

$$E(Y)=\frac{\nu }{\mu +\nu }\times \frac{1-p}{p}=\frac{\nu }{\mu +\nu }\times \frac{\kappa }{\mu +\delta +\sigma }.$$

$$R_0=\frac{\nu }{\mu +\nu }\times \frac{\kappa }{\mu +\delta +\sigma }.$$

$$\alpha >0$$

$$⟺\frac{(\nu -\delta -\sigma {)}^2}{4}+\kappa \nu >{\left(\mu +\frac{\nu +\delta +\sigma }{2}\right)}^2$$

$$⟺4\kappa \nu >(2\mu +\nu +\delta +\sigma {)}^2-(\nu -\delta -\sigma {)}^2$$

$$⟺4\kappa \nu >(2\mu +2\nu )(2\mu +2\delta +2\sigma )$$

$$⟺\frac{\kappa \nu }{(\mu +\nu )(\mu +\delta +\sigma )}>1$$

$$⟺R_0>1.$$

$$R_0=\kappa \times \frac{\nu }{\mu +\nu }\times \frac{1}{\mu +\delta +\sigma }.$$

$$g(z)$$

$$=\frac{\mu }{\mu +\nu }+\frac{\nu }{\mu +\nu }\frac{\delta +\mu +\sigma }{\kappa +\delta +\mu +\sigma }+\sum _{k=1}^{\infty }\frac{\nu }{\mu +\nu }\frac{\delta +\mu +\sigma }{\kappa +\delta +\mu +\sigma }{\left(\frac{\kappa }{\kappa +\delta +\mu +\sigma }\right)}^k{z}^k$$

$$=a+\frac{(1-a)b}{1-(1-b)z},\text{ with }a=\frac{\mu }{\mu +\nu }\text{ and }b=\frac{\delta +\mu +\sigma }{\kappa +\delta +\mu +\sigma }.$$

$$z=a+\frac{(1-a)b}{1-(1-b)z}.$$

$$z_0=1\text{ and }{z}_1=a+\frac{b}{1-b}=\frac{\mu }{\mu +\nu }+\frac{\delta +\mu +\sigma }{\kappa }=\frac{\mu }{\mu +\nu }+\frac{\nu }{\nu +\mu }\frac{1}{R_0}.$$

$$\displaystyle\displaystyle\pi=\left\{\begin{array}[]{ll}\displaystyle 1&\textrm{ if }R_{0}\leq 1,\\ \\ \displaystyle\frac{\mu}{\mu+\nu}+\frac{\nu}{\mu+\nu}\frac{1}{R_{0}}&\textrm{ if }R_{0}>1.\end{array}\right.$$

$$\displaystyle\displaystyle\pi_{(m,k)}=\left\{\begin{array}[]{ll}\displaystyle 1&\textrm{ if }R_{0}\leq 1,\\ \\ \displaystyle\left(\frac{\mu}{\mu+\nu}+\frac{\nu}{\mu+\nu}\frac{1}{R_{0}}\right)^{m}\left(\frac{1}{R_{0}}\right)^{k}&\textrm{ if }R_{0}>1.\end{array}\right.$$

$$\displaystyle\displaystyle\pi_{(0,1)}=\left\{\begin{array}[]{ll}\displaystyle 1&\textrm{ if }R_{0}\leq 1,\\ \\ \displaystyle\frac{1}{R_{0}}&\textrm{ if }R_{0}>1.\end{array}\right.$$

$$\begin{array}{rl}& \frac{dS}{dt}=\lambda N+\rho R-\kappa S\frac{I}{N}-\mu S,\\ & \frac{dE}{dt}=\kappa S\frac{I}{N}-(\nu +\mu )E,\\ & \frac{dI}{dt}=\nu E-(\delta +\mu +\sigma )I,\\ & \frac{dR}{dt}=\delta I-(\rho +\mu )R,\\ & N=S+E+I+R,\\ & S(0)>0,E(0)>0,I(0)>0,R(0)\ge 0.\end{array}$$

$$\begin{array}{rl}(\rho +\mu )\nu \kappa \sigma \lambda +\rho \kappa \nu \delta (\lambda -\mu )& \\ -(\rho +\mu )(\nu +\mu )(\delta +\mu +\sigma )[\kappa (\lambda -\mu )+\mu \sigma ]& =0,\end{array}$$

$$\displaystyle\left\{\begin{array}[]{l}S(0)=(\nu\kappa)^{-1}(\nu+\mu)(\delta+\mu+\sigma)N(0),\\ \\ E(0)=(\nu\sigma)^{-1}(\delta+\mu+\sigma)(\lambda-\mu)N(0),\\ \\ I(0)=\sigma^{-1}(\lambda-\mu)N(0),\\ \\ R(0)=(\sigma(\rho+\mu))^{-1}\delta(\lambda-\mu)N(0),\\ \\ \textrm{ with }N(0)>0.\end{array}\right.$$

$$\begin{array}{rl}& \frac{ds}{dt}=\lambda -\lambda s+\rho r+(\sigma -\kappa )si,\\ & \frac{de}{dt}=\kappa si-(\lambda +\nu )e+\sigma ei,\\ & \frac{di}{dt}=\nu e-(\lambda +\delta +\sigma )i+\sigma i^2,\\ & \frac{dr}{dt}=-(\lambda +\rho )r+\delta i+\sigma ri,\\ & s+e+i+r=1.\end{array}$$

$$\begin{array}{rl}& \frac{ds}{dt}=\lambda +\rho -(\lambda +\rho )s-\rho e-\rho i+(\sigma -\kappa )si,\\ & \frac{de}{dt}=\kappa si-(\lambda +\nu )e+\sigma ei,\\ & \frac{di}{dt}=\nu e-(\lambda +\delta +\sigma )i+\sigma i^2,\end{array}$$

$$D=\{(s,e,i);s\ge 0,e\ge 0,i\ge 0,s+e+i\le 1\}.$$

$$\begin{array}{rl}& \frac{de}{dt}=\kappa si-(\lambda +\nu )e+\sigma ei,\\ & \frac{di}{dt}=\nu e-(\lambda +\delta +\sigma )i+\sigma i^2.\end{array}$$

$$\begin{array}{rl}& \frac{de}{dt}=\kappa i-(\lambda +\nu )e,\\ & \frac{di}{dt}=\nu e-(\lambda +\delta +\sigma )i.\end{array}$$

Full text

1 Introduction

Infectious diseases remain a threat for developing countries as well as for developed countries. Many mathematicians focus their efforts to understand the dynamics of infectious diseases, in order to find the conditions to eradicate them. In mathematical modelling of infectious disease epidemics, the population in which the disease is spreading is partitioned in several compartments according to the status of the individuals, related to the disease. Every epidemic model has at least, the compartment $I$ of the infectious individuals who are infected and able to transmit the disease to others through contact, and the compartment $S$ of the susceptible individuals (those who are not infected but may be infected if they contact an infectious individual). Two other compartments often used are the compartment $E$ of the exposed or latent individuals who are already infected but not yet able to transmit the disease to others, and the compartment $R$ of the recovered or removed individuals (those who are healed from the disease with a permanent or non-permanent immunity). In a $SEIR$ epidemic, a susceptible individual infected through a contact with an infectious, becomes infected and latent; at the end of the latent period he/she becomes infectious and at the end of the infectious period he/she recovers with a life-long immunity. An $SEIRS$ epidemic is almost the same as the preceding, the only difference is that a recovered individual has a non-permanent immunity (He/she can lose his immunity, becoming susceptible again). Diphtheria, influenza and pneumonia are examples of diseases with latent period and non-permanent immunity [11].

In [5], Britton and Trapman studied stochastic SIR and SEIR models in a growing population. They derived the basic reproduction number and the Malthusian parameter of the epidemic, stated results for the initial phase and showed that the stochastic proportions process converges to a deterministic process.

In [11], Greenhalgh studied an SEIRS deterministic model with vaccination and found that under some conditions, the solution has Hopf bifurcations.

The aim of this paper is to study the dynamic of a stochastic SEIRS epidemic model with disease induced mortality, in an exponentially growing population. As in [5] , we assume that without the disease, the population has a birth rate $\lambda$, and a natural death rate $\mu$, such that $\lambda>\mu$. That is, initially the population process is a super-critical linear birth and death process. An SEIRS epidemic is introduced by infecting one individual. With the disease, the population is divided in four compartments according to the status of the individuals, related to the disease. The compartment S of the susceptible individuals, the compartment E of the latent or exposed individuals , the compartment I of the infectious individuals , and the compartment R of the removed individuals (those who are healed of the disease with a non-permanent immunity). The process is initiated by setting $(S(0),E(0),I(0),R(0))=(n-1,1,0,0)$. The transfer diagram of the model is given by Figure 1.

We derive the Malthusian parameter $\alpha$, the basic reproduction number $R_{0}$ and the probability of minor outbreak $\pi$ of the epidemic. If $R_{0}\leq 1$, then the disease cannot invade the population, that remains a super critical process. If $R_{0}>1$, then the epidemic has a positive probability $1-\pi$ of taking off, with the remaining probability $\pi$, it dies out. If the epidemic takes off, another threshold parameter $R_{1}$ determines the behavior of the proportion of infected individuals. If $R_{1}\leq 1$, then the fraction of infected stays small; while it persists when $R_{1}>1$. If $R_{1}>1$, or equivalently $\alpha>\lambda-\mu$, then the number of the infected grows initially quicker than the population, the disease fatalities diminish the growth rate of the population. In this case the asymptotic behavior of the population rely on a third threshold quantity $R_{2}$. If $R_{2}>1$, then the population goes on growing, while it becomes a sub-critical process when $R_{2}\leq 1$. In the latter case, when the number of individuals become low, the population should vanish with the disease, or regrows after the extinction of the epidemic.

We start by defining the stochastic model in Section $2$. Then, in Section $3$, we study the initial phase of the epidemic, thereafter we consider the deterministic model in Section $4$. Afterward, we give some illustrations by simulating different scenarios of epidemics in Section $5$. In Section $6$, we conclude the paper and discuss some perspectives.

2 The model

2.1 The initial dynamic of the population

Initially (before the introduction of the disease), the population model is a linear birth and death (B-D) process with individual birth rate $\lambda$ and individual death rate $\mu$. We assume that $\lambda>\mu$, that is the process is super-critical. $N(t)$ denotes the number of individuals in the population at time $t$.

2.2 The Markovian SEIRS epidemic model

Now, we define a uniformly mixing Markovian epidemic model on the population described above, implying that individuals give birth at rate $\lambda$ and die from other causes at rate $\mu$, irrespective of disease status. Initially (at $t=0$), the population consists of $n$ susceptible individuals. At this time, an SEIRS infectious disease is introduced by infecting one individual. The disease spreading is modelled as follow. An infected individual remains latent (infected but not yet infectious) for exponential time with rate $\nu$. After this period, he/she becomes infectious (unless he/she dies before). An infectious individual remains infectious for an exponential time with rate $\delta$ (unless he/she dies before). The disease induces an additional death rate $\sigma$ for the infectious individuals. After the period of infectiousness, an infectious individual recovers with a temporary immunity. A recovered person remains immune for an exponential time with rate $\rho$ unless he/she dies before. After the period of immunity he/she becomes susceptible again. During the infectious period, the infective has infectious contacts randomly in time according to a homogeneous Poisson process with rate $\kappa$, each time with a uniformly selected random individual. If the contacted person is susceptible, he/she becomes infected and latent (not yet infectious), otherwise the contact has no effect. There is no vertical transmission, that is all the new born individuals are susceptible.

Let $Z(t)=(S(t),E(t),I(t),R(t))$ respectively denote the number of susceptible, latent, infectious and immune individuals at time $t$. Therefore, the population size at time $t$ is $N(t)=S(t)+E(t)+I(t)+R(t)$ (Wherever $n$ is important an $n$-index is added). The population is initiated at $Z_{n}(0)=(S_{n}(0),E_{n}(0),I_{n}(0),R_{n}(0))=(n-1,1,0,0)$. However, later we also derive the probability of minor outbreak for an epidemic starting with $m$ latent individuals and $j$ infectious individuals with a very small infected fraction ($m+j\ll n$).

In short, the population model has two parameters, the birth rate $\lambda$ and the death rate $\mu$; and the disease model has five parameters, the transmission rate $\kappa$, the end of latency rate $\nu$, the recovery rate $\delta$, the disease death rate $\sigma$ and the immunity waning rate $\rho$. The possible events and their rates, when currently in state $Z(t)=(S(t),E(t),I(t),R(t))=(u,v,x,y)=z$, are given in Table 1.

The SI, SIS, SIR, SIRS, SEI, SEIS, SEIR models are special cases of the SEIRS model defined above. If $\rho=0$, then an infected individual cannot go back to the susceptible state and we get an SEIR model. If $\rho\longrightarrow\infty$, then the recovered state vanishes and we get an SEIS model. If $\nu\longrightarrow\infty$, then the latent state vanishes giving an SIRS model. If $\delta=0$, then we have lifelong infectivity and hence an SEI model. If $\nu\longrightarrow\infty$ and $\rho\longrightarrow\infty$, then the latent state and the recovered state vanish and we have an SIS model. If $\nu\longrightarrow\infty$ and $\delta=\rho=0$, then the latent state vanishes and we have lifelong infectivity, giving hence an SI model. Therefore from the results of the SEIRS model, one can deduce those of the others.

In this section we have presented the stochastic SEIRS model studied in this paper.

3 Results for the initial phase

3.1 The dynamics of the population size N(t)

Without the disease, the population process is a linear super-critical birth and death process with individual birth rate $\lambda$ and individual death rate $\mu$. But, with the introduction of the disease that induces an extra death rate $\sigma$ for the infectious individuals, we have two possible events. Birth with rate $\lambda N(t)$ and death with rate $\mu N(t)+\sigma I(t)=[\mu+\sigma I(t)/N(t)]N(t)$, where $N(t)$ is the current total number of individuals in the population and $I(t)$ the current number of infectious individuals. Then, the population process is no longer a linear birth and death process (unless $\sigma=0$). However, at the beginning of the epidemic when the fraction of the infectious individuals is very small ($I(t)/N(t)\approx 0$), the population will behave almost as a linear birth and death process with birth rate $\lambda$ and death rate $\mu$. Thus for the initial phase of the epidemic, we assume that the population size is a linear super-critical birth and death process with birth rate $\lambda$ and death rate $\mu$. On the other hand, if the epidemic takes off after the initial phase, whenever the fraction of the infectious individuals remains below $(\lambda-\mu)/\sigma$, the population process will be a super-critical process, having a positive probability to grow to $\infty$. But, if the infectious fraction grows beyond $(\lambda-\mu)/\sigma$, then the population process becomes a sub-critical process.

Remark 3.1**.**

If $(\lambda-\mu)/\sigma\geq 1$ i.e. $\lambda\geq\mu+\sigma$, then N(t) is always a super-critical process although the epidemic.

3.2 Approximation of the initial phase of the epidemic

Now, we consider the beginning of the epidemic, when the fraction of the infected individuals is still small. $E(t)$ increases with rate $\kappa I(t)S(t)/N(t)$ due to infection and decreases by rate $(\nu+\mu)E(t)$ due to death or end of latency. $I(t)$ increases with rate $\nu E(t)$ due to end of latency, and decreases with rate $(\delta+\mu+\sigma)I(t)$ due to recovery or death. As we start with $n-1$ susceptible and one infected, assuming that $n$ is large, at the beginning of the epidemic $S(t)/N(t)$ is very close to $1$. Then, the number of exposed individuals increases almost with rate $\kappa I(t)$. Let $L_{n}(t)=E_{n}(t)+I_{n}(t)$ be the number of infected (latent or infectious) individuals at time $t$. $L_{n}(t)$ can be approximated by a branching process, $L_{\infty}(t)=E_{\infty}(t)+I_{\infty}(t)$ with two stages: the childhood (or latent) stage $E_{\infty}$ and the adult (or infectious) stage $I_{\infty}$ [9, p. 54]. Initiated at $(E_{\infty}(0),I_{\infty}(0))=(1,0)$, where $E_{\infty}(t)$ increases with rate $\kappa I_{\infty}(t)$ and decreases with rate $(\mu+\nu)E_{\infty}(t)$, and $I_{\infty}(t)$ increases with rate $\nu E_{\infty}(t)$ (end of childhood) and decreases with rate $(\delta+\mu+\sigma)I_{\infty}(t)$ (end of adult stage).

Theorem 3.2**.**

Let $L_{n}(t)$ be the epidemic process and $L_{\infty}(t)$ be the branching process defined above. Then, $L_{n}(t)$ converges weakly to $L_{\infty}(t)$ ($L_{n}\stackrel{{\scriptstyle w}}{{\Longrightarrow}}L_{\infty}$), as $n\longrightarrow\infty$, on any finite interval $[0,t_{1}]$.

Proof.

Like initial phase of other epidemics [9, p. 54], when $n$ tends to infinity, the transition probabilities of the epidemic converge to that of the branching process. ∎

The results below are for the branching process $L_{\infty}$, but since $L_{n}\stackrel{{\scriptstyle w}}{{\Longrightarrow}}L_{\infty}$, they apply to the epidemic as $n\longrightarrow\infty$.

3.3 Thresholds

In this subsection, we derive the Malthusian parameter $\alpha$ and the basic reproduction number $R_{0}$ of the limiting branching process $L_{\infty}$.

The Malthusian parameter $\alpha$ is defined as the exponential growth/decay rate the epidemic branching process has. It is the solution of

[TABLE]

where $c(t)$ is the expected rate at which an individual gives birth (has infectious contacts) $t$ time units after it was infected [16, page 10].

Theorem 3.3**.**

The Malthusian parameter of the epidemic is given by

[TABLE]

Proof.

In the SEIRS model, the contact rate is [math] during the latent period and $\kappa$ during the infectious period. By conditioning on when the latent period ends, it follows that

[TABLE]

Thus,

[TABLE]

Inserting this into Equation (3.1), one gets

[TABLE]

Then,

[TABLE]

∎

The basic reproduction number $R_{0}$, is the expected number of secondary cases per primary case in a virgin population [9, page 4].

Theorem 3.4**.**

The basic reproduction number of the epidemic is

[TABLE]

Proof.

Let $Y$ be the number of infectious contacts that an individual has during the infectious period. Then,

[TABLE]

The first term is the probability that the individual dies during the latent period, and the second term the probability that the individual does not die during the latent period but leaves the infectious compartment by death or recovery without an infectious contact. And for all positive integer $k$,

[TABLE]

Therefore, $Z$ has a zero modified geometric distribution [13, page 16], with parameter $p=(\delta+\mu+\sigma)/(\kappa+\delta+\mu+\sigma)$ . Then the expected value of $Y$ is

[TABLE]

Thus,

[TABLE]

∎

Remark 3.5**.**

The first factor is the probability that the individual does not die during the latent period, and the second is the expected number of infectious contacts while being infectious. An alternative way to derive $R_{0}$ is by the relation $R_{0}=\int_{0}^{\infty}c(t)dt$ [13, page 69]. This formula gives the same result as above.

As mentioned above, the SI, SIS, SIR, SIRS, SEI, SEIS, SEIR models are all sub-models of the SEIRS model. Then, from the results obtained for the SEIRS model, we deduce that of the others. In Table 2, we have the values of the basic reproduction number $R_{0}$ and the Malthusian parameter $\alpha$ of the models listed above. The fourth column gives the changes to get the corresponding model. Further by setting $\sigma=0$, one gets the corresponding results for a disease without an additional death rate for the infectious.

Remark 3.6**.**

The basic reproduction number $R_{0}$ and the Malthusian parameter $\alpha$ are identical for the SEIR, SEIS and SEIRS models.

From Equations (3.2) and (3.3), we have

[TABLE]

It follows that the basic reproduction number $R_{0}$ exceeds $1$ if and only if, the Malthusian parameter $\alpha$ exceeds [math]. That is, the sign relation $sign(\alpha)=sign(R_{0}-1)$ is verified.

Remark 3.7**.**

It is well known that to surely prevent the disease to invade the population, $R_{0}$ must be less than $1$. To control the epidemic one need then to diminish $R_{0}$. The basic reproduction number $R_{0}$ can be written in the following form.

[TABLE]

So, it is clear that $R_{0}$ increases with $\kappa$ and $\nu$, and decreases with $\mu,\delta$ and $\sigma$. The contact rate $\kappa$ can be reduced by hospitalization or quarantine of the infectious individuals. The recovery rate $\delta$ can be increased by medication. In the case of an epizootic, $\sigma$ the disease related death rate, can be increased by culling infectious animals.

We have derived two thresholds (the Malthusian parameter $\alpha$ and the basic reproduction number $R_{0}$) of the SEIRS epidemic branching process, deduced the corresponding thresholds for the sub-models of the SEIRS epidemic model and established that $sign(\alpha)=sign(R_{0}-1)$.

3.4 Main result for the initial phase

In this subsection, we derive the probability for a minor outbreak of the epidemic branching process, and state the main result for the initial phase of the epidemic.

Let $\pi=P(\lim_{t}L_{\infty}(t)=0)$ be the probability of a minor outbreak of $L_{\infty}$ and $Y$ be the number of infectious contacts that an individual has during the infectious period. As we start with one latent individual, $\pi$ is the smallest positive solution of the equation $z=g(z)$, where $g$ is the probability generating function (pgf) of $Y$ [13, page 113]. We have

[TABLE]

Then, $\pi$ is the smallest solution in $[0,1]$ of the following equation.

[TABLE]

Equation (3.4) has two solutions,

[TABLE]

Then, we have the following theorem.

Theorem 3.8**.**

Let $\pi$ be the probability of a minor outbreak of the epidemic when started with one latent individual. Then,

[TABLE]

Corollary 3.9**.**

Let $\pi_{(m,k)}$ be the probability of a minor outbreak when the epidemic started with $m$ latent individuals and $k$ infectious individuals. Then,

[TABLE]

Proof.

Let $\pi_{(m,k)}$ be the probability of a minor outbreak when the epidemic starts with $m$ latent and $k$ infectious individuals. With this notation, we have

$\pi_{(1,0)}=P(\textrm{minor outbreak}|E(0)=1,I(0)=0)=\pi$, and

$\pi_{(0,1)}=P(\textrm{minor outbreak}|E(0)=0,I(0)=1)$. Thus, $\pi=\mu/(\mu+\nu)+(\nu/(\mu+\nu))\pi_{(0,1)}$.

Therefore, $\pi_{(0,1)}=((\mu+\nu)/\nu)\pi-\mu/\nu$. Thus, by Equation (3.8) one gets

[TABLE]

We have $\pi_{(m,k)}=\left(\pi_{(1,0)}\right)^{m}\left(\pi_{(0,1)}\right)^{k}$, since all the $m+k$ independent epidemics must die out [13, page 112]. Thus, by Equations (3.8) and (3.16), we get Equation (3.12). ∎

Remark 3.10**.**

This result is the same as that of the SEIR stochastic model studied by Allen and Lahodny in [1].

As noted above, due to the additive death rate in the infectious compartment, if the epidemic takes off, the population process can be turned to a sub-critical process. In this case, the population may go extinct. As we start the process with $n$ individuals, assuming that $n$ is large, we define the probability of minor outbreak of the epidemic, as the probability that the number of infected individuals does not exceed $\sqrt{n}$, that is $P(L_{n}(t)<\sqrt{n},\forall t\geq 0)$ [9, page 55]. Let us state now the main result for the initial phase.

Theorem 3.11**.**

*Consider the uniform SEIRS epidemic model defined above, with

$(S_{n}(0),E_{n}(0),I_{n}(0),R_{n}(0))=(n-1,1,0,0),L_{n}(t)=E_{n}(t)+I_{n}(t),N_{n}(t)=S_{n}(t)+E_{n}(t)+I_{n}(t)+R_{n}(t)$, and let $L_{\infty}(t)$ denote the birth and death process defined above, $\alpha$ its Malthusian parameter, $\pi$ the probability of a minor outbreak of $L_{\infty}$, and $\pi_{n}:=P(L_{n}(t)<\sqrt{n},\forall t\geq 0)$ denote the probability of a minor outbreak of the epidemic. Then as $n\rightarrow\infty$, we have the following results:*

1. i.

If $\alpha\leq 0$, then for any $n,L_{n}(t)\rightarrow 0$ as $t\rightarrow\infty$ with probability $1$. 2. ii.

If $0<\alpha<\lambda-\mu$, then $\pi_{n}\rightarrow\pi=\mu/(\mu+\nu)+\nu/[(\mu+\nu)R_{0}]$. With the remaining probability $(1-\pi_{n})\rightarrow 1-\pi$, $L_{n}$ grows exponentially: $L_{n}(t)\sim e^{\alpha t}$, but $L_{n}(t)/N_{n}(t)\rightarrow 0$ as $t\rightarrow\infty$. 3. iii.

If $\alpha>\lambda-\mu$ , then $\pi_{n}\rightarrow\pi=\mu/(\mu+\nu)+\nu/[(\mu+\nu)R_{0}]$. With the remaining probability $(1-\pi_{n})\rightarrow 1-\pi$, during the initial phase of the epidemic, $L_{n}$ grows exponentially with rate $\alpha$.

Proof.

1. i.

If $\alpha<0$, then $L_{\infty}$ is sub-critical and dies out with probability $1$. If $\alpha=0$, then $L_{\infty}$ is critical and dies out with probability $1$, since $P(Y=1)\neq 1$ [13]. 2. ii.

If $\alpha>0$, then $L_{\infty}$ is super-critical. It dies out with probability $\pi$ and with the remaining probability $1-\pi$, it grows exponentially with rate $\alpha$, that is $L_{\infty}\sim e^{\alpha t}$. As $N_{n}(t)\sim e^{(\lambda-\mu)t}$, if $\alpha<\lambda-\mu$, then $L_{n}(t)/N_{n}(t)\longrightarrow 0$, when $t\longrightarrow\infty$. 3. iii.

If $\alpha>\lambda-\mu$, then $\alpha>0$ since $\lambda-\mu>0$. Thus, $L_{\infty}$ is super-critical. It dies out with probability $\pi$ and with the remaining probability $1-\pi$, it grows exponentially with rate $\alpha$. As $L_{n}\Longrightarrow L_{\infty}$, the same applies to $L_{n}$.

∎

Remark 3.12**.**

If $\alpha>\lambda-\mu$, then $L_{\infty}$ is super-critical. Thus the epidemic may take off. If it does, the number of infected individuals grows initially with the rate $\alpha$ that is larger than the initial growth rate of the population ($\lambda-\mu$). After the initial phase, in the case of a major outbreak, several scenarios are possible. i)The population goes on growing exponentially, eventually with a lower rate; ii) due to the additive death rate of the infectious, the effective death rate of the population becomes larger than its birth rate, and the population process becomes a sub-critical process. The different scenarios are treated in Section 4 and illustrated by simulations in Section 5.

We have derived the probability of a minor outbreak of the epidemic branching process, stated and shown the main result of the initial phase of the epidemic.

4 The deterministic SEIRS model

Now we consider the corresponding deterministic model of the stochastic model studied above. As the population size is varying, we study first the fractions system and then deduce the asymptotic behavior of the compartments sizes. In this section the deterministic sizes of the compartments S, E, I, R and the population size at the time $t$ are denoted $S(t),E(t),I(t),R(t)$ and $N(t)$ respectively.

The corresponding deterministic SEIRS model of the model above is given by the following system of ordinary differential equations (ODE).

[TABLE]

Remark 4.1**.**

System (4.1) is the same as System (2) studied by Greenhalgh in [11], with a constant contact rate ($\beta(N)=\beta$), a constant death rate ($f(N)=\mu$) and without vaccination ($p=q=0$). But Greenhalgh assumed that the death rate $f(N)$ is a strictly monotone increasing continuously differentiable function of N. So, the model that we study is not a sub-model of that of Greenhalgh since we consider a constant death rate.

From System (4.1), we have $dN/dt=(\lambda-\mu)N-\sigma I=(\lambda-\mu-\sigma i)N$, where $i$ is the fraction of the infectious individuals. Thus, the population should grow if $\lambda>\mu+\sigma i$, stabilize if $\lambda=\mu+\sigma i$ and decrease if $\lambda<\mu+\sigma i$.

Theorem 4.2**.**

$N(t)$* is constant and positive ($N(t)=N(0)>0,\forall t>0$), if and only if the parameters verify the following equality:*

[TABLE]

and the initial values verify

[TABLE]

Proof.

By using successively the derivatives of $N,I,E,R$ and $S$, one gets that $N(t)$ is constant and positive, if and only if $I,E,R$ and $S$ are constant, the parameters verify Equation (4.2) and the initial values verify System (4.3). ∎

Generally, Equation (4.2) and System (4.3) are not verified, thus $N(t)$ is not constant. Therefore, we consider the fractions $s=S/N,e=E/N,i=I/N$ and $r=R/N$. By System (4.1), one gets

[TABLE]

This system is equivalent to System (2.2) in [11], with $p=q=0$.

Remark 4.3**.**

The natural death rate $\mu$ does not intervene in the derivatives of the fractions. This is logic, since this rate is the same for all the compartments, it has no effect on the fractions.

Since $r=1-s-e-i$, it is enough to study the system

[TABLE]

in the domain

[TABLE]

Theorem 4.4**.**

The domain D is positively invariant for System 4.5.

Proof.

If $s=0$, then $ds/dt=\lambda+\rho(1-e-i)>0$. If $e=0$, then $de/dt=\kappa si\geq 0$. If $i=0$, then $di/dt=\nu e\geq 0$. If $s+e+i=1$, then $d(s+e+i)/dt=-\delta i\leq 0$. Then every solution of System 4.5 starting in $D$, remains there for all $t>0$. That is $D$ is positively invariant for System 4.5 ∎

In the following, we use the next generation matrix (NGM) described in [10] to derive a threshold parameter, with threshold value $1$ for System (4.5).

By setting $ds/dt=de/dt=di/dt=0$ with $e=i=0$ in System (4.5), we get $s=1$. Then, $(1,0,0)$ is the unique disease free equilibrium (DFE) of System (4.5). $e$ and $i$ are the infected fractions of the model. Thus, the infected subsystem is

[TABLE]

Hence, the linearized infected subsystem at the DFE is

[TABLE]

Let $x=(e,i)^{t}$ be the vector of the infected fractions. Thus, System (4.7) is equivalent to

[TABLE]

with

[TABLE]

$T$ is the transmissions matrix and $\Sigma$ is the transitions matrix. Then, the next generation matrix with large domain is

[TABLE]

Let $R_{1}$ be the spectral radius of $K_{L}$. We have

[TABLE]

An equilibrium is said to be stable if nearby solutions stay nearby for all future time [14, p. 175]. More precisely an equilibrium $x^{*}$ is said to be stable, if for every neighborhood $V$ of $x^{*}$ there is a neighborhood $V_{1}$ of $x^{*}$, such that every solution starting in $V_{1}$ remains in $V$ for all $t>0$. If $V_{1}$ can be chosen such that $\lim_{t\rightarrow\infty}x(t)=x^{*}$, then $x^{*}$ is said to be asymptotically stable. An equilibrium is said to be unstable, when it is not stable. An equilibrium $x^{*}$ is said to be globally asymptotically stable (GAS) in an invariant set $D$, $(x^{*}\in D)$, if it is locally stable and $\lim_{t\rightarrow\infty}x(t)=x^{*}$, for every solution $x(t)$ starting in $D$.

Theorem 4.5**.**

The disease free equilibrium (DFE) of System (4.5) is globally asymptotically stable (GAS) in $D$, if $R_{1}\leq 1$, and unstable if $R_{1}>1$.

Proof.

Let $f(s,e,i)$ be the Rhs of System (4.5). Then,

[TABLE]

The Jacobian of $f$ at the disease free equilibrium is

[TABLE]

The characteristic polynomial of $Df(1,0,0)$ is

[TABLE]

$-(\lambda+\rho)$ is an evident negative root of $P(x)$. Thus, by the Routh-Hurwitz criterion [20, page 11], all the roots of $P(x)$ has negative real part if and only if $(\lambda+\nu)(\lambda+\delta+\sigma)-\nu\kappa>0$. And we have $(\lambda+\nu)(\lambda+\delta+\sigma)-\nu\kappa>0\Longleftrightarrow R_{1}<1$. Thus, the disease free equilibrium is locally asymptotically stable if $R_{1}<1$, and unstable if $R_{1}>1$.

Let $V$ denote the function defined on $D$ by $V(s,e,i)=\nu e+(\lambda+\nu)i$. Then,

[TABLE]

with $L(s,e,i)=\nu\kappa s+\nu\sigma e+(\lambda+\nu)\sigma i-(\lambda+\nu)(\lambda+\delta+\sigma)$. The affinity of $L$ implies that it achieves its maximum at the extreme points of the boundary of the closed set $D$. But $L(0,0,0)=-(\lambda+\nu)(\lambda+\delta+\sigma),L(0,0,1)=-(\lambda+\nu)(\lambda+\delta),L(0,1,0)=-\lambda\sigma-(\lambda+\nu)(\lambda+\delta)$ and $L(1,0,0)=\nu\kappa-(\lambda+\nu)(\lambda+\delta+\sigma)=(\lambda+\nu)(\lambda+\delta+\sigma)(R_{1}-1)$. Thus, $\stackrel{{\scriptstyle.}}{{V}}\leq 0$ in $D$ if $R_{1}\leq 1$. Then, $V$ is a Lyapunov function of System (4.5). The only invariant subset of the set with $\stackrel{{\scriptstyle.}}{{V}}=0$ is $\{(1,0,0)\}$. It follows from LaSalle’s Invariance Principle [14, p. 200], that the disease free equilibrium (DFE) is globally asymptotically stable (GAS) in $D$, when $R_{1}\leq 1$. ∎

By Equations (4.8) and (3.2), we have

[TABLE]

It follows that the fraction’s threshold $R_{1}$ exceeds $1$ if and only if the Malthusian parameter $\alpha$, exceeds the initial growth rate of the population $\lambda-\mu$. That is, we have the sign relation $sign(R_{1}-1)=sign(\alpha-(\lambda-\mu))$. Thus, the global stability of the disease free equilibrium of the fraction’s system when $R_{1}<1$, confirms that if $\alpha<\lambda-\mu$, then the infected fraction vanishes even if the epidemic takes off (Theorem (3.11) (ii)).

Greenhalgh has shown [11, Theorem 2.3] that if $R_{1}>1$, then System (4.5) has at least one endemic equilibrium, and that this equilibrium is unique and locally asymptotically stable (LAS) when the average duration of immunity exceeds both the average infectious and incubation periods, that is $\delta>\rho$ and $\nu>\rho$. We have not proved, but we strongly believe that if $R_{1}>1$, then System (4.5) has one and only one endemic equilibrium, and that this equilibrium is globally asymptotically stable in the interior of $D$. The simulations that we made support this conjecture (Figure 4 (c) and (d)).

Theorem 4.6**.**

Let $(S(t),E(t),I(t),R(t))$ be a solution of System (4.1) and $R_{0}$ denoted the basic reproduction number given by Equation (3.3).

1. i.

If $R_{0}<1$, then $(S(t),E(t),I(t),R(t))\longrightarrow(\infty,0,0,0)$; 2. ii.

if $R_{0}=1$, then $(S(t),E(t),I(t),R(t))\longrightarrow(\infty,E^{*},I^{*},R^{*})$,

with $E^{*}>0,I^{*}>0,R^{*}>0$; 3. iii.

if $R_{0}>1\geq R_{1}$, then $(S(t),E(t),I(t),R(t))\longrightarrow(\infty,\infty,\infty,\infty)$.

Remark 4.7**.**

The case $R_{1}>1$ is treated in Theorem 4.9.

Proof.

We have

[TABLE]

Then, $R_{0}>R_{1}$, since $\lambda>\mu$. Therefore, in the three cases of the Theorem 4.6, we have $R_{1}\leq 1$. Let us assume that $R_{1}\leq 1$. Thus, by Theorem (4.5), $(s,e,i,r)\longrightarrow(1,0,0,0)$ when $t\longrightarrow\infty$. $dN/dt=(\lambda-\mu)N-\sigma I=(\lambda-\mu-\sigma i)N$. Then, $dN/dt\longrightarrow(\lambda-\mu)N$, when $t\longrightarrow\infty$. Thus, $N\longrightarrow\infty$, when $t\longrightarrow\infty$ because $\lambda>\mu$. Therefore, $S\longrightarrow\infty$, when $t\longrightarrow\infty$, since $S/N\longrightarrow 1$. By using the derivatives of $E$ and $I$, one gets

[TABLE]

[TABLE]

Thus $E/I$ can be approximate by a solution of the following equation, when $t\longrightarrow\infty$.

[TABLE]

Equation (4.9) is a Riccati’s equation [12]. By solving it, one gets

[TABLE]

Then, $E/I\longrightarrow(\delta+\sigma-\nu+\sqrt{\Delta})/(2\nu)$, when $t\longrightarrow\infty$.

We have $dI/dt=\nu E-(\delta+\mu+\sigma)I=[\nu(E/I)-(\delta+\mu+\sigma)]I$. Thus, by substituting $E/I$ by its asymptotic value, one gets

[TABLE]

and

[TABLE]

Therefore, $dI/dt\longrightarrow\alpha I$, when $t\longrightarrow\infty$.

As $dE/dt=\kappa SI/N-(\nu+\mu)E$, $S/N\longrightarrow 1$ and $E/I\longrightarrow(\delta+\sigma-\nu+\sqrt{\Delta})/(2\nu)$, after some algebra, one gets $dE/dt\longrightarrow\alpha E$, when $t\longrightarrow\infty$. As $sign(\alpha)=sign(R_{0}-1)$,

[TABLE]

For the number of the recovered $R$, as $dR/dt=\delta I-(\rho+\mu)R$, it is obvious that $R$ has the same asymptotic behavior as $I$. ∎

Remark 4.8**.**

In the proof, we have shown that if $R_{1}\leq 1$, then the Malthusian parameter $\alpha$ of the stochastic model, is also the common asymptotic growth rate of the compartments E and I, and $\lambda-\mu$ is the asymptotic growth rate of the population. Since $sign(R_{1}-1)=sign(\alpha-(\lambda-\mu))$, this is coherent with Theorem 3.11 (ii).

Theorem 4.6 gives the asymptotic behavior of the compartments sizes, when $R_{1}\leq 1$. If $R_{1}>1$, then the fraction disease free equilibrium is unstable. Therefore, the disease will remain endemic in the population in term of the fraction infected. The following theorem gives the asymptotic behavior of the compartments sizes, when $R_{1}>1$, assuming that the fraction system admits an endemic equilibrium that is globally asymptotically stable in the interior of the feasible region $D$.

Theorem 4.9**.**

Assume that $R_{1}>1$ and that System (4.5) has a unique endemic equilibrium $(s^{*},e^{*},i^{*})$ that is globally asymptotically stable in $\stackrel{{\scriptstyle o}}{{D}}$, and set

[TABLE]

1. i.

If $R_{2}>1$, then $(S,E,I,R)\longrightarrow(\infty,\infty,\infty,\infty)$. 2. ii.

*If $R_{2}=1$, then $(S,E,I,R)\longrightarrow(S^{*},E^{*},I^{*},R^{*})$, *

with $S^{*}>0,E^{*}>0,I^{*}>0$, $R^{*}>0$. 3. iii.

If $R_{2}<1$, then $(S,E,I,R)\longrightarrow(0,0,0,0)$.

Proof.

Let us assume that there is an endemic equilibrium $(s^{*},e^{*},i^{*})$ of the fraction System (4.5) and that it is globally asymptotically stable in $\stackrel{{\scriptstyle o}}{{D}}$. Then,

[TABLE]

Asymptotically, the population would increase with rate $\lambda$, and decrease with rate $\mu+\sigma i^{*}$. As the fraction system admits an endemic equilibrium, that is globally asymptotically stable in the interior of the feasible region, all the compartments have the same asymptotic behavior as the population. Let us set $\alpha_{2}=\lambda-\mu-\sigma i^{*}$. The quantity $\alpha_{2}$ is the common asymptotic exponential growth/decay rate of all the compartments S, E, I and R. We have $sign(\alpha_{2})=sign(R_{2}-1)$. Thus, the results follow. ∎

In this section we have studied the corresponding deterministic SEIRS model of the previous stochastic model. We derived a threshold quantity $R_{1}$ for the fraction model. If $R_{1}\leq 1$, then the fraction’s disease free equilibrium is globally asymptotically stable in the feasible region $D$, otherwise it is unstable. When $R_{1}\leq 1$, the behavior of the number of infected is determined by the basic reproduction number $R_{0}$. If $R_{0}<1$, then the number of infected vanishes. If $R_{0}=1$, then the number of infected stabilizes to a positive value; when $R_{1}\leq 1<R_{0}$, then the number of infected grows exponentially, but at a lower rate than the population. If $R_{1}>1$, then the number of infected grows initially quicker than the population and the asymptotic behavior of the population is governed by the threshold quantity $R_{2}$. If $R_{2}<1$, then the population vanishes; if $R_{2}=1$, then the population stabilizes; if $R_{2}>1$, then the population grows, but with a lower rate than its initial growth rate.

5 Simulations

In this section, we use the software R to illustrate and confirm the results found in the previous sections. In the following, we set $\mu=1$, that is the time unit is the life expectancy, except for the simulations of influenza epidemics in Burkina Faso where we set one year as the time unit. The other parameters and the initial values are chosen arbitrary, unless otherwise stated.

5.1 Simulations of the initial phase

In this subsection, we give some examples of simulations of epidemics starting by one latent individual, and using different values of the parameters. The population is initiated with $1$ latent individual and $999$ susceptible individuals.

In Figure 2 (a) and (b), where $R_{0}=0.41$ and $R_{0}=0.73$ respectively, all the $10$ epidemics die without any major outbreak. In (a) the maximum of infected individuals is $2$, while it is $6$ in (b).

In Figure 3, where $R_{0}=2$, four simulated epidemics out of $10$ have a major outbreak. The other $6$ epidemics die out without many getting infected. For the epidemics with major outbreak, the number of the infected individuals grow exponentially but the time where the exponential growth starts varies.

Now we estimate the probability of a minor outbreak $\pi$ by simulating $1000$ epidemics and setting $\hat{\pi}_{n}=n_{0}/1000$, where $n_{0}$ is the number of minor epidemics. We set $\lambda=3,\mu=1,\nu=50,\delta=10,\sigma=4,\rho=3$, and set successively $\kappa=10,20,30,50,100$ to get different values of $\pi$. Table 3 gives the different values of $\pi$ and the estimate $\hat{\pi}_{n}$ for $n=1000$ and $n=2000$ respectively. These results confirm that the probability of extinction of the branching process $L_{\infty}$ is a good approximation of the probability of a minor outbreak of the epidemic starting with one latent individual in a population of size $n$, when $n$ is large.

These simulations confirm that when $R_{0}$ is less than $1$, the disease cannot invade the population, and that if $R_{0}$ is larger than $1$, then with a positive probability $(1-\pi_{n})\longrightarrow(1-\pi)$, the disease can invade the population, and in this case the number of the infected individuals grows exponentially during the initial phase.

5.2 Simulations of major outbreaks

In this subsection, we show some simulations of major outbreaks, where the epidemic starts with a positive fraction of infected individuals. So the simulations illustrate what happens once the number of infected has reached a small but positive fraction of the community. We use the blue color for the susceptible, green for the exposed, red for the infectious and black for the recovered.

We start by simulating the deterministic fraction’s system. In Figure 4, we have four cases with $R_{1}=0.5,1,1.76,5.33$ respectively. In each case we have ten solutions paths of System (4.4) with different initial values. In (a) as in (b), the $10$ solutions of System (4.4) approach the disease free equilibrium, confirming that if $R_{1}\leq 1$, then the disease free equilibrium of the deterministic fraction system is globally asymptotically stable in the feasible region. In (c), all the $10$ solutions approach the same endemic equilibrium $(s^{*},e^{*},i^{*},r^{*})\approx(0.50,0.13,0.14,0.24)$. In (d), all the $10$ solutions approach the same endemic equilibrium $(s^{*},e^{*},i^{*},r^{*})\approx(0.12,0.49,0.31,0.08)$. The results of (c) and (d) confirm that when $R_{1}>1$, the disease free equilibrium is unstable and that there is an endemic equilibrium that is globally asymptotically stable in the interior of the feasible region. That is the endemic level is independent of the starting value. However the endemic equilibrium of (d) is different of that of (c), hence the endemic equilibrium vary with the parameters values.

In the following we simulate both the stochastic and the deterministic models. In each case, we have simulated the stochastic epidemic, as well as integrated numerically the deterministic system (4.1), both starting at the same values. From the numbers, we got the fractions by setting $s=S/N,e=E/N,i=I/N,r=R/N,\textrm{ with }N=S+E+I+R$. One distinguishes the stochastic solutions from the deterministic by the fact that the deterministic solutions are represented by smooth lines, while the solutions of the stochastic solutions are represented by broken lines.

In Figure 5, where $R_{0}=1.27>1$ and $R_{1}=0.95<1$, the numbers of latent, infectious and recovered grow exponentially as the population, but the population growth rate is even larger and the fractions of the infected compartments go to zero. This means that, in term of the number of infected individuals, the disease is endemic, but the disease dies out in term of the fractions.

In Figure 6, initially the number of the infected grows quicker than the number of the susceptible. After that, the number of the susceptible declines. Afterward all the compartments grow exponentially, but with a rate $\alpha_{2}\approx 0.16$ that is lower than the initial growth rate of the population $\lambda-\mu=1$. The fractions go to an endemic equilibrium.

For Figure 7, the parameters are chosen such that Equation (4.2) is verified, allowing then the existence of endemic equilibrium for the deterministic System 4.1. The asymptotic reproduction number of the population $R_{2}=1$, then for the deterministic solution, the population stabilizes, when $t\longrightarrow\infty$. The stochastic numbers fluctuate around the deterministic numbers. The fractions have the same behavior as that of the numbers.

In Figure 8, for the deterministic model, the epidemic turned the population exponential growth to an exponential decay due to the disease induced death rate, while the fractions go to an endemic equilibrium. For the stochastic model, the population size first decreases but at some instance when there are only few remaining individuals, the disease goes extinct and then the population size starts growing again. The deterministic model suggests that the population will go extinct, whereas in the stochastic model the disease first dies out and then the population becomes super critical again, thus regrowing. In the stochastic setting, what happens when the numbers become low is random. Both the disease and the population could die out, or the population can grow again after the extinction of the epidemic. Only analyzing the deterministic fraction model would give a misleading conclusion since the fractions seem to stabilize, whereas in what really happens is that all numbers in the deterministic model tend to [math].

In Figure 9, we have the case where $R_{0}=1$. For the stochastic model the disease goes extinct, while it persists in the deterministic one. In both models the population goes on growing exponentially. In (b) we made a zoom to see the dynamics of $E,I$ and $R$.

These simulations show the different possible scenarios in the case of a major outbreak. They confirm the theoretical results and show the similarities and differences between the stochastic model and the deterministic model.

5.3 Simulation of influenza epidemics in Burkina Faso

Now we simulate two influenza epidemics in Burkina Faso with two different basic reproduction numbers. Burkina Faso is an inland country of West Africa. Its population in $2016$ is about $16~{}000~{}000$. The individual annual birth rate and death rate are estimated to $\lambda=0.046$ and $\mu=0.0118$ respectively [15].

According to the World Health Organization (WHO) [21], influenza is caused by a virus that attacks mainly the respiratory tract: the nose, the throat, the bronchi and rarely also the lungs. The infection usually lasts for about a week. It is characterized by sudden onset of high fever, myalgia, headache and severe malaise, non-productive cough, sore throat, and rhinitis. Most people recover within one to two weeks without requiring any medical treatment. The virus is easily passed from person to person through the air by droplets and small particles excreted when infected individuals cough or sneeze. The influenza virus enters the body through the nose or the throat. It then takes between one and four days for the person to develop symptoms. Someone suffering from influenza can be infectious from the day before he/she develops symptoms until seven days afterwards. The disease spreads very quickly among the population especially in crowded circumstances. Cold and dry weather enables the virus to survive longer outside the body than in other conditions and, as a consequence, seasonal epidemics in temperate areas appear in winter. Much less is known about the impact of influenza in the developing world. However, influenza outbreaks in the tropics where viral transmission normally continues year-round tend to have high attack and case-fatality rates. Therefore by setting one year as the time unit, we make the following estimation for the influenza parameters. The latent period is approximately $2.5$ days, the infectious period is approximately $7$ days and the immunity period is approximately one year, thus $\nu=365/2.5$, $\delta=365/7$, $\rho=1$. We set the influenza case fatality rate (CFR) to $0.1\%$ and deduce the influenza related death rate $\sigma\approx 0.0522$.

The reproduction number of the $1918$ pandemic influenza is estimated to be between $2$ and $3$ [19]. Thus we set $R_{0}=2.5$ and deduce the contact number $\kappa=130$ from the other parameters and Equation (3.3). We simulate the epidemic for a period of $10$ years starting in $2016$ by integrating numerically the deterministic System (4.1). Figure 10 gives the evolution of the numbers of susceptible, latent, infectious and recovered individuals and the corresponding fractions during the $10$ years. Figures 10 (c) and (d), show respectively the dynamics of $E$ and $I$ and that of the fractions $e$ and $i$. We have a peak with $2~{}783~{}834$ infectious individuals at the $11^{th}$ week of the epidemic. After that the number of infectious declines because the number of susceptible is low. Afterward we have a minor peak every year due to the immunity waning and the newborns that increase the number of susceptible. The number of recovered individuals grow quickly and reach its maximum $13~{}139~{}592$ at the $17^{th}$ week. The fractions go to an endemic equilibrium through damped oscillations. The population in $2026$ is estimated to $22~{}370~{}000$ individuals with $8~{}960~{}000$ susceptible, $94~{}000$ latent individuals, $260~{}000$ infectious individuals and $13~{}060~{}000$ recovered with non-permanent immunity. The number of recovered individuals is larger than the number of susceptible individuals. The fractions go to an endemic equilibrium with more than $1\%$ of the population infected at every time. As the influenza last about one week and we have $52$ weeks within a year, more than $50\%$ of the population should be infected during the year $2026$. That will have a very important negative impact on the economy of the country. We have $R_{1}\approx 2.50$ and $R_{2}\approx 3.71$, thus $R_{1}$ and $R_{2}$ are both larger than $1$ and according to Theorem 4.9 the fractions should go to an endemic equilibrium and all the compartments should grow exponentially. Therefore the simulations agree with the theoretical results.

The basic reproduction number for the novel influenza A (H1N1) has been estimated to be between $1.4$ and $1.6$ [8]. Thus, we set now $R_{0}=1.5$ and deduce the contact number from the other parameters. The results of this simulation are shown in Figure 11. In this case the epidemic and the population have globally the same dynamics as above. But the impact of the epidemic is fewer. The major peak of the epidemic happens later, at the $25^{th}$ week, with $794~{}004$ infectious individuals. Contrary to the preceding case, the number of the recovered individuals is below the number of the susceptible individuals.

In [8], Coburn, Wagner and Blower simulated an influenza epidemic using a SIR model with demography. Their result for the number of infectious individuals [8, Figure 2 (a)] is similar to that of Figure 10 (c) and Figure 11 (c). We assume no seasonal effects. Adding seasonality should make seasonal effects remain [8, Figure 2 (b)]. The simulations of influenza in Burkina Faso show that in spite of the epidemic the population should go on growing. But the number of infected individuals will grow also. Furthermore the peak of the epidemic in the first year show that the emergence of a new strain of influenza virus will be a very serious threat for the world. The Global Influenza Program (GIP) of the World Health Organization (WHO), provides Member States with strategic guidance, technical support and coordination of activities essential to make their health systems better prepared against this threat.

In this section we have illustrated and validated the theoretical results of the previous sections by simulations.

6 Conclusion and discussions

In this paper we have studied a stochastic SEIRS epidemic model, with a disease related death, in a population which grows exponentially without the disease. We assumed that initially, the population process is a super-critical linear birth and death process with birth rate $\lambda$ and death rate $\mu$.

We have derived, the basic reproduction number $R_{0}$, the Malthusian parameter $\alpha$ and the probability of minor outbreak $\pi$ assuming that the initial size $n$ of the population tends to infinity. Considering the deterministic model, we derived the threshold quantity $R_{1}$ for the fractions.

If $R_{0}\leq 1$, then the disease dies out with probability $1$. That is, there is no possibility of major outbreak if $R_{0}\leq 1$. In this case, the population remains a super-critical process.

If $R_{0}>1$ then, with a positive probability, the epidemic can take off. If the epidemic takes off, then the number of the infected (exposed or infectious) individuals grows exponentially with rate $\alpha$. If $0<\alpha\leq\lambda-\mu$, then the sizes of all the compartments grow exponentially while the fraction of the infected individuals goes to [math]. The number of infected people grows, but at a lower rate than the population, implying that the fraction infected becomes negligible. If $\alpha>\lambda-\mu$, then the number of infected individuals grows initially with a rate that is larger than the population growth rate, and different scenarios are possible. Due to the additional death rate $\sigma$ in the infectious compartment, the population will go on growing but with a lower rate, or the population will become a sub-critical branching process and thus have a decreasing size. In the latter case the population vanishes in the deterministic model while in the stochastic one, what happens when the numbers become low is random. Both the disease and the population could die out, or the population could grow again after the extinction of the epidemic.

We have illustrated and validated the theoretical results by simulations. These simulations show the similarities and the differences between the stochastic model and the corresponding deterministic model. If $R_{0}>1$, then in the deterministic model, the epidemic will invade the population surely; while in the stochastic model, with a positive probability $\pi$ the disease vanishes. When $R_{2}<1$, the population vanishes surely in the deterministic model; whereas in the stochastic model, it can regrow exponentially after the extinction of the epidemic. One need to remember that an epidemic is always a stochastic process, and the deterministic model fits only when we have a large community with a large number of infectious individuals [6].

For some diseases (e.g. influenza, measles), the susceptibility and the infectiousness vary with the age of the individuals. Therefore it is more realistic to consider an heterogeneous population as in [18]. Furthermore adding seasonal forcing will increase realism for seasonal diseases. For other diseases like sexually transmitted diseases (STD), a close and/or long contact is required for transmission. Then the dynamics of the epidemic are linked to the special network in the host population [4, 2, 18]. Due to the development of the migration of populations, an epidemic starting in one location can be exported to another one quickly, then a meta-population model is convenient to find the conditions for a global health security [3, 7]. Under the leadership of the World Health Organization (WHO) many policies (vaccination, medication, quarantine, etc.) are implemented to prevent major outbreak of epidemics. Nevertheless the infectious diseases remain a serious threat for Humanity. Another step towards realism is hence to extend this model by adding vaccination [6, 11] and treatment [17]. This should allow to find the optimal control strategy to realize the herd immunity, that is to prevent the major outbreaks of infectious diseases.

Acknowledgments

We are grateful to the International Science Program (ISP) of Uppsala University and the Swedish International Development Agency (SIDA) for their financial support.