SEIRS epidemics in growing populations
SSEIRS epidemics in growing populations
Tom Britton a , Désiré Ouédraogo b, ∗ a Department of Mathematics, Stockholm University, SE-106 91 Stockholm, Sweden b Laboratoire de Mathématiques et Informatique (LAMI), EDST, Université Ouaga I Pr. JosephKi-Zerbo 03 B.P.7021 Ouagadougou 03, Burkina Faso
Abstract
An SEIRS epidemic with disease fatalities is introduced in a growing population (modelledas a super-critical linear birth and death process). The study of the initial phase ofthe 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) Thedisease dies out quickly, only infecting few; ii) the epidemic takes off, the number ofinfected individuals grows exponentially, but the fraction of infected individuals remainsnegligible; iii) the epidemic takes off, the number of infected grows initially quicker thanthe population, the disease fatalities diminish the growth rate of the population, but itremains super critical, and the fraction of infected go to an endemic equilibrium; iv) theepidemic takes off, the number of infected individuals grows initially quicker than thepopulation, the diseases fatalities turn the exponential growth of the population to anexponential decay.
Keywords:
SEIRS epidemic, threshold quantities, initial growth, endemic level
1. Introduction
Infectious diseases remain a threat for developing countries as well as for developed coun-tries. Many mathematicians focus their efforts to understand the dynamics of infectiousdiseases, in order to find the conditions to eradicate them. In mathematical modellingof infectious disease epidemics, the population in which the disease is spreading is parti-tioned in several compartments according to the status of the individuals, related to thedisease. Every epidemic model has at least, the compartment I of the infectious indi-viduals who are infected and able to transmit the disease to others through contact, andthe compartment S of the susceptible individuals (those who are not infected but maybe infected if they contact an infectious individual). Two other compartments often usedare the compartment E of the exposed or latent individuals who are already infected butnot yet able to transmit the disease to others, and the compartment R of the recoveredor removed individuals (those who are healed from the disease with a permanent or non-permanent immunity). In a SEIR epidemic, a susceptible individual infected througha contact with an infectious, becomes infected and latent; at the end of the latent pe-riod he/she becomes infectious and at the end of the infectious period he/she recovers ∗ Corresponding author
Email addresses: [email protected] (Tom Britton), [email protected] (Désiré Ouédraogo)
Preprint submitted to Elsevier September 27, 2018 a r X i v : . [ q - b i o . P E ] M a r ith a life-long immunity. An SEIRS epidemic is almost the same as the preceding, theonly difference is that a recovered individual has a non-permanent immunity (He/she canlose his immunity, becoming susceptible again). Diphtheria, influenza and pneumonia areexamples of diseases with latent period and non-permanent immunity [11].In [5], Britton and Trapman studied stochastic SIR and SEIR models in a growing popu-lation. They derived the basic reproduction number and the Malthusian parameter of theepidemic, stated results for the initial phase and showed that the stochastic proportionsprocess converges to a deterministic process.In [11], Greenhalgh studied an SEIRS deterministic model with vaccination and foundthat under some conditions, the solution has Hopf bifurcations.The aim of this paper is to study the dynamic of a stochastic SEIRS epidemic model withdisease induced mortality, in an exponentially growing population. As in [5] , we assumethat without the disease, the population has a birth rate λ , and a natural death rate µ ,such that λ > µ . That is, initially the population process is a super-critical linear birthand death process. An SEIRS epidemic is introduced by infecting one individual. Withthe disease, the population is divided in four compartments according to the status of theindividuals, related to the disease. The compartment S of the susceptible individuals, thecompartment E of the latent or exposed individuals , the compartment I of the infectiousindividuals , and the compartment R of the removed individuals (those who are healedof the disease with a non-permanent immunity). The process is initiated by setting ( S (0) , E (0) , I (0) , R (0)) = ( n − , , , . The transfer diagram of the model is given byFigure 1.We derive the Malthusian parameter α , the basic reproduction number R and the prob-ability of minor outbreak π of the epidemic. If R ≤ , then the disease cannot invadethe population, that remains a super critical process. If R > , then the epidemic hasa positive probability − π of taking off, with the remaining probability π , it dies out.If the epidemic takes off, another threshold parameter R determines the behavior of theproportion of infected individuals. If R ≤ , then the fraction of infected stays small;while it persists when R > . If R > , or equivalently α > λ − µ , then the number of the infected grows initially quicker than the population, the disease fatalities diminishthe growth rate of the population. In this case the asymptotic behavior of the populationrely on a third threshold quantity R . If R > , then the population goes on growing,while it becomes a sub-critical process when R ≤ . In the latter case, when the numberof individuals become low, the population should vanish with the disease, or regrows afterthe extinction of the epidemic.We start by defining the stochastic model in Section . Then, in Section , we study theinitial phase of the epidemic, thereafter we consider the deterministic model in Section . Afterward, we give some illustrations by simulating different scenarios of epidemics inSection . In Section , we conclude the paper and discuss some perspectives.
2. The model
Initially (before the introduction of the disease), the population model is a linear birthand death (B-D) process with individual birth rate λ and individual death rate µ . We2 igure 1: The transfer diagram of the SEIRS model assume that λ > µ , that is the process is super-critical. N ( t ) denotes the number ofindividuals in the population at time t . Now, we define a uniformly mixing Markovian epidemic model on the population describedabove, implying that individuals give birth at rate λ and die from other causes at rate µ ,irrespective of disease status. Initially (at t = 0 ), the population consists of n susceptibleindividuals. At this time, an SEIRS infectious disease is introduced by infecting oneindividual. The disease spreading is modelled as follow. An infected individual remainslatent (infected but not yet infectious) for exponential time with rate ν . After this period,he/she becomes infectious (unless he/she dies before). An infectious individual remainsinfectious for an exponential time with rate δ (unless he/she dies before). The diseaseinduces an additional death rate σ for the infectious individuals. After the period ofinfectiousness, an infectious individual recovers with a temporary immunity. A recoveredperson remains immune for an exponential time with rate ρ unless he/she dies before.After the period of immunity he/she becomes susceptible again. During the infectiousperiod, the infective has infectious contacts randomly in time according to a homogeneousPoisson process with rate κ , each time with a uniformly selected random individual. If thecontacted 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 newborn 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). Thepopulation is initiated at Z n (0) = ( S n (0) , E n (0) , I n (0) , R n (0)) = ( n − , , , . 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 (cid:28) n ).In short, the population model has two parameters, the birth rate λ and the death rate3 ; and the disease model has five parameters, the transmission rate κ , the end of latencyrate ν , the recovery rate δ , the disease death rate σ and the immunity waning rate ρ . Thepossible 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 modeldefined above. If ρ = 0 , then an infected individual cannot go back to the susceptiblestate and we get an SEIR model. If ρ −→ ∞ , then the recovered state vanishes and weget an SEIS model. If ν −→ ∞ , then the latent state vanishes giving an SIRS model.If δ = 0 , then we have lifelong infectivity and hence an SEI model. If ν −→ ∞ and ρ −→ ∞ , then the latent state and the recovered state vanish and we have an SIS model.If ν −→ ∞ and δ = ρ = 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 deducethose of the others. Event State change l RateBirth (1 , , , λN Death of susceptible ( − , , , µS Death of exposed (0 , − , , µE Death of infective (0 , , − ,
0) ( µ + σ ) I Death of recovered (0 , , , − µR Infection ( − , , , κSI/N End of latency (0 , − , , νE Recovery (0 , , − , δI Loss of immunity (1 , , , − ρR Table 1: The uniform Markovian dynamic SEIRS epidemic: type of events, their state change l (The oldstate z = ( u, v, x, y ) is hence changed to z + l ) and their rates. In this section we have presented the stochastic SEIRS model studied in this paper.
3. Results for the initial phase
Without the disease, the population process is a linear super-critical birth and deathprocess with individual birth rate λ and individual death rate µ . But, with the introduc-tion of the disease that induces an extra death rate σ for the infectious individuals, wehave two possible events. Birth with rate λN ( t ) and death with rate µN ( t ) + σI ( t ) =[ µ + σI ( t ) /N ( t )] N ( t ) , where N ( t ) is the current total number of individuals in the popula-tion and I ( t ) the current number of infectious individuals. Then, the population processis no longer a linear birth and death process (unless σ = 0 ). However, at the beginning ofthe epidemic when the fraction of the infectious individuals is very small ( I ( t ) /N ( t ) ≈ ),the population will behave almost as a linear birth and death process with birth rate λ and death rate µ . Thus for the initial phase of the epidemic, we assume that the popu-lation size is a linear super-critical birth and death process with birth rate λ and deathrate µ . On the other hand, if the epidemic takes off after the initial phase, whenever the fraction of the infectious individuals remains below ( λ − µ ) /σ , the population process will4e a super-critical process, having a positive probability to grow to ∞ . But, if the infec-tious fraction grows beyond ( λ − µ ) /σ , then the population process becomes a sub-criticalprocess. Remark 3.1. If ( λ − µ ) /σ ≥ i.e. λ ≥ µ + σ , then N(t) is always a super-critical processalthough 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 individ-uals is still small. E ( t ) increases with rate κI ( t ) S ( t ) /N ( t ) due to infection and decreasesby rate ( ν + µ ) E ( t ) due to death or end of latency. I ( t ) increases with rate νE ( t ) due toend of latency, and decreases with rate ( δ + µ + σ ) I ( t ) due to recovery or death. As westart with n − susceptible and one infected, assuming that n is large, at the beginningof the epidemic S ( t ) /N ( t ) is very close to . Then, the number of exposed individualsincreases almost with rate κ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 branchingprocess, L ∞ ( t ) = E ∞ ( t ) + I ∞ ( t ) with two stages: the childhood (or latent) stage E ∞ and the adult (or infectious) stage I ∞ [9, p. 54]. Initiated at ( E ∞ (0) , I ∞ (0)) = (1 , ,where E ∞ ( t ) increases with rate κI ∞ ( t ) and decreases with rate ( µ + ν ) E ∞ ( t ) , and I ∞ ( t ) increases with rate νE ∞ ( t ) (end of childhood) and decreases with rate ( δ + µ + σ ) I ∞ ( t ) (end of adult stage). Theorem 3.2.
Let L n ( t ) be the epidemic process and L ∞ ( t ) be the branching processdefined above. Then, L n ( t ) converges weakly to L ∞ ( t ) ( L n w = ⇒ L ∞ ), as n −→ ∞ , on anyfinite interval [0 , t ] .Proof. Like initial phase of other epidemics [9, p. 54], when n tends to infinity, thetransition probabilities of the epidemic converge to that of the branching process.The results below are for the branching process L ∞ , but since L n w = ⇒ L ∞ , they apply tothe epidemic as n −→ ∞ . In this subsection, we derive the Malthusian parameter α and the basic reproductionnumber R of the limiting branching process L ∞ .The Malthusian parameter α is defined as the exponential growth/decay rate the epidemicbranching process has. It is the solution of (cid:90) ∞ e − αt c ( t ) dt = 1 , (3.1)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 α = − (cid:18) µ + ν + δ + σ (cid:19) + (cid:114) ( ν − δ − σ ) κν. (3.2)5 roof. In the SEIRS model, the contact rate is during the latent period and κ duringthe infectious period. By conditioning on when the latent period ends, it follows that c ( t ) = κe − µt (cid:90) t νe − νs e − ( δ + σ )( t − s ) ds = κνe − ( µ + δ + σ ) t (cid:90) t e − ( ν − ( δ + σ )) s ds. Thus, c ( t ) = κνν − ( δ + σ ) (cid:0) e − ( µ + δ + σ ) t − e − ( ν + µ ) t (cid:1) , if ν (cid:54) = δ + σ,κνte − ( µ + δ + σ ) t , if ν = δ + σ. Inserting this into Equation (3.1), one gets α = − (cid:18) µ + ν + δ + σ (cid:19) + (cid:114) ( ν − δ − σ ) κν , if ν (cid:54) = δ + σ, √ κν − ( µ + δ + σ ) , if ν = δ + σ. Then, α = − (cid:18) µ + ν + δ + σ (cid:19) + (cid:114) ( ν − δ − σ ) κν. The basic reproduction number R , is the expected number of secondary cases per primarycase in a virgin population [9, page 4]. Theorem 3.4.
The basic reproduction number of the epidemic is R = νκ ( µ + ν )( δ + µ + σ ) . (3.3) Proof.
Let Y be the number of infectious contacts that an individual has during theinfectious period. Then, P ( Y = 0) = µµ + ν + νν + µ × δ + µ + σκ + δ + µ + σ . The first term is the probability that the individual dies during the latent period, andthe second term the probability that the individual does not die during the latent periodbut leaves the infectious compartment by death or recovery without an infectious contact.And for all positive integer k , P ( Y = k ) = νµ + ν × (cid:18) κκ + δ + µ + σ (cid:19) k × δ + µ + σκ + δ + µ + σ . Therefore, Z has a zero modified geometric distribution [13, page 16], with parameter p = ( δ + µ + σ ) / ( κ + δ + µ + σ ) . Then the expected value of Y is E ( Y ) = νµ + ν × − pp = νµ + ν × κµ + δ + σ . R = νµ + ν × κµ + δ + σ . Remark 3.5.
The first factor is the probability that the individual does not die duringthe latent period, and the second is the expected number of infectious contacts while beinginfectious. An alternative way to derive R is by the relation R = (cid:82) ∞ c ( t ) dt [13, page69]. This formula gives the same result as above. As mentioned above, the SI, SIS, SIR, SIRS, SEI, SEIS, SEIR models are all sub-models ofthe SEIRS model. Then, from the results obtained for the SEIRS model, we deduce thatof the others. In Table 2, we have the values of the basic reproduction number R and theMalthusian parameter α of the models listed above. The fourth column gives the changesto get the corresponding model. Further by setting σ = 0 , one gets the correspondingresults for a disease without an additional death rate for the infectious.Model R α Parameters change
SEIRS κν ( µ + ν )( µ + δ + σ ) − (cid:0) µ + ν + δ + σ (cid:1) + (cid:113) ( ν − δ − σ ) + κνSEIR κν ( µ + ν )( µ + δ + σ ) − (cid:0) µ + ν + δ + σ (cid:1) + (cid:113) ( ν − δ − σ ) + κν ρ = 0 SEIS κν ( µ + ν )( µ + δ + σ ) − (cid:0) µ + ν + δ + σ (cid:1) + (cid:113) ( ν − δ − σ ) + κν ρ → ∞ SEI κν ( µ + ν )( µ + σ ) − (cid:0) µ + ν + σ (cid:1) + (cid:113) ( ν − σ ) + κν δ = 0 , ρ = 0 SIRS κµ + δ + σ κ − ( µ + δ + σ ) ν → ∞ SIR κµ + δ + σ κ − ( µ + δ + σ ) ν → ∞ , ρ = 0 SIS κµ + δ + σ κ − ( µ + δ + σ ) ν → ∞ , ρ → ∞ SI κµ + σ κ − ( µ + σ ) ν → ∞ , δ = 0 , ρ = 0 Table 2: Thresholds for sub-models of the SEIRS model
Remark 3.6.
The basic reproduction number R and the Malthusian parameter α areidentical for the SEIR, SEIS and SEIRS models. From Equations (3.2) and (3.3), we have α > ⇐⇒ − (cid:18) µ + ν + δ + σ (cid:19) + (cid:114) ( ν − δ − σ ) κν > ⇐⇒ ( ν − δ − σ ) κν > (cid:18) µ + ν + δ + σ (cid:19) ⇐⇒ κν > (2 µ + ν + δ + σ ) − ( ν − δ − σ ) ⇐⇒ κν > (2 µ + 2 ν )(2 µ + 2 δ + 2 σ ) ⇐⇒ κν ( µ + ν )( µ + δ + σ ) > ⇐⇒ R > . It follows that the basic reproduction number R exceeds if and only if, the Malthusianparameter α exceeds . That is, the sign relation sign ( α ) = sign ( R − is verified.7 emark 3.7. It is well known that to surely prevent the disease to invade the population, R must be less than . To control the epidemic one need then to diminish R . The basicreproduction number R can be written in the following form. R = κ × νµ + ν × µ + δ + σ . So, it is clear that R increases with κ and ν , and decreases with µ, δ and σ . The contactrate κ can be reduced by hospitalization or quarantine of the infectious individuals. Therecovery rate δ can be increased by medication. In the case of an epizootic, σ the diseaserelated death rate, can be increased by culling infectious animals. We have derived two thresholds (the Malthusian parameter α and the basic reproduc-tion number R ) of the SEIRS epidemic branching process, deduced the correspond-ing thresholds for the sub-models of the SEIRS epidemic model and established that sign ( α ) = sign ( R − . In this subsection, we derive the probability for a minor outbreak of the epidemic branch-ing process, and state the main result for the initial phase of the epidemic.Let π = P (lim t L ∞ ( t ) = 0) be the probability of a minor outbreak of L ∞ and Y be the number of infectious contacts that an individual has during the infectious period. Aswe start with one latent individual, π is the smallest positive solution of the equation z = g ( z ) , where g is the probability generating function (pgf) of Y [13, page 113]. Wehave g ( z ) = ∞ (cid:88) k =0 P ( Y = k ) z k = µµ + ν + νµ + ν δ + µ + σκ + δ + µ + σ + ∞ (cid:88) k =1 νµ + ν δ + µ + σκ + δ + µ + σ (cid:18) κκ + δ + µ + σ (cid:19) k z k = a + (1 − a ) b − (1 − b ) z , with a = µµ + ν and b = δ + µ + σκ + δ + µ + σ . Then, π is the smallest solution in [0 , of the following equation. z = a + (1 − a ) b − (1 − b ) z . (3.4)Equation (3.4) has two solutions, z = 1 and z = a + b − b = µµ + ν + δ + µ + σκ = µµ + ν + νν + µ R . Then, we have the following theorem.
Theorem 3.8.
Let π be the probability of a minor outbreak of the epidemic when startedwith one latent individual. Then, π = if R ≤ ,µµ + ν + νµ + ν R if R > . (3.5)8 orollary 3.9. Let π ( m,k ) be the probability of a minor outbreak when the epidemic startedwith m latent individuals and k infectious individuals. Then, π ( m,k ) = if R ≤ , (cid:18) µµ + ν + νµ + ν R (cid:19) m (cid:18) R (cid:19) k if R > . (3.6) Proof.
Let π ( 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 π (1 , = P ( minor outbreak | E (0) = 1 , I (0) = 0) = π , and π (0 , = P ( minor outbreak | E (0) = 0 , I (0) = 1) . Thus, π = µ/ ( µ + ν ) + ( ν/ ( µ + ν )) π (0 , .Therefore, π (0 , = (( µ + ν ) /ν ) π − µ/ν . Thus, by Equation (3.5) one gets π (0 , = if R ≤ , R if R > . (3.7)We have π ( m,k ) = (cid:0) π (1 , (cid:1) m (cid:0) π (0 , (cid:1) k , since all the m + k independent epidemics must dieout [13, page 112]. Thus, by Equations (3.5) and (3.7), we get Equation (3.6). Remark 3.10.
This result is the same as that of the SEIR stochastic model studied byAllen and Lahodny in [1].
As noted above, due to the additive death rate in the infectious compartment, if theepidemic takes off, the population process can be turned to a sub-critical process. Inthis 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 √ n , that is P ( L n ( t ) < √ n, ∀ t ≥ [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 − , , , , 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 ∞ ( t ) denote the birth and death process defined above, α its Malthusian parameter, π the probability of a minor outbreak of L ∞ , and π n := P ( L n ( t ) < √ n, ∀ t ≥ denote the probability of a minor outbreak of the epidemic. Thenas n → ∞ , we have the following results: i. If α ≤ , then for any n, L n ( t ) → as t → ∞ with probability . ii. If < α < λ − µ , then π n → π = µ/ ( µ + ν ) + ν/ [( µ + ν ) R ] . With the remainingprobability (1 − π n ) → − π , L n grows exponentially: L n ( t ) ∼ e αt , but L n ( t ) /N n ( t ) → as t → ∞ . iii. If α > λ − µ , then π n → π = µ/ ( µ + ν ) + ν/ [( µ + ν ) R ] . With the remainingprobability (1 − π n ) → − π , during the initial phase of the epidemic, L n growsexponentially with rate α .Proof. i. If α < , then L ∞ is sub-critical and dies out with probability . If α = 0 ,then L ∞ is critical and dies out with probability , since P ( Y = 1) (cid:54) = 1 [13].9i. If α > , then L ∞ is super-critical. It dies out with probability π and with theremaining probability − π , it grows exponentially with rate α , that is L ∞ ∼ e αt .As N n ( t ) ∼ e ( λ − µ ) t , if α < λ − µ , then L n ( t ) /N n ( t ) −→ , when t −→ ∞ .iii. If α > λ − µ , then α > since λ − µ > . Thus, L ∞ is super-critical. It dies outwith probability π and with the remaining probability − π , it grows exponentiallywith rate α . As L n = ⇒ L ∞ , the same applies to L n . Remark 3.12. If α > λ − µ , then L ∞ is super-critical. Thus the epidemic may takeoff. If it does, the number of infected individuals grows initially with the rate α that islarger than the initial growth rate of the population ( λ − µ ). After the initial phase, inthe case of a major outbreak, several scenarios are possible. i)The population goes ongrowing exponentially, eventually with a lower rate; ii) due to the additive death rate ofthe 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 aretreated 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 studiedabove. As the population size is varying, we study first the fractions system and then de-duce the asymptotic behavior of the compartments sizes. In this section the deterministicsizes 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 followingsystem of ordinary differential equations (ODE). dSdt = λN + ρR − κS IN − µS,dEdt = κS IN − ( ν + µ ) E,dIdt = νE − ( δ + µ + σ ) I,dRdt = δI − ( ρ + µ ) R,N = S + E + I + R,S (0) > , E (0) > , I (0) > , R (0) ≥ . (4.1) Remark 4.1.
System (4.1) is the same as System (2) studied by Greenhalgh in [11],with a constant contact rate ( β ( N ) = β ), a constant death rate ( f ( N ) = µ ) and withoutvaccination ( p = q = 0 ). But Greenhalgh assumed that the death rate f ( N ) is a strictlymonotone increasing continuously differentiable function of N. So, the model that we studyis not a sub-model of that of Greenhalgh since we consider a constant death rate. dN/dt = ( λ − µ ) N − σI = ( λ − µ − σi ) N , where i is the fraction of the infectious individuals. Thus, the population should grow if λ > µ + σi ,stabilize if λ = µ + σi and decrease if λ < µ + σi . Theorem 4.2. N ( t ) is constant and positive ( N ( t ) = N (0) > , ∀ t > ), if and only ifthe parameters verify the following equality: ( ρ + µ ) νκσλ + ρκνδ ( λ − µ ) − ( ρ + µ )( ν + µ )( δ + µ + σ )[ κ ( λ − µ ) + µσ ] = 0 , (4.2) and the initial values verify S (0) = ( νκ ) − ( ν + µ )( δ + µ + σ ) N (0) ,E (0) = ( νσ ) − ( δ + µ + σ )( λ − µ ) N (0) ,I (0) = σ − ( λ − µ ) N (0) ,R (0) = ( σ ( ρ + µ )) − δ ( λ − µ ) N (0) , with N (0) > . (4.3) Proof.
By using successively the derivatives of
N, I, E, R and S , one gets that N ( t ) isconstant and positive, if and only if I, E, R and S are constant, the parameters verifyEquation (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 . BySystem (4.1), one gets dsdt = λ − λs + ρr + ( σ − κ ) si,dedt = κsi − ( λ + ν ) e + σei,didt = νe − ( λ + δ + σ ) i + σi ,drdt = − ( λ + ρ ) r + δi + σri,s + e + i + r = 1 . (4.4)This system is equivalent to System (2.2) in [11], with p = q = 0 . Remark 4.3.
The natural death rate µ does not intervene in the derivatives of the frac-tions. This is logic, since this rate is the same for all the compartments, it has no effecton the fractions. r = 1 − s − e − i , it is enough to study the system dsdt = λ + ρ − ( λ + ρ ) s − ρe − ρi + ( σ − κ ) si,dedt = κsi − ( λ + ν ) e + σei,didt = νe − ( λ + δ + σ ) i + σi , (4.5)in the domain D = { ( s, e, i ); s ≥ , e ≥ , i ≥ , s + e + i ≤ } . (4.6) Theorem 4.4.
The domain D is positively invariant for System 4.5.Proof. If s = 0 , then ds/dt = λ + ρ (1 − e − i ) > . If e = 0 , then de/dt = κsi ≥ . If i = 0 , then di/dt = νe ≥ . If s + e + i = 1 , then d ( s + e + i ) /dt = − δi ≤ . Then everysolution of System 4.5 starting in D , remains there for all t > . That is D is positivelyinvariant for System 4.5In the following, we use the next generation matrix (NGM) described in [10] to derive athreshold parameter, with threshold value 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 , , is the unique disease free equilibrium (DFE) of System (4.5). e and i arethe infected fractions of the model. Thus, the infected subsystem is dedt = κsi − ( λ + ν ) e + σei,didt = νe − ( λ + δ + σ ) i + σi . Hence, the linearized infected subsystem at the DFE is dedt = κi − ( λ + ν ) e,didt = νe − ( λ + δ + σ ) i. (4.7)Let x = ( e, i ) t be the vector of the infected fractions . Thus, System (4.7) is equivalent to . x = ( T + Σ) x, with T = (cid:18) κ (cid:19) and Σ = (cid:18) − ( λ + ν ) 0 ν − ( λ + δ + σ ) (cid:19) .T is the transmissions matrix and Σ is the transitions matrix. Then, the next generationmatrix with large domain is K L = − T Σ − . Let R be the spectral radius of K L . We have R = ρ ( K L ) = κν ( λ + ν )( λ + δ + σ ) . (4.8)12n 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 of x ∗ , such that every solution starting in V remainsin V for all t > . If V can be chosen such that lim t →∞ x ( t ) = x ∗ , then x ∗ is said to beasymptotically stable. An equilibrium is said to be unstable, when it is not stable. Anequilibrium x ∗ is said to be globally asymptotically stable (GAS) in an invariant set D , ( x ∗ ∈ D ) , if it is locally stable and lim t →∞ 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 asymptot-ically stable (GAS) in D , if R ≤ , and unstable if R > .Proof. Let f ( s, e, i ) be the Rhs of System (4.5). Then, f ( s, e, i ) = λ + ρ − ( λ + ρ ) s − ρe − ρi + ( σ − κ ) siκsi − ( λ + ν ) e + σeiνe − ( λ + δ + σ ) i + σi . The Jacobian of f at the disease free equilibrium is Df (1 , ,
0) = − ( λ + ρ ) − ρ σ − κ − ρ − ( λ + ν ) κ ν − ( λ + δ + σ ) . The characteristic polynomial of Df (1 , , is P ( x ) = ( − λ − ρ − x )[ x + (2 λ + ν + δ + σ ) x + ( λ + ν )( λ + δ + σ ) − νκ ] . − ( λ + ρ ) 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 ( λ + ν )( λ + δ + σ ) − νκ > .And we have ( λ + ν )( λ + δ + σ ) − νκ > ⇐⇒ R < . Thus, the disease free equilibriumis locally asymptotically stable if R < , and unstable if R > .Let V denote the function defined on D by V ( s, e, i ) = νe + ( λ + ν ) i . Then, . V ( s, e, i ) = ν dedt + ( λ + ν ) didt = ν [ κsi − ( λ + ν ) e + σei ] + ( λ + ν )[ νe − ( λ + δ + σ ) i + σi ]= i [ νκs + νσe + ( λ + ν ) σi − ( λ + ν )( λ + δ + σ )]= iL ( s, e, i ) , with L ( s, e, i ) = νκs + νσe + ( λ + ν ) σi − ( λ + ν )( λ + δ + σ ) . The affinity of L implies thatit achieves its maximum at the extreme points of the boundary of the closed set D . But L (0 , ,
0) = − ( λ + ν )( λ + δ + σ ) , L (0 , ,
1) = − ( λ + ν )( λ + δ ) , L (0 , ,
0) = − λσ − ( λ + ν )( λ + δ ) and L (1 , ,
0) = νκ − ( λ + ν )( λ + δ + σ ) = ( λ + ν )( λ + δ + σ )( R − . Thus, . V ≤ in D if R ≤ . Then, V is a Lyapunov function of System (4.5). The only invariant subsetof the set with . V = 0 is { (1 , , } . 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 ≤ . 13y Equations (4.8) and (3.2), we have R > ⇐⇒ κν ( λ + ν )( λ + δ + σ ) > ⇐⇒ κν > ( λ + ν )( λ + δ + σ ) ⇐⇒ κν > (cid:0) (2 λ + ν + δ + σ ) − ( ν − δ − σ ) (cid:1) ⇐⇒ ( ν − δ − σ ) κν > (2 λ + ν + δ + σ ) ⇐⇒ (cid:114) ( ν − δ − σ ) κν > λ + ν + δ + σ ⇐⇒ − (cid:18) µ + ν + δ + σ (cid:19) + (cid:114) ( ν − δ − σ ) κν > λ − µ ⇐⇒ α > λ − µ. It follows that the fraction’s threshold R exceeds if and only if the Malthusian parameter α , exceeds the initial growth rate of the population λ − µ . That is, we have the sign relation sign ( R −
1) = sign ( α − ( λ − µ )) . Thus, the global stability of the disease free equilibriumof the fraction’s system when R < , confirms that if α < λ − µ , 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 > , then System (4.5) has at least oneendemic equilibrium, and that this equilibrium is unique and locally asymptotically stable(LAS) when the average duration of immunity exceeds both the average infectious andincubation periods, that is δ > ρ and ν > ρ . We have not proved, but we strongly believethat if R > , then System (4.5) has one and only one endemic equilibrium, and thatthis equilibrium is globally asymptotically stable in the interior of D . The simulationsthat 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 denotedthe basic reproduction number given by Equation (3.3). i. If R < , then ( S ( t ) , E ( t ) , I ( t ) , R ( t )) −→ ( ∞ , , , ; ii. if R = 1 , then ( S ( t ) , E ( t ) , I ( t ) , R ( t )) −→ ( ∞ , E ∗ , I ∗ , R ∗ ) ,with E ∗ > , I ∗ > , R ∗ > ; iii. if R > ≥ R , then ( S ( t ) , E ( t ) , I ( t ) , R ( t )) −→ ( ∞ , ∞ , ∞ , ∞ ) . Remark 4.7.
The case R > is treated in Theorem 4.9.Proof. We have R = νκ ( µ + ν )( µ + δ + σ ) and R = νκ ( λ + ν )( λ + δ + σ ) . Then, R > R , since λ > µ . Therefore, in the three cases of the Theorem 4.6, we have R ≤ . Let us assume that R ≤ . Thus, by Theorem (4.5), ( s, e, i, r ) −→ (1 , , , when t −→ ∞ . dN/dt = ( λ − µ ) N − σI = ( λ − µ − σi ) N . Then, dN/dt −→ ( λ − µ ) N ,14hen t −→ ∞ . Thus, N −→ ∞ , when t −→ ∞ because λ > µ . Therefore, S −→ ∞ ,when t −→ ∞ , since S/N −→ . By using the derivatives of E and I , one gets (cid:18) EI (cid:19) (cid:48) = κs + ( δ + σ − ν ) EI − ν (cid:18) EI (cid:19) . Where the prime denotes the derivative . Then, (cid:18) EI (cid:19) (cid:48) −→ κ + ( δ + σ − ν ) EI − ν (cid:18) EI (cid:19) , when t −→ ∞ . Thus
E/I can be approximate by a solution of the following equation, when t −→ ∞ . y (cid:48) = κ + ( δ + σ − ν ) y − νy (4.9)Equation (4.9) is a Riccati’s equation [12]. By solving it, one gets y : t (cid:55)−→ (cid:18) Ce √ ∆ t − ν √ ∆ (cid:19) − + δ + σ − ν + √ ∆2 ν , with C > , where ∆ = ( δ + σ − ν ) +4 νκ. Then,
E/I −→ ( δ + σ − ν + √ ∆) / (2 ν ) , when t −→ ∞ .We have dI/dt = νE − ( δ + µ + σ ) I = [ ν ( E/I ) − ( δ + µ + σ )] I . Thus, by substituting E/I by its asymptotic value, one gets dIdt −→ (cid:34) δ + σ − ν + √ ∆2 − ( δ + µ + σ ) (cid:35) I , when t −→ ∞ ; and δ + σ − ν + √ ∆2 − ( δ + µ + σ ) = − (cid:18) µ + δ + σ + ν (cid:19) + √ ∆2= − (cid:18) µ + δ + σ + ν (cid:19) + (cid:114) ( δ + σ − ν ) κν = α , the Malthusian parameter given by Equation (3.2) . Therefore, dI/dt −→ αI , when t −→ ∞ .As dE/dt = κSI/N − ( ν + µ ) E , S/N −→ and E/I −→ ( δ + σ − ν + √ ∆) / (2 ν ) , aftersome algebra, one gets dE/dt −→ αE , when t −→ ∞ . As sign ( α ) = sign ( R − , ( E, I ) −→ (0 , if R < E, I ) −→ ( E ∗ , I ∗ ) , with E ∗ > and I ∗ > , if R = 1;( E, I ) −→ ( ∞ , ∞ ) if R > . For the number of the recovered R , as dR/dt = δI − ( ρ + µ ) R , it is obvious that R hasthe same asymptotic behavior as I . Remark 4.8.
In the proof, we have shown that if R ≤ , then the Malthusian parameter α of the stochastic model, is also the common asymptotic growth rate of the compartmentsE and I, and λ − µ is the asymptotic growth rate of the population. Since sign ( R −
1) = sign ( α − ( λ − µ )) , this is coherent with Theorem 3.11 (ii). R ≤ . If R > , then the fraction disease free equilibrium is unstable. Therefore, the disease willremain endemic in the population in term of the fraction infected. The following theoremgives the asymptotic behavior of the compartments sizes, when R > , assuming thatthe fraction system admits an endemic equilibrium that is globally asymptotically stablein the interior of the feasible region D . Theorem 4.9.
Assume that R > and that System (4.5) has a unique endemic equilib-rium ( s ∗ , e ∗ , i ∗ ) that is globally asymptotically stable in o D , and set R = λµ + σi ∗ . (4.10)i. If R > , then ( S, E, I, R ) −→ ( ∞ , ∞ , ∞ , ∞ ) . ii. If R = 1 , then ( S, E, I, R ) −→ ( S ∗ , E ∗ , I ∗ , R ∗ ) ,with S ∗ > , E ∗ > , I ∗ > , R ∗ > . iii. If R < , then ( S, E, I, R ) −→ (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 o D . Then, dNdt −→ ( λ − µ − σi ∗ ) N , when t −→ ∞ . Asymptotically, the population would increase with rate λ , and decrease with rate µ + σi ∗ .As the fraction system admits an endemic equilibrium, that is globally asymptoticallystable in the interior of the feasible region, all the compartments have the same asymptoticbehavior as the population. Let us set α = λ − µ − σi ∗ . The quantity α is the commonasymptotic exponential growth/decay rate of all the compartments S, E, I and R. Wehave sign ( α ) = sign ( R − . Thus, the results follow.In this section we have studied the corresponding deterministic SEIRS model of the pre-vious stochastic model. We derived a threshold quantity R for the fraction model. If R ≤ , then the fraction ’s disease free equilibrium is globally asymptotically stable inthe feasible region D , otherwise it is unstable. When R ≤ , the behavior of the number of infected is determined by the basic reproduction number R . If R < , then the num-ber of infected vanishes. If R = 1 , then the number of infected stabilizes to a positivevalue; when R ≤ < R , then the number of infected grows exponentially, but at alower rate than the population. If R > , then the number of infected grows initiallyquicker than the population and the asymptotic behavior of the population is governedby the threshold quantity R . If R < , then the population vanishes; if R = 1 , then thepopulation stabilizes; if R > , then the population grows, but with a lower rate thanits initial growth rate.
5. Simulations
In this section, we use the software R to illustrate and confirm the results found in theprevious sections. In the following, we set µ = 1 , that is the time unit is the life expectancy,16xcept for the simulations of influenza epidemics in Burkina Faso where we set one yearas the time unit. The other parameters and the initial values are chosen arbitrary, unlessotherwise stated. In this subsection, we give some examples of simulations of epidemics starting by one latentindividual, and using different values of the parameters. The population is initiated with latent individual and susceptible individuals.In Figure 2 (a) and (b), where R = 0 . and R = 0 . respectively, all the epidemicsdie without any major outbreak. In (a) the maximum of infected individuals is , whileit is in (b). (a) R = 0 . (b) R = 0 . Figure 2: (a) 10 SEIRS simulations with λ = 3 , µ = 1 , σ = 7 , δ = 15 , κ = 10 , ν = 20 , ρ = 5 , thatgives R = 0 . , α = − . , π = 1 . (b) 10 SEIRS simulations with λ = 1 . , µ = 1 , σ = 7 , δ = 5 , κ =10 , ν = 20 , ρ = 5 , that gives R = 0 . , α = − . , π = 1 . In both cases, all the epidemics die outwithout a major outbreak. However, they die out quicker and the number of infected is fewer in (a) thanin (b). In Figure 3, where R = 2 , four simulated epidemics out of have a major outbreak.The other epidemics die out without many getting infected. For the epidemics withmajor outbreak, the number of the infected individuals grow exponentially but the timewhere the exponential growth starts varies.Now we estimate the probability of a minor outbreak π by simulating epidemics andsetting ˆ π n = n / , where n is the number of minor epidemics. We set λ = 3 , µ =1 , ν = 50 , δ = 10 , σ = 4 , ρ = 3 , and set successively κ = 10 , , , , to get differentvalues of π . Table 3 gives the different values of π and the estimate ˆ π n for n = 1000 and n = 2000 respectively. These results confirm that the probability of extinction of thebranching process L ∞ is a good approximation of the probability of a minor outbreak ofthe epidemic starting with one latent individual in a population of size n , when n is large.17 a) (b)Figure 3: 10 SEIRS simulations ( dying quickly) with λ = 3 , µ = 1 , σ = 4 , δ = 5 , κ = 21 , ν = 20 , ρ =5 , n = 1000 , R = 2 , R = 1 . , α = 5 . , π = 0 . . In (b) we made a zoom so that the six minorepidemics can be seen. n κ
10 10 20 20 30 30 50 50 100 100 R .
654 0 .
654 1 .
307 1 .
307 1 .
961 1 .
961 3 .
268 3 .
268 6 .
536 6 . π .
000 0 .
770 0 .
770 0 .
520 0 .
520 0 .
320 0 .
320 0 .
170 0 . π n .
000 1 .
000 0 .
788 0 .
782 0 .
513 0 .
520 0 .
303 0 .
329 0 .
165 0 . Table 3: Estimation of the probability of minor outbreak π n of the epidemic starting with one latentindividual and n − susceptible individuals. The theoretical result is π and the simulated result is ˆ π n from simulations. These simulations confirm that when R is less than , the disease cannot invade thepopulation, and that if R is larger than , then with a positive probability (1 − π n ) −→ (1 − π ) , the disease can invade the population, and in this case the number of the infectedindividuals grows exponentially during the initial phase. In this subsection, we show some simulations of major outbreaks, where the epidemicstarts with a positive fraction of infected individuals. So the simulations illustrate whathappens once the number of infected has reached a small but positive fraction of thecommunity. We use the blue color for the susceptible, green for the exposed, red for theinfectious and black for the recovered.We start by simulating the deterministic fraction ’s system. In Figure 4, we have fourcases with R = 0 . , , . , . respectively. In each case we have ten solutions pathsof System (4.4) with different initial values. In (a) as in (b), the solutions of System(4.4) approach the disease free equilibrium, confirming that if R ≤ , then the diseasefree equilibrium of the deterministic fraction system is globally asymptotically stable inthe feasible region. In (c), all the solutions approach the same endemic equilibrium18 s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . . In (d), all the solutions approach the sameendemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . . The results of (c) and (d)confirm that when R > , the disease free equilibrium is unstable and that there is anendemic equilibrium that is globally asymptotically stable in the interior of the feasibleregion. That is the endemic level is independent of the starting value. However theendemic equilibrium of (d) is different of that of (c), hence the endemic equilibrium varywith the parameters values. (a) R = 0 . (b) R = 1 (c) R = 1 . (d) R = 5 . Figure 4: Simulations of of System (4.4). In each case we have 10 solutions paths of System (4.4)with different initial values. For (a) we have λ = 2 , µ = 1 , σ = 8 , δ = 5 , κ = 9 , ν = 10 , ρ = 10 thatgives R = 0 . . For (b) the parameters have the same values as in (a) except that we set κ = 18 to get R = 1 . In (a) and in (b) all the solutions approach the disease free equilibrium (1 , , , .The time scale is longer in (b), so the epidemic takes longer time to die out when R is close to .For (c) we have λ = 2 , µ = 1 , σ = 8 , δ = 5 , κ = 30 , ν = 15 , ρ = 2 , that gives R = 1 . ; all thesolutions approach the same endemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . . For (d), wehave λ = 2 , µ = 1 , σ = 8 , δ = 5 , κ = 100 , ν = 8 , ρ = 20 that gives R = 5 . ; all the solutions approachthe same endemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . .
19n the following we simulate both the stochastic and the deterministic models. In eachcase, we have simulated the stochastic epidemic, as well as integrated numerically thedeterministic system (4.1), both starting at the same values. From the numbers , we gotthe fractions by setting s = S/N, e = E/N, i = I/N, r = R/N, with N = S + E + I + R . One distinguishes the stochastic solutions from the deterministic by the fact thatthe deterministic solutions are represented by smooth lines, while the solutions of thestochastic solutions are represented by broken lines.In Figure 5, where R = 1 . > and R = 0 . < , the numbers of latent, infectiousand recovered grow exponentially as the population, but the population growth rate iseven 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 diesout in term of the fractions . (a) numbers (b) Fractions
Figure 5: SEIRS curves with λ = 3 , µ = 1 , σ = 3 , δ = 5 , κ = 12 , ν = 20 , ρ = 3 that gives R =1 . , R = 0 . , π = 0 . , α = 1 . The initial values are N (0) = 1000 with ( S (0) , E (0) , I (0) , R (0)) =(800 , , , . All the numbers grow exponentially. But the population grow faster, and the fractionsgo to the disease free equilibrium. One distinguishes hardly the deterministic curves from the stochasticbecause they are very close. In Figure 6, initially the number of the infected grows quicker than the number of thesusceptible. After that, the number of the susceptible declines. Afterward all the com-partments grow exponentially, but with a rate α ≈ . that is lower than the initialgrowth rate of the population λ − µ = 1 . The fractions go to an endemic equilibrium.For Figure 7, the parameters are chosen such that Equation (4.2) is verified, allowingthen the existence of endemic equilibrium for the deterministic System 4.1. The asymp-totic reproduction number of the population R = 1 , then for the deterministic solution,the population stabilizes, when t −→ ∞ . The stochastic numbers fluctuate around thedeterministic numbers . The fractions have the same behavior as that of the numbers .In Figure 8, for the deterministic model, the epidemic turned the population exponentialgrowth to an exponential decay due to the disease induced death rate, while the frac-tions go to an endemic equilibrium. For the stochastic model, the population size first20 a) Numbers (b) Fractions
Figure 6: SEIRS curves with λ = 2 , µ = 1 , σ = 4 , δ = 5 , κ = 21 , ν = 20 , ρ = 5 that gives R =2 , R = 1 . , π = 0 . , α = 5 . . The initial values are ( S (0) , E (0) , I (0) , R (0)) = (980 , , , . Allthe numbers grow exponentially with rate α ≈ . , while the fractions go to an endemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . .(a) Numbers (b)
Fractions
Figure 7: SEIRS curves with λ = 1 . , µ = 1 , σ = 5 , δ = 6 , κ ≈ . , ν = 20 , ρ = 3 that gives R = 1 . , R = 1 . , α = 5 . , The initial values are ( S (0) , E (0) , I (0) , R (0)) = (900 , , , . Theasymptotic growth rate of the population is α = 0 , its asymptotic reproduction number R = 1 .Thus for the deterministic solutions the numbers approach an endemic equilibrium ( S ∗ , E ∗ , I ∗ , R ∗ ) ≈ (581 , , , and the fractions go to an endemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . .The stochastic solutions fluctuate around the deterministic. decreases but at some instance when there are only few remaining individuals, the diseasegoes extinct and then the population size starts growing again. The deterministic modelsuggests that the population will go extinct, whereas in the stochastic model the diseasefirst dies out and then the population becomes super critical again, thus regrowing. In21he stochastic setting, what happens when the numbers become low is random. Boththe disease and the population could die out, or the population can grow again afterthe extinction of the epidemic. Only analyzing the deterministic fraction model wouldgive a misleading conclusion since the fractions seem to stabilize, whereas in what reallyhappens is that all numbers in the deterministic model tend to . (a) Numbers (b)
Fractions
Figure 8: SEIRS curves with λ = 2 , µ = 1 , σ = 8 , δ = 5 , κ = 30 , ν = 20 , ρ = 3 that gives R =2 . , R = 1 . , π = 0 . , α = 7 . , the initial values are ( S (0) , E (0) , I (0) , R (0)) = (980 , , , . Theasymptotic decay rate of the population is α ≈ − . , its asymptotic reproduction number R ≈ . . R > and R < thus, in the deterministic solutions, the numbers vanish, while the fractions go to anendemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . . In the stochastic solution, all the numbersvanish except the number of the susceptible that decreases until the extinction of the disease and regrowexponentially after that. In Figure 9, we have the case where R = 1 . For the stochastic model the disease goesextinct, while it persists in the deterministic one. In both models the population goes ongrowing 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 thestochastic model and the deterministic model. Now we simulate two influenza epidemics in Burkina Faso with two different basic repro-duction numbers. Burkina Faso is an inland country of West Africa. Its population in is about
16 000 000 . The individual annual birth rate and death rate are estimatedto λ = 0 . and µ = 0 . respectively [15].According to the World Health Organization (WHO) [21], influenza is caused by a virusthat attacks mainly the respiratory tract: the nose, the throat, the bronchi and rarely alsothe lungs. The infection usually lasts for about a week. It is characterized by sudden onsetof 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 medical22 a) Numbers R = 1 (b) Numbers R = 1 (Zoom-in) Figure 9: SEIRS curves with λ = 2 , µ = 1 , σ = 8 , δ = 9 , κ = 19 . , ν = 15 , ρ = 3 that gives R =1 , the initial values are ( S (0) , E (0) , I (0) , R (0)) = (90 , , , . The population grows exponentially inthe stochastic and in the deterministic model. In the stochastic model the latent, the infectious andthe recovered vanish, while in the deterministic model, they stabilize to positive values ( E ∗ , I ∗ , R ∗ ) ≈ (1 . , . , . . In (b) the scale is chosen to show the dynamics of E, I and R ; S ( t ) is much larger andnot shown. treatment. The virus is easily passed from person to person through the air by dropletsand small particles excreted when infected individuals cough or sneeze. The influenzavirus enters the body through the nose or the throat. It then takes between one andfour days for the person to develop symptoms. Someone suffering from influenza can beinfectious from the day before he/she develops symptoms until seven days afterwards. Thedisease 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 otherconditions 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-roundtend to have high attack and case-fatality rates. Therefore by setting one year as thetime unit, we make the following estimation for the influenza parameters. The latentperiod is approximately . days, the infectious period is approximately days and theimmunity period is approximately one year, thus ν = 365 / . , δ = 365 / , ρ = 1 . Weset the influenza case fatality rate (CFR) to . and deduce the influenza related deathrate σ ≈ . .The reproduction number of the pandemic influenza is estimated to be between and [19]. Thus we set R = 2 . and deduce the contact number κ = 130 from the otherparameters and Equation (3.3). We simulate the epidemic for a period of years startingin by integrating numerically the deterministic System (4.1). Figure 10 gives theevolution of the numbers of susceptible, latent, infectious and recovered individuals andthe corresponding fractions during the years. Figures 10 (c) and (d), show respectivelythe dynamics of E and I and that of the fractions e and i . We have a peak with infectious individuals at the th week of the epidemic. After that the number of infec-23ious declines because the number of susceptible is low. Afterward we have a minor peakevery year due to the immunity waning and the newborns that increase the number ofsusceptible. The number of recovered individuals grow quickly and reach its maximum
13 139 592 at the th week. The fractions go to an endemic equilibrium through dampedoscillations. The population in is estimated to
22 370 000 individuals with susceptible,
94 000 latent individuals,
260 000 infectious individuals and
13 060 000 recov-ered with non-permanent immunity. The number of recovered individuals is larger thanthe number of susceptible individuals. The fractions go to an endemic equilibrium withmore than of the population infected at every time. As the influenza last about oneweek and we have weeks within a year, more than of the population should beinfected during the year . That will have a very important negative impact on theeconomy of the country. We have R ≈ . and R ≈ . , thus R and R are bothlarger than and according to Theorem 4.9 the fractions should go to an endemic equi-librium and all the compartments should grow exponentially. Therefore the simulationsagree with the theoretical results.The basic reproduction number for the novel influenza A (H1N1) has been estimated tobe between . and . [8]. Thus, we set now R = 1 . and deduce the contact numberfrom the other parameters. The results of this simulation are shown in Figure 11. In thiscase the epidemic and the population have globally the same dynamics as above. Butthe impact of the epidemic is fewer. The major peak of the epidemic happens later, atthe th week, with
794 004 infectious individuals. Contrary to the preceding case, thenumber 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 modelwith 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 simulationsof influenza in Burkina Faso show that in spite of the epidemic the population should goon growing. But the number of infected individuals will grow also. Furthermore the peakof the epidemic in the first year show that the emergence of a new strain of influenzavirus will be a very serious threat for the world. The Global Influenza Program (GIP) ofthe World Health Organization (WHO), provides Member States with strategic guidance,technical support and coordination of activities essential to make their health systemsbetter prepared against this threat.In this section we have illustrated and validated the theoretical results of the previoussections by simulations.
6. Conclusion and discussions
In this paper we have studied a stochastic SEIRS epidemic model, with a disease relateddeath, in a population which grows exponentially without the disease. We assumed thatinitially, the population process is a super-critical linear birth and death process withbirth rate λ and death rate µ .We have derived, the basic reproduction number R , the Malthusian parameter α and theprobability of minor outbreak π assuming that the initial size n of the population tendsto infinity. Considering the deterministic model, we derived the threshold quantity R forthe fractions . 24 a) Numbers (b)
Fractions (c)
Numbers of infected (d)
Fractions of infected
Figure 10: Simulation of influenza in Burkina Faso with R = 2 . . The parameters values are λ = 0 . , µ = 0 . , σ = 0 . , δ = 365 / , κ = 130 . , ν = 365 / . , ρ = 1 that gives R = 2 . , R = 2 . . The initial values are ( S (0) , E (0) , I (0) , R (0)) = (16 000 000 , , , . By dampedoscillations, the fractions approach an endemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . .The initial growth rate of the population is λ − µ = 0 . , with the epidemic its asymptotic growth rateis α ≈ . and its asymptotic reproduction number rate is R ≈ . . (c) show the dynamics of E and I ; (d) show the dynamics of e and i . If R ≤ , then the disease dies out with probability . That is, there is no possibility ofmajor outbreak if R ≤ . In this case, the population remains a super-critical process.If R > then, with a positive probability, the epidemic can take off. If the epidemic takesoff, then the number of the infected (exposed or infectious) individuals grows exponentiallywith rate α . If < α ≤ λ − µ , then the sizes of all the compartments grow exponentiallywhile the fraction of the infected individuals goes to . The number of infected peoplegrows, but at a lower rate than the population, implying that the fraction infected becomesnegligible. If α > λ − µ , then the number of infected individuals grows initially with a25 a) Numbers (b)
Fractions (c)
Numbers (d)
Fractions
Figure 11: Simulation of influenza in Burkina Faso with R = 1 . . The parameters values are λ =0 . , µ = 0 . , σ = 0 . , δ = 365 / , κ = 78 . , ν = 365 / . , ρ = 1 that gives R = 1 . , R = 1 . .The initial values are ( S (0) , E (0) , I (0) , R (0)) = (16 000 000 , , , . By damped oscillations, thefractions approach an endemic equilibrium ( s ∗ , e ∗ , i ∗ , r ∗ ) ≈ (0 . , . , . , . . The initial growthrate of the population is λ − µ = 0 . , with the epidemic its asymptotic growth rate is α ≈ . and its asymptotic reproduction number is R ≈ . . (c) show the dynamics of E and I ; (d) show thedynamics of e and i . rate that is larger than the population growth rate, and different scenarios are possible.Due to the additional death rate σ in the infectious compartment, the population will goon growing but with a lower rate, or the population will become a sub-critical branchingprocess and thus have a decreasing size. In the latter case the population vanishes in thedeterministic model while in the stochastic one, what happens when the numbers becomelow is random. Both the disease and the population could die out, or the populationcould grow again after the extinction of the epidemic.We have illustrated and validated the theoretical results by simulations. These simulations26how the similarities and the differences between the stochastic model and the correspond-ing deterministic model. If R > , then in the deterministic model, the epidemic willinvade the population surely; while in the stochastic model, with a positive probability π the disease vanishes. When R < , the population vanishes surely in the deterministicmodel; whereas in the stochastic model, it can regrow exponentially after the extinction ofthe epidemic. One need to remember that an epidemic is always a stochastic process, andthe deterministic model fits only when we have a large community with a large numberof infectious individuals [6].For some diseases (e.g. influenza, measles), the susceptibility and the infectiousness varywith the age of the individuals. Therefore it is more realistic to consider an heteroge-neous population as in [18]. Furthermore adding seasonal forcing will increase realismfor seasonal diseases. For other diseases like sexually transmitted diseases (STD), a closeand/or long contact is required for transmission. Then the dynamics of the epidemic arelinked to the special network in the host population [4, 2, 18]. Due to the developmentof the migration of populations, an epidemic starting in one location can be exported toanother one quickly, then a meta-population model is convenient to find the conditionsfor a global health security [3, 7]. Under the leadership of the World Health Organiza-tion (WHO) many policies (vaccination, medication, quarantine, etc.) are implementedto prevent major outbreak of epidemics. Nevertheless the infectious diseases remain aserious threat for Humanity. Another step towards realism is hence to extend this modelby adding vaccination [6, 11] and treatment [17]. This should allow to find the optimalcontrol strategy to realize the herd immunity, that is to prevent the major outbreaks ofinfectious diseases. Acknowledgments
We are grateful to the International Science Program (ISP) of Uppsala University andthe Swedish International Development Agency (SIDA) for their financial support.
References [1] L. J. Allen and GE. Jr. Lahodny. ( ). Extinction Thresholds in Deterministicand Stochastic Epidemic Models.
Journal of Biological Dynamics , , − .https://doi.org/10.1080/17513758.2012.665502[2] F. Ball, T. Britton and D. Sirl. ( ). A Network with Tunable Clustering, De-gree Correlation and Degree Distribution, and an Epidemic Thereon. Journal ofMathematical Biology , Vol. , Issue , pp − . http://dx.doi:10.1007/s00285-012-0609-7[3] F. Ball, T. Britton, T. House, V. Isham, D. Mollisone, L. Pellis and G. S. Tomba.( ). Seven Challenges for Metapopulation Models of Epidemics, Including House-holds Models. Epidemics
10 63 − . http://dx.doi.org/10.1016/j.epidem.2014.08.001[4] F. Brauer. ( ). An Introduction to Networks in Epidemic Modeling, in Mathemat-ical Epidemiology,
Lecture Notes in Mathematics. Vol 1945 pp − . Springer-Verlag, Berlin. 275] T. Britton and P. Trapman. ( ). Stochastic Epidemics in Growing Populations. Bull. Math. Biol.
76: 985-996. http://dx.doi:10.1007/s11538-014-9942-x[6] T. Britton. ( ). Stochastic Epidemic Models : A survey.
Mathematical Biosciences , − . http://dx.doi.org/10.1016/j.mbs.2010.01.006[7] V. Colizza and A. Vespignani. ( ) Epidemic Modeling in Metapopulation Systemwith Heterogeneous Coupling Pattern: Theory and Simulations. Journal of Theoret-ical Biology
251 450 − . http://dx.doi.org/10.1016/j.jtbi.2007.11.028[8] B. J. Coburn, B. G. Wagner and S. Blower. ( ). Modeling Influenza Epidemicsand Pandemics: Insights into the Future of Swine Flu (H1N1). BMC Medicine , , . http://dx.doi.org/10.1186/1741-7015-7-30[9] O. Diekmann, H. Heesterbeek and T. Britton. ( ). Mathematical Tools for Un-derstanding Infectious Disease Dynamics . Princeton University Press.[10] O. Diekmann, H. Heesterbeek and M. G. Roberts. ( ). The Construction of Next-Generation Matrices for Compartmental Epidemic Models.
J. R. Soc. Interface , − . http://dx.doi.org/10.1098/rsif.2009.0386[11] D. Greenhalgh. ( ). Hopf Bifurcation in Epidemic Models with a Latent Periodand Non permanent Immunity. Mathematical and Computer Modelling , . : − . http://dx.doi.org/10.1016/S0895-7177(97)00009-5[12] D. R. Haaneim and F. M. Stein. ( ). Methods of Solution of the RiccatiDifferential Equation. Mathematics Magasine vol. , No. , pp. − .http://dx.doi.org/10.2307/2688697[13] P. Haccou, P. Jagers and V. A. Valutin. ( ). Branching Processes: Variation,Growth, and Extinction of Populations . Cambridge University Press, Cambridge.[14] M. W. Hirsh, S. Smale and R. L. Devaney. ( ). Differential Equations, DynamicalSystems and an Introduction to Chaos . Elsivier Academic Press.[15]
Le Burkina en Chiffres, Ed. . Institut National de la Statistique et de laDémographie.[16] P. Jagers. ( ). Branching Processes with Biological Applications . New York:Wiley.[17] A. Kumar, P. K. Srivastava. ( ). Vaccination and Treatment as Con-trol Interventions in an Infectious Disease Model with their Cost Optimization.
Communications in Nonlinear Science and Numerical Simulation − .http://dx.doi.org/10.1016/j.cnsns.2016.08.005[18] J. C. Miller. ( ). Epidemic Size and Probability in Populations with Het-erogeneous Infectivity and Susceptibility. Physical Review E. : R ) .https://doi.org/10.1103/PhysRevE.76.010101[19] C. E. Mills, J. M. Robins and M. Lipsitch. ( ). Transmissibility of PandemicInfluenza.
Nature , Vol. . http://dx.doi.org/10.1038/nature030632820] E. J. Routh. ( ). A Treatise on the Stability of a Given State of Motion:Particularly Steady Motion . Macmillan and Company.[21] World Health Organization (WHO) (2003).
Influenza.
Fact sheet N o RevisedMarch2003