EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 12, No. 3, 2019, 834-845 ISSN 1307-5543 – www.ejpam.com Published by New York Business Global Stochastic approach of epidemic model using the SEIRS Hay Yoba Talkibing1, Barro Diakarya2,∗, Ouoba Fabrice1 1 UFR-SEA, Université Ouaga I, Pr Joseph Ki-Zerbo, Ouagadougou, Burkina Faso 2 UFR-SEG, Université Ouaga 2, Ouagadougou, Burkina Faso Abstract. The study of infectious diseases represents one of the oldest and richest sectors of biomathematics. The transmission dynamics of these diseases are still a major problem in math- ematical epidemiology. In this work, we propose a stochastic version of a SEIRS epidemiological model for infectious diseases evolving in a random environment for the propagation of infectious diseases. This random model takes into account the rates of immigration and mortality in each compartment and the spread of these diseases follows a four-state Markovian process. We first study the stability of the model and then estimate the marginal parameters (means, variances and covariates) of each disease state over time. Real measles data are applied to the model. 2010 Mathematics Subject Classifications: 60G05, 92D30, 60J10 Key Words and Phrases: Stochastic processes, epidemiology, Markov chain, SEIRS 1. Introduction Infectious diseases are of increasing concern to authorities and public health officials [2]. In the face of increasing bacterial resistance, the emergence of new pathogens and the rapid spread of the epidemic make surveillance and prevention of disease transmission par- ticularly important and indispensable. Even with permanent and continuous surveillance of infectious diseases, it must be noted that their etiologies are still largely unknown. To combat this random phenomenon, in recent decades mathematical models have been devel- oped and implemented to study the spread of infectious diseases since the early twentieth century in the field of mathematical epidemiology through the work of Daniel Bernoulli [1] . Stochastic and deterministic models of epidemics allow researchers to gain valuable knowledge about many infectious diseases and to study strategies to combat them. The first stochastic epidemic models were presented by L. Reed and W, H. Frost in [3] to accurately describe the discrete-time spread of disease in a population. These models developed in 1928 and 1931[3] are based on the notion of probability and statistics. Later ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v12i3.3400 Email addresses: talkibingfils@yahoo.fr (H. Y. Talkibing), dbarro2@gmail.com (B. Diakarya), ddifab@yahoo.fr (O. Fabrice) http://www.ejpam.com 834 c© 2019 EJPAM All rights reserved. H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 835 Bartlett [M. Costa 2011] proposed in 1949, the continuous time stochastic SIR model which serves as the basis for stochastic models. These last yeas, we can find in the liter- ature several works but much remains to be done with the emergence of some of the new pathogens. In this work, we are interested in the SEIRS model described in [7] for infectious diseases. We propose a stochastic version of this model using birth and death processes to describe the spread of infectious diseases in a random environment. 2. Stochastic probability of Birth and Mortality In this section, we propose a stochastic version of the SEIRS model (fig 1) reviewed above. So, we assume that the spread of the disease follows a linear process of birth and death in continuous time and in four states in a variable environment. Moreover, we as- sume that the switching between environments follows a homogeneous Markov chain in continuous time. A model of the dynamics of infectious disease transition can be formulated using a stochastic process {Xn, n = 1, 2, 3, ...} as a Markov string with continuous time T and discrete state space X0, X1, ..., Xn satisfying the following Markov property : P (Xn+1 = jn+1/X0 = j0, X1 = j1, ..., Xn = jn) = P (Xn+1 = jn+1/Xn = jn). (1) That is, the evolution of the process at time n+1 depends only on the state of the process at time n but not on the past states 0, 1, 2, ..., n − 1 of the process. In biological modeling, a disease transmission process described with the states X0, ..., Xn can be seen as a Markov chain with a transition matrix Q = (pij) where the probabilities of the transitions are such as : pij = P (Xn+1 = j/Xn = i) for all i, j = 1, 2, 3, ..., n. (2) The population of interest is divided into four compartments where each of them rep- resents a specific stage of the epidemic. We have X1 for susceptible or holy individuals, X2 for lectures, X3 for infected individuals and X4 for healers. Add to this model the apparitions rates σi for i = 1, 2, 3, 4 independent of lambda and the natural death rates µi for all i = 1, 2, 3, 4 independent of the disease in the classes of persons susceptible to X1(t), exposed X2(t), infected X3(t) and healed X4(t) respectively at time t. Moreover, let Z(t) = X1(t) +X2(t) +X3(t) +X4(t) be the population size at a given time t. Assume that all rates are not time-dependent t and disease transmission occurs only by contact between susceptible and infectious individuals in a relatively small time in- terval dt. Once the vaccination is administered to the population, the immunized subjects against the disease and healed subjects enter the compartment X4 and lose their immunity H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 836 rate significance λ the rate of newborns α the rate of attacks of the disease β The infection rate γ recovery rate Θλ the rate of newborns vaccinated and immunized against the disease τ the tau of individuals healed but having lost temporary immunity η the tau of holy individuals vaccinated and immunized against the disease Θ the tau of unvaccinated newborns σi, i = 1.4 rates of immigrant individuals in different compartments µi, i = 1.4 mortality rates of the different compartments independent of the disease Xi, i = 1, 4 The steps in the process against this disease after a given latent period to become susceptible again. Newborn un- vaccinated enter the susceptible compartment X1 with a rate θ = 1− λ and after contact with an infected subject first expose themselves to the disease and become infected, healed and susceptible again according to the diagram 1). In the case of infectious diseases, the ideal objective is to eradicate these diseases completely by preventive measures or by set- ting up a more effective mass vaccination program. Following the SEIR and SEIRS models with vital dynamics, Anderson and May studied vaccinations applied to the population [7]. Figure 1: The propagation cycle of the SEIRS epidemic model Transitions intensities are defined as a probability limit of transitions between times t and t + dt, when dt tends to zero [4]. Under the assumption of the existence of these limits, let’s determine the infinitesimal probability laws of each process over a small time interval dt. • Birth and emigration of S : Once the population are vaccinated, the compartment S contains unimmunized indi- H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 837 viduals who are therefore susceptible to the disease. Among these individuals are immi- grants with a rate of σ1 from another locality and newborns not immune to the disease with a rate of λ. The probability that there will be a birth and then immigration of a susceptible person in the n process state is given by : P  S(t+ dt) = n/S(t) = n− 1 E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′  = [λ(1− θ)(n− 1) + σ1]dt+ 0(dt). • Mortality of S : The mortality rate of compartment S, independent of the cause of the disease, is µ1. The probability of deaths in state n of the process is given by: P  S(t+ dt) = n/S(t) = n+ 1 E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′  = µ1dt+ 0(dt). • Transition from S to E: The individuals in compartment S come into contact with infected people at a rate of α and expose themselves to the disease. The probability of switching from the n state to the m state is given by: P  S(t+ dt) = n/S(t) = n+ 1 E(t+ dt) = m/E(t) = m− 1 I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′  = α(n+ 1)dt+ 0(dt) • Mortality of E : The probability of deaths in the m state of the process is given by : P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m− 1 I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′  = µ2dt+ 0(dt) • Mortality of I: The probability that there will be deaths regardless of the cause of the disease among the infected in the n ′ state of the process is given by: H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 838 P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m+ 1 I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′  = µ3dt+ 0(dt). • Mortality of R : The probability of deaths among people restored to the m ′ state of the process is given by: P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′ + 1  = µ4dt+ 0(dt). 3. Stochastic probability of migration and Transition • Immigration of E: The compartment E contains people exposed to the disease after contact with the patients. The probability of immigration at the m state of the process is therefore: P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m− 1 I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′  = σ2dt+ 0(dt) • Transition from E to I: The exposed individuals E, become infectious with a rate β and the probability of passing from the state m and to the state n ′ is given by: P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m+ 1 I(t+ dt) = n ′ /S(t) = n ′ − 1 R(t+ dt) = m ′ /S(t) = m ′  = β(m+ 1)dt+ 0(dt) • Immigration of I: The probability of immigration to the state n ′ of the process in compartment I is: P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ − 1 R(t+ dt) = m ′ /S(t) = m ′  = σ3dt+ 0(dt) H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 839 • Transition from I to R: Infected individuals I heal with a rate of gammama. The probability of switching from the state n ′ and to the state m ′ is given by : P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ + 1 R(t+ dt) = m ′ /S(t) = m ′ − 1  = γ(n3 + 1)dt+ 0(dt) • Immigration of R: The probability of immigration in compartment R in the m ′ state of the process is: P  S(t+ dt) = n/S(t) = n E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′ − 1  = [σ4 + θλ]dt+ 0(dt) • Transition from R to S: Recovered people become susceptible again with a rate of tau. The probability of switching from the state m ′ and to the state n ′ is therefore: P  S(t+ dt) = n/S(t) = n− 1 E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /R(t) = m ′ + 1  = τ(n4 + 1)dt+ 0(dt). • Transition from S to R: Individuals immunized against the disease with a rate of τ after direct vaccination coverage integrate the R compartments. The probability of transition from the n state to the m ′ state is therefore: P  S(t+ dt) = n/S(t) = n+ 1 E(t+ dt) = m/E(t) = m I(t+ dt) = n ′ /S(t) = n ′ R(t+ dt) = m ′ /S(t) = m ′ − 1  = η(n+ 3)dt+ 0(dt). 4. Stochastic and Linear properties of the transition The Markov property for the law of a multistate process is that the events {X(t) = j}, j = 0, ..., n are independent on the past Xn−1 given X(s). For Markov processes, the transition probability matrices satisfy the Chapman-Kolmogorov property [6]. The set of H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 840 transition intensities forms the matrix of transition intensity Q. The sum of the terms of each line of these matrices is equal to zero. We have the matrix of transmission probabilities Q such as . Q =  −M1 + µ1 + α α 0 η 0 −M2 + µ2 β 0 0 0 −M3 + µ3 γ τ 0 0 −M4 + τ + µ4  Where ⇒  M1 = λ(1− θ) + σ1 + η M2 = α+ β + σ2 M3 = β + γ + σ3 M4 = γ + η + θλ+ σ4 (3) M1, M2, M3 et M4 are positive quantity. The Kolmogorov equation gives the matrix of transition probabilities Q(t) as a function of the exponential of a matrix : P (t) = exp {Qt} . The exponential of a matrix is calculated relatively easily, particularly when the ma- trix can be diagonalized [4]. If Q has eigenvalues which are all different, it can be written as Q = PDP−1, where D is the diagonal matrix of eigen-values, and P is the matrix of eigenvectors. In this case, the probability P (t) at stable condition gives P (t) = PetDP−1 where the exponential of a diagonal matrix D is a diagonal matrix whose diagonal terms are the exponential of terms of D. Using the cumulants generating function [5] , we then calculate the means, variances and covariances of the disease states as a function of time. For this purpose, the probability generating function is defined by : Φ(x1, x2, x3, x4, t) = ∑ n,n ′ ,m,m ′>0 p(n,m, n ′ ,m ′ , t)xn1x m 2 x n ′ 3 x m ′ 4 . H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 841 Using A(t) = (Ai,j(t))16i,j64 = (tAi,j) [i 6, 4j 6], the matrix of births with positive or null coefficients, such as : Aij =  λ(1− θ) + σ1 0 0 0 0 σ2 0 0 0 0 σ3 0 0 0 0 σ4  ; and B(t) = (Bi,j(t)) the matrix of deaths given by: Bij =  µ1 + α+ η −α 0 −η 0 µ2 + β −β 0 0 0 µ3 + γ −γ −τ 0 0 µ4 + τ  ; then, the joint generating function of the moments is defined by: M(a1, a2, a3, a4, t) = ∞∑ n,n′ ,m,m′=0 p(n, n ′ ,m,m ′ , t)ea1nea2n ′ ea3mea4m ′ . (4) By setting xk = eak , we get ∂xk = eak∂ak for all k = 1, 2, 3, 4 Let us define the generating function of the cumulants by: K(a1, a2, a3, a4, t) = ∞∑ n,n′ ,m,m′=1 an1 n! an ′ 2 n′ ! am3 m! am ′ 4 m′ ! knn′mm′ (a1, a2, a3, a4, t) (5) where knn′mm′ (a1, a2, a3, a4, t) are cumulative: knn′mm′ (a1, a2, a3, a4, t) = drK(a1, a2, a3, a4, t) dtr |t=0. (6) Furthermore, the generating function of the cumulants is computed as follows: K(a1, a2, a3, a4, t) = logM(a1, a2, a3, a4, t) (7) And ∂M(a1, a2, a3, a4, t) = eK∂K(a1, a2, a3, a4, t). (8) H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 842 By substituting ea = ∑∞ i=0 ai i! in (5) and then, by derivating the functionsK(a1, a2, a3, a4, t), it comes that: ∂K(a1, a2, a3, a4, t) ∂t = [ ε1 ∞∑ n=0 an1 n! + ε3 ∞∑ n=0 (−a1)n n! − ε2 ]∂knn′mm′ (a1, a2, a3, a4, t) ∂a1 (9) + [ σ2 ∞∑ n′=0 an ′ 2 n′ ! + α ∞∑ n=0 an1 n! ∞∑ n ′ =0 (−a2)n ′ n′ ! − ε4 − ε5 ]∂knn′mm′ (a1, a2, a3, a4, t) ∂a2 + [ σ3 ∞∑ m=0 am3 m! + β ∞∑ n′=0 an ′ 2 n′ ! ∞∑ m=0 (−a3)m m! − ε7 ∞∑ m=0 (−a3)m m! − ε6 ]∂knn′mm′ (a1, a2, a3, a4, t) ∂a3 + [ σ4 ∞∑ m′=0 am ′ 4 m′ ! + η ∞∑ n=0 an1 n! ∞∑ m=0 (−a4)m ′ m′ ! + γ ∞∑ m=0 am3 m! ∞∑ m′=0 (−a4)m ′ m′ ! − ε9 ∞∑ m ′ =0 (−a4)m ′ m′ ! − ε8 ]∂knn′mm′ (a1, a2, a3, a4, t) ∂a4 . Moreover, knowing that: E(X) = ∂knn′mm ′ (a1, a2, a3, a4, t) V (X) = ∂2knn′mm′ (a1, a2, a3, a4, t). . (10) Which gives us: E(Xi(t)) = dknn′mm′ (a1, a2, a3, a4, t) dai |t=0 V (Xi(t)) = d2knn′mm′ (a1, a2, a3, a4, t) da2i |t=0 pour i = 1, 2, 3, 4 (11) Using the moment generating function (5), we can determine the average and variance of each state. The derivative of 5 gives: ∂K(a1, a2, a3, a4, t) ∂t = ∞∑ n,n′ ,m,m′=0 an1 n! an ′ 2 n′ ! am3 m! am ′ 4 m′ ! k ′ n,n′ ,m,m′ (a1, a2, a3, a4, t). (12) H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 843 By equalizing the two previous expressions (12) and (5) and using the relationship (11), we obtain a system of linear differential equations of the first degree whose variables are means of each state of the process defined by the relationship (13). dE(X1) dt = ρ1E(X1) dE(X2) dt = ρ2E(X2) dE(X3) dt = ρ3E(X3) + σ4E(X2X3) + ρ5E(X2) dE(X4) dt = ρ6E(X4) + ρ7E(X2X4)− ε5E(X3) . (13) With ρi, i = 1, ..., 8, the strictly positive constants defined by the system (14) : ρ1 = σ4 + σ2 + η + ε8 − ε5 − ε4 + ε3 − ε2 + γ + α ρ2 = σ4 + η + ε8 − ε4 + ε3 − ε2 + γ + α+ β ρ3 = σ4 + σ3 + σ2 + η − ε9 + ε8 − ε7 + ε6 −−ε5 − ε4 + ε3 − ε2 + γ + α+ β ρ4 = σ3 − ε7 + ε6 − β ρ5 = σ3 + σ7 − β ρ6 = σ4 + σ3 + σ2 + η − ε9 + ε8 − σ7 + ε6 − ε5 − ε4 + ε3 − ε2 + γ + α+ β ρ7 = σ3 +−ε7 + ε6 + β (14) The dispersion in turn of the average of each state of the process is given by the relationship (15). For a statistic with little dispersion, the observations are close to each other, and therefore to their average.  dV (X1) dt = 1 2(σ2 + η + ε2 + ε3 + γ + α)V (X1) dV (X2) dt = 1 2(σ4 + σ2 + η + ε6 + ε3 + ε2 + γ + β)V (X2) dV (X3) dt = 1 2(σ4 + σ3 + σ2 + η + ε9 + ε3 + ε2 + γ + β)V (X3) dV (X4) dt = 1 2(σ4 + σ3 + σ2 + η + ε9 + ε8 + ε7 + +ε6 + +ε5 + ε3 + α)V (X4) . (15) H. Y. Talkibing, B. Diakarya, O. Fabrice / Eur. J. Pure Appl. Math, 12 (3) (2019), 834-845 844 The degree of the links between the different states of the process can be estimated from the co-variables by solving differential equations of the relationship (16).  dCOV (X1X2) dt = (σ4 + σ2 + η − ε9 + ε8 + ε3 − ε2 + γ + α)COV (X1X2) dCOV (X1X3) dt = (σ4 + σ3 + σ2 + η − ε9 + ε8 − ε7 + ε6 − ε5 − ε4 + ε3 − ε2 + γ + α+ β)COV (X1X3) dCOV (X2X3) dt = (σ4 + σ3 + σ2 + η − ε9 + ε8 − ε5 − ε4 + ε3 − ε2 + γ + α+ β)COV (X2X3) dCOV (X2X4) dt = (σ4 + σ2 + η − ε9 − ε5 − ε4 + ε3 − ε2 + γ + α+ β)COV (X2X4) dCOV (X3X4) dt = (σ4 + σ3 + σ2 − ε9 + ε8 + ε7 + ε6 − ε5 − ε3 − ε2 + α+ β)COV (X3X4) dCOV (X1X4) dt = (σ4 + σ2 + η − ε9 − ε5 − ε4 + ε3 − ε2 + γ + α)COV (X1X4) (16) The resolution of the written equation system by the relationships (13), (15) and (16) have as initial conditions:  S(0) = S0 > 0, E(0) = E0 ≥ 0, I(0) = I0 > 0, R(0) = R0 ≥ 0 Φ(x1, x2, x3, x4, 0) = xS0 1 x E0 2 xI03 x R0 4 ,M(a1, a2, a3, a4, 0) = ea1S0ea2E0ea3I0ea4R0 K(a1, a2, a3, a4, 0) = a1S0 + a2E0 + a3I0 + a4R0 E(X1(0)) = S0, E(X2(0)) = E0, E(X3(0)) = I0, E(X4(0)) = R0 . (17) 5. Conclusions In this study, we have proposed the stochastic approach of the deterministic SEIRS model with a vital dynamic and temporary immunity of [7]. whose balance and stability have been studied. The stochastic approach can be used to accurately describe the spread of infectious diseases in a random environment. This SEIRS model generalizes the SI, SIR and SEIR models and is therefore a valuable tool for stochastic modelling of epidemiolog- ical problems. The resulting systems of equations incorporate complexities and assume that randomness is important and explicitly include it in the behaviour of the system com- pared to the deterministic model that assumes that randomness has a negligible effect and considers only the average behaviour of the system. For good monitoring and prevention of transmission of epidemics, the results obtained provide more precision on the random REFERENCES 845 evolution of the disease over time. The main disadvantage of this approach is that in practice, it is often difficult to have data adapted to the model to be implemented. It is also clear that most of the results remain valid not only for a finite number of Markovian environments, but also for more general ergodic environments. References [1] J. Gani D. J. Daley. Epidemic modelling. In Cambridge University Press, editor, An introduction, page 14, 2001. [2] Bibliothque de l’OMS. Modlisation mathmatique pour le contrle et llimination du paludisme. novembre 2010. [3] J. Wu Fred Brauer, Pauline van den Driessche. Lecteure note in mathematics. In Oxford Mathematical Biosciences Subseries: P.K. Maini, editor, Mathematical epi- demiology, 1945. [4] Forrest .W. Crawford Vladimir N. Minin Marc A. Suchard Lam Si Tung Ho, Jasen Ho. Birth/birth-death processes and their cumputable transition probabilities with biolog- ical applications. Cornell University Library, 7 Aug 2017. [5] Abdelkarim Ed-Darraz Nicolas Bacaer. un environnement alatoire. Journal of Math- ematical Biology, Springer Verlag (Germany), 2014. [6] Christelle Verg Pierre Del Moral. Modeles et methodes stochastiques. In Springer, editor, Une introduction avec applications. Comit de Lecture 20122015/Editorial Board 20122015, 2014. [7] Marek B. Trawicki. Cdeterministic seirs epidemic model for modeling vital dynamics, vaccinations, and temporary immunity. Molecular Diversity Preservation Interna- tional, January 2017.