# Global Stability Analysis and Parameter Estimation for a Diphtheria Model: A Case Study of an Epidemic in Rohingya Refugee Camp in Bangladesh

Authors: Zahurul Islam; Shohel Ahmed; M. M. Rahman; M. F. Karim; M. R. Amin
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

In this article, we have developed a deterministic Susceptible-Latent-Infectious-Recovered (SLIR) model for diphtheria outbreaks. Here, we have studied a case of the diphtheria outbreak in the Rohingya refugee camp in Bangladesh to trace the disease dynamics and find out the peak value of the infection. Both analytical and numerical investigations have been performed on the model to find several remarkable behaviors like the positive and bounded solution, basic reproductive ratio, and equilibria such as disease extinction equilibrium and disease persistence equilibrium which are characterized depending on the basic reproductive ratio and global stability of the model using Lyapunov function for both equilibria. Parameter estimation has been performed to determine the values of the parameter from the daily case data using numerical technique and determined the value of the basic reproductive number for the outbreak asR0=5.86.


---

## Body

## 1. Introduction

Diphtheria is a rapidly spreading disease which is generated byCorynebacterium diphtheriae. Diphtheria transmits in the populations, usually through respiratory droplets, like coughing or sneezing [1, 2]. When the bacteria release the poison or toxin into the body, then the actual disease appears. Fever and throat bruises are the initial symptoms of diphtheria. Besides, a thick grey layer induces the “croup,” which can block the airway and cause a barking cough. Anyone can be infected by diphtheria, but 5-7-aged children who did not receive the appropriate vaccine are usually infected [3–6]. During 1990-1995, above cases 140,000 and 4000 mortalities have been recorded worldwide through the Regional Office of World Health Organization (WHO) for Europe [7–9]. Nowadays, diphtheria is a rare outbreak in the developed world. However, in 2017, several diphtheria outbreaks occurred in Yemen and refugee camps in Bangladesh [10, 11]. In the Rohingya refugee camp in Cox’s Bazar, Bangladesh, a massive-scale diphtheria pestilence was reported. Until December 26, 2017, there were an aggregate number of 2,526 cases and 27 mortalities [12]. There are diphtheria antitoxins in diphtheria treatments to stop poisons from the bacteria and antitoxins to kill the bacteria. The best way to repel diphtheria is through vaccinations [3, 4, 13]. The three shots of the diphtheria-tetanus-pertussis (DTP) vaccine were applied in massive levels to children to control the diphtheria outbreak. To break the transmission chains of the diphtheria outbreak in the Rohingya refugee camp in Bangladesh, emergency vaccination has been applied to children since December 12, 2017, and at the end of 2017, above 90% overall coverage [12].Many researchers have studied epidemic or pandemic disease using mathematical techniques such as Wu and Zhao [14] who have mathematically analyzed an age-structured epidemic model of HIV/AIDS with HAART and spatial diffusion. In a discrete-time SIVS model with saturation incidence rate, Parsamanesh and Erfanian [15] investigated the stability and bifurcations. To examine the impact of an environmental toxin on the spread of infectious illnesses in the population, Saha and Samanta [16] used a toxin-dependent dynamical model. In a discrete time epidemic model including vaccination and vital dynamics, Parsamanesh et al. [17] investigated the stability and bifurcations. Kabir et al. [18] have analyzed the effect of border enforcement measures and socioeconomic cost in export-importation epidemic dynamics using game theory. In a random environment, Samanta and Bera [19] looked at a dynamical model of Chlamydia illness with changing total population size, bilinear incidence rate, and pulse vaccination approach. Parsamanesh and Erfanian [20] looked at the global dynamics of a model with a standard incidence rate and immunization approach. Shahrear et al. [21] have predicted and mathematically analyzed the COVID-19 outbreak in Bangladeshi scenario, and Maugeri et al. [22] have analyzed the transmission of the COVID-19 pandemic in Saudi Arabia and Indonesia. By eliciting behavioural reactions in the community, Saha et al. [23] explored an epidemic model of the COVID-19 outbreak. Liu and Zhang [24] have analyzed global stability for a tuberculosis model. Gao and Huang [25] have investigated a tuberculosis model with optimal control.Some of the researchers have also analyzed the diphtheria epidemic, such as Vitek and Wharton [26] who have studied the potential of the reemergence of diphtheria and other vaccine-preventable infections. Zakikhany and Efstratiou [27] analyzed the current problems and new challenges of diphtheria in Europe. Torrea et al. [28] have studied the diphtheria outbreak with the SIRM model. Ilahi and Widiana [29] have developed an SEIR model for the diphtheria outbreak and analyze vaccination’s effectiveness against the outbreak. Matsuyama et al. [30] have analyzed the sensitivity and ambiguity based on the basic reproductive ratio R0 of the diphtheria epidemic in the Rohingya refugee camp in Bangladesh.Due to the vulnerability of diphtheria epidemics in a confined area, we propose a controlled Susceptible-Latent-Infectious-Recovered (SLIR) model, which is an extension of the simple Susceptible-Infectious-Recovered (SIR) model by adjoining a compartment (L) that tracks the latent people in the cohort. Analytical analysis of the proposed model is performed to prove the existence, uniqueness, positivity, and bounds of the solution. Equilibria of the system and the basic reproductive ratio are also evaluated, and the global stability of the model is proven depending on the basic reproductive ratio. To illustrate the disease dynamics, parameter values are estimated from the daily case data of the outbreak in the Rohingya refugee camp in Bangladesh and found to be the equilibria of the system.


## 2. Mathematical Model

In this section, a mathematical model [31] is developed for the expanse of diphtheria into the populations, which is shown diagrammatically in Figure 1. The entire population at time t is indicated by Nt that is partitioned into four groups: susceptible St, latent (asymptotic) stage Lt, individual affected by diphtheria in the acutely infected stage It, and recovered individuals affected by diphtheria Rt; here, we suppose that the recovered people are not further contagious. Here, λ is a constant that signifies all recruitment that enters the susceptible class, and μ is the natural mortality rate that leaves all classes. The infectious state has an extra mortality rate due to diseases by α, and δ is that rate in which latent infection in people becomes an acute infection. Thus, the people move to state I from state L at a rate of δL. Infectious people are successfully treated with a fixed rate γ, listing to the recovered state. Susceptible people acquire diphtheria infection among active diphtheria at rate βSI, where β signifies the infection transmission coefficient. Moreover,lsignifies a fraction of susceptible people that earn diphtheria infection and migrate to the latent diphtheria stateL, at ratelβSI, and the residual portion,1−l, departs to the active diphtheria stateI. Here, the individuals of the latent class are assumed not to transmit infection.Figure 1
Diagram interaction of each compartment.Assembling all the aforenamed suppositions, the model concerning the transmission dynamics of diphtheria is presented by the subsequent system of differential equations:(1)dStdt=λ−βStIt−μSt,dLtdt=lβStIt−μ+δLt,dItdt=1−lβStIt+δLt−μ+γ+αIt,dRtdt=γIt−μRt,Nt=St+Lt+It+Rt,with following subsidiary conditions:
(2)S0=S0>0,(3)L0=L0≥0,(4)I0=I0≥0,(5)R0=R0≥0.


## 3. Some Basic Characteristic of the Model

To retain the model’s biological efficacy, we want to show the existence, positivity, and boundedness of the solutions to the differential equations for all time.Theorem 1 (existence of unique solution).
Suppose thatS0,L0,I0,R0∈ℝ. Then, there exists continuous differentiable functions S,,L,I,R:0,t0⟶ℝ for positive time t0>0 such that the 4-tuple S,L,I,R covers (1) and S,L,I,R0=S0,L0,I0,R0.Proof of Theorem 1.
By Picard-Lindelöf theorem, it is narrated that the initial value problemy′t=gyt,yt0=y0 has a unique solution yt for locally Lipschitz and continuous function g in time t∈t0−ϵ,t0+ϵ, where ϵ>0. As the system (1) is autonomous, it is enough to prove that the function g:ℝ4⟶ℝ4 is locally Lipschitz in y. Here, g is defined as
(6)gy=λ−βSI−μSlβSI−μ+δL1−lβSI+δL−μ+γ+αIγI−μR.
The Jacobian matrix ofg is obtained as
(7)∇gy=−βI−μ0−βS0lβI−μ+δlβS01−lβIδ1−lβS−μ+γ+α000γ−μ.
This Jacobian is linear inℝ4. Thus, ∇gy satisfies the continuity and differentiability for an interval I∈ℝ4. According to the mean value theorem,
(8)gy1−gy2y1−y2≤∇gy∗,where y∗∈I1. By assuming ∇gy∗=M, we obtain gy1−gy2≤My1−y2 for y1,y2∈I1 and thus, gy is bounded locally for each y∈ℝ4. Therefore, for all compact subset of ℝ4, the derivative of g is continuous and bounded and thus, g is locally Lipschitz. Hence, according to the Picard-Lindelöf theorem, the initial value problem y′t=gyt,y0=y0 for t0>0 has a unique solution yt.Theorem 2.
The proposed model (1) is invariant in the nonnegative orthant ℝ+4.Proof.
LetY=S,L,I,RT; then, model (1) will take the form
(9)dYtdt=LY+C,where
(10)L=−βIt+μ000lβIt−μ+δ001−lβItδ−μ+γ+α000γ−μ,and
(11)C=λ000.
Here,C≥0 and in matrix L, all off-diagonal elements are greater than or equal zero. Hence, L is a Metzler matrix and the system (1) is positive invariant in ℝ+4 [32].Theorem 3.
Fort>0, any solution S,L,I,R of the model (1) with condition (2) is positive.Proof.
The R.H.S. of the model (1) is differentiable; therefore, connecting it with Cauchy problem covenants that there exists a unique maximal solution. The solution of the first equation of system (1) can be figured out alternatively as
(12)dStdt+βIt+μSt=λ.
The solution of Equation (12) is
(13)St=S0e−μt+∫0tβIxdx+e−μt+∫0tβIxdx×∫0tλeμy+∫0tβIududy,
for allt>0. Hence, the R.H.S. of Equation (13) is greater than or equal to zero, i.e., St>0 for all t>0. In the same way, the solution of the second, third, and fourth equations of model (1) is of the form
(14)Lt=L0e−μ+δt+e−μ+δt×∫0tlβSyIyeμ+δydy,It=I0e−μ+γ+αt−∫0t1−lβSxdx+e−μ+γ+αt−∫0t1−lβSxdx×∫0tδLyeμ+γ+αy−∫0t1−lβSududy,Rt=R0e−μt+e−μt×∫0tγIyeμydy,respectively. Those solutions show that all Lt,It, and Rt are greater than or equal zero ∀t>0.Theorem 4 (boundedness).
Suppose the model (1) satisfies S0>0,L0>0,I0>0, and R0>0 and has a unique solution on 0,t0 for some t0>0 by Theorem 1; then, the state functions St,Lt,It, and Rt will be bounded and be positive ∀t∈0,t0.Proof.
Initially, suppose that the values ofSt,Lt,It, and Rt are positive. From Theorem 1, for t>0, there exists a solution on 0,t. Now, denote the largest time by T∗ at which all the populations are positive, or
(15)T∗=supt>0:Ss,Ls,Is,Rs>0,∀s∈0,t.
Since all initial conditions are nonnegative and the solutions are continuous, hence, the solutions must be positive on an interval which is denoted asT∗>0. Therefore, we calculate each term on 0,T∗: instantly, the lower bounds on L,I, and R can be placed.
(16)dLtdt=lβStIt−μ+δLt≥−μ+δLt,as the reduction expressions are linear; this achieves
(17)dLtLt≥−μ+δdt,or
(18)lnLt+lnC≥−μ+δt,or
(19)Lt≥Ce−μ+δt.
Applying initial condition, we get(20)L0≥C,⇒Lt≥L0e−μ+δt>0,for t∈0,T∗.
Again,(21)dItdt=1−lβStIt+δLt−μ+γ+αIt≥−μ+γ+αIt,as the reduction expressions are linear; this achieves
(22)It≥I0e−μ+γ+αt>0,∀t∈0,T∗.
Further,(23)dRtdt=γIt−μRt≥−μRt,i.e.,
(24)Rt≥R0e−μt>0,∀t∈0,T∗.
Similarly, by placing the upper bound ondS/dt, we get
(25)dStdt=λ−βStIt−μSt≤λ,i.e.,
(26)St≤S0+λt≤C1+t,where C is an arbitrary constant which is depending on the upper bound of S0 and λ. Now, by adding the equations forL,I,andRand placing the bounds on this sum and by the positivity of these functions, for the upper bound ofSt, we get
(27)ddtL+I+R=βStIt−μLt−μ+αIt−μRt≤βC1+tIt+μLt+μ+αIt+μRt≤C11+tL+I+R,where C1≥maxβC,μ,μ+α, i.e.,
(28)L+I+Rt≤C2et2,where the constant C2>0 for t∈0,T∗ that only depends on L0,I0,R0, and C1. For the positivity of Lt,It, and Rt are positive, an upper bound can be placed on both L,I, and R by
(29)C2et2≥L+I+Rt≥Lt,C2et2≥L+I+Rt≥It,C2et2≥L+I+Rt≥Rt.
Now,St can be bounded from below using
(30)dSdt=λ−βSI−μS≥−βSI−μS≥−μS−βC2et2S,≥−C31+et2S,where C3≥maxβC2,μ,⇒dS/dt+C31+et2S≥0, i.e.,
(31)St≥S0e−C3∫0t1+eκ2dκ>0.
Therefore,S,L,I, and R remain rigorously positive ∀t∈0,T∗. Hence, there exists a t>T∗ for the continuity, at which the state variables St,Lt,It, and Rt are still positive, which contradicts with the definition of T∗ and specifies that St,Lt,It, and Rt are rigorously positive on the whole interval 0,t. Moreover, all functions remain bounded with this interval; thus, the existing interval can be further extended. Actually, the bounds on S,L,I, and R obtained earlier exist on each compact time interval. For the extension of the time interval to 0,t∀t>0 at which the solution endures and of the above discussion, the solutions continue both positive and bounded on 0,t.


## 4. Equilibria of the System

In this section, we trace the presence of steady states for the dynamical system of nonlinear ODEs (1), describing the Diphtheria disease dynamics. These steady states can be obtained by placing the R.H.S. of (1) to zero; we obtain
(32)λ−βSI−μS=0,(33)lβSI−μ+δL=0,(34)1−lβSI+δL−μ+γ+αI=0,(35)γI−μR=0.Moreover, by solving the above equations, we have found two biologically meaningful equilibrium points. We can classify these two points to be while the infection is either terminated from populations, i.e.,L=I=R=0, or insists in the populationsL≠0,I≠0,R≠0astgrows large.We start to determine the equilibria from the nonlinear intercommunicated terms into Equations (33), (34), and (35) that give
(36)I1−lμ+δβS+δlβS−μ+δμ+γ+α=0.Thus, eitherI=0 or S=μ+δμ+γ+α/1−lμ+δβ. Using I=0 in Equations (33), (34), and (35), we get the disease extinction equilibrium point as
(37)E0=S0,L0,I0,R0=λμ,0,0,0.By settingS=μ+δμ+γ+α/1−lμ+δβ into Equations (32) and (35) yields the infectious persistence equilibrium that exists at the point
(38)E∗=S,L,I,R=μ+δμ+γ+α1−lμ+δβ,lλμ+δ−lμμ+γ+αβδ+1−lμ,λδ+1−lμμ+δμ+γ+α−μβ,γλδ+1−lμμμ+δμ+γ+α−γβ.In the biological sense,E0 is defined as a disease extinction equilibrium point in which an infection survives for a short time and then is naturally dispelled from the populations. The infection is not insisted. The other case, in which the system incline towards E∗, denoted that the populations are impotent to remove the disease spontaneously. If it closes up this remaining fact, then after a particular period, the diphtheria disease model fails its pertinency as it gets broader to keep up the populations.


## 5. Basic Reproductive Ratio

The basic reproductive ratio is also called basic reproductive rate or basic reproduction number and is denoted byR0. It is a significant threshold value generated in epidemiology to mathematically identify the doubt of an infectious disease. This quantity represents the average number of infected persons generated by one infected person introduced into an entirely uninfected susceptible population. We use the next-generation method [33, 34] to obtain the basic reproductive ratio R0.Using the next-generation matrix method on the model (1), we get
(39)F=0lβS001−lβS0,V=μ+δ0−δμ+γ+α.Therefore, we have,(40)FV−1=lδβS0μ+γ+αμ+δlβS0μ+γ+α1−lδβS0μ+γ+αμ+δ1−lβS0μ+γ+αμ+δ.Thus, the spectral radius of the matrixFV−1and the basic reproductive ratioR0are obtained [35].
(41)R0=βδ+μ−lμS0α+γ+μδ+μ.PuttingS0=λ/μ, we obtain,
(42)R0=λβδ+1−lμμα+γ+μδ+μ.This expression ofR0 represents the basic reproductive ratio for the model (1).Remark 5.
The infectious equilibrium point with the expression of basic reproduction numberR0(43)S∗,L∗,I∗,R∗=λμR0,lλR0−1μ+δR0,μβR0−1,γβR0−1.


## 6. Global Stability Analysis

### 6.1. Global Stability at Infectious Extinction Equilibrium

For disease extinction equilibriumE0=S0,L0,I0,R0=λ/μ,0,0,0, we assume the following Lyapunov function:
(44)Ut=S0StS0−1−lnStS0+δ1−lμ+δLt+μ+δ1−lμ+δIt.By differentiation, we get(45)dUdt=1−S0SS′+δ1−lμ+δL′+μ+δ1−lμ+δI′.Substituting the values ofS′,L′, and I′ in the above equation, we have
(46)dUdt=1−S0Sλ−βSI−μS+δ1−lμ+δlβSI−μ+δL+μ+δ1−lμ+δ1−lβSI+δL−μ+γ+αI=λ−μS1−S0S+βS0I−μ+δμ+γ+αμ1−l+δI.After substituting the value ofS0=λ/μ, we are left with
(47)dUdt=−λ−μS2μS+μ+δμ+γ+αμ1−l+δR0−1.At the disease extinction equilibriumE0, the basic reproductive ratio R0≤1, and for all positive values of S,L,I, and R, it is clear that dU/dt≤0. Hence, using LaSalle’s Invariance Principle [36], it is concluded that the model (1) is globally asymptotically stable.Lemma 6.
The infectious extinction equilibriumE0 of the model (1) is globally asymptotically stable when R0≤1, and the disease is naturally dispelled from the populations.


### 6.2. Global Stability at Infectious Persistence Equilibrium

Since none of the state variables are zero at the infectious persistence equilibriumE∗=S∗,L∗,I∗,R∗, thus a Lyapunov function is assumed as
(48)Ut=S−S∗−S∗lnSS∗+B1L−L∗−L∗lnLL∗+B2I−I∗−I∗lnII∗+B3R−R∗−R∗lnRR∗,where B1,B2, and B3 are all nonnegative constants to be obtained. This kind of Lyapunov function has been studied in [37–40].The infectious persistence equilibriumE∗=S∗,L∗,I∗,R∗ satisfies the following equations:
(49)λ=βS∗I∗+μS∗,(50)μ+δL∗=lβS∗I∗,(51)μ+γ+α=1−lβS∗I∗+δL∗,(52)μR∗=γI∗.Now, differentiateU with respect to time t,
(53)U′=1−S∗SS′+B11−L∗LL′+B21−I∗II′+B31−R∗RR′=1−S∗SβS∗I∗+μS∗−βSI−μS+B11−L∗LlβSI−B1μ+δL+B1μ+δL∗+B21−I∗I1−lβSI+δL−B2μ+γ+αI+B2μ+γ+αI∗+B31−R∗RγI−B3μR+B3μR∗,which can be further simplified to
(54)U′=−μS−S∗2S+βS∗I∗1−S∗S+SI−β+B1lβ+B21−lβ+I−B2μ+γ+α+B3γ+βS∗+L−B1μ+δ+B2δ+R−B3μ−B1lβSIL∗L+B1lβS∗I∗−B21−lβSI∗−B2δLI∗I+B21−lβS∗I∗+B2δL∗−B3γIR∗R+B3γI∗.For the positive constantsB1,B2, and B3, the coefficients of SI,I,L, and R must be zero, that is,
(55)−β+B1lβ+B21−lβ=0,(56)−B2μ+γ+α+B3γ+βS∗=0,(57)−B1μ+δ+B2δ=0,(58)−B3μ=0.By solving, the above equation (55) yields
(59)B1=δμ+δB2,B2=μ+δ1−lμ+δ,B3=0.For advantage, we set up new variablesx=S/S∗,y=L/L∗,z=I/I∗, and u=R/R∗ to seek S,L,I, and R and setting the expressions of B1,B2, and B3 in Equation (54), we have
(60)U′=−μS−S∗2S+B21−lβS∗I∗2−1x−x+B1lβS∗I∗2−1x−xzy+B2δL∗1−yz+B3γI∗1−zu.Multiplying byB1 to the 2nd equation of (49) and the 3rd equation of (55) by L∗ yields
(61)B1μ+δL∗=B1lβS∗I∗,B1μ+δL∗=B2δL∗.Hence, it follows that(62)−B1lβS∗I∗+B2δL∗=0.Multiplying byF1Xto the last equation, whereF1Xis considered as a general function that will be determined later andX=x,y,z,u, yields
(63)−B1lβS∗I∗F1X+B2δL∗F1X=0Multiplying the4th equation of (49) by B3 and the 4th equation of (55) by R∗ yields
(64)B3μR∗=B3γI∗,B3μR∗=0.Hence, it follows that(65)B3γI∗=0.Multiplying byF2Xto the last equation, whereF2Xis considered as a general function that will be determined later andX=x,y,z,u, yields
(66)B3γI∗F2X=0.From (54) using (63) and (66) yields
(67)U′=−μS−S∗2S+B21−lβS∗I∗2−1x−x+B1lβS∗I∗2−1x−xzy−F1X+B2δL∗1−yz+F1X+B3γI∗1−zu+F2X.Now, the functionsF1X and F2X are taken so that the coefficients of L∗ and I∗ are zero. For these cases, we get
(68)F1X=yz−1,and
(69)F2X=zu−1.Then, Equation (67) becomes
(70)U′=−μS−S∗2S+B21−lβS∗I∗2−1x−x+B1lβS∗I∗2−1x−xzy−yz+1=−μS−S∗2S+B21−lβS∗I∗2−x−1x+B1lβS∗I∗3−1x−yz−xzy.By the arithmetic mean-geometric mean inequality, for equality, if and only ifS=S∗andy=z=u, the last expression must be less than or equal to zero. Thus, we have U′≤0 with equality if and only if S=S∗ and L/L∗=I/I∗=R/R∗. By LaSalle’s Invariance Principle [36], for each solution, the omega-limit set remains in an invariant set that is contained in Ω=S,L,I,R:S=S∗,L/L∗=I/I∗=R/R∗. Since S must be in S∗,S′ turns zero, which implies that I=I∗, L=L∗, and R=R∗. Thus, there is only invariant set in Ω which is singleton E1. For each solution that intersects, ℝ+04L=I=R=0 limits to E1, which concludes that the disease persistence equilibrium E∗ of (1) is globally asymptotically stable in ℝ+04L=I=R=0 [24].Lemma 7.
The infectious persistence equilibriumE∗ of the model (1) is globally asymptotically stable when R0>1, and the disease persists in the populations for a long time.


## 6.1. Global Stability at Infectious Extinction Equilibrium

For disease extinction equilibriumE0=S0,L0,I0,R0=λ/μ,0,0,0, we assume the following Lyapunov function:
(44)Ut=S0StS0−1−lnStS0+δ1−lμ+δLt+μ+δ1−lμ+δIt.By differentiation, we get(45)dUdt=1−S0SS′+δ1−lμ+δL′+μ+δ1−lμ+δI′.Substituting the values ofS′,L′, and I′ in the above equation, we have
(46)dUdt=1−S0Sλ−βSI−μS+δ1−lμ+δlβSI−μ+δL+μ+δ1−lμ+δ1−lβSI+δL−μ+γ+αI=λ−μS1−S0S+βS0I−μ+δμ+γ+αμ1−l+δI.After substituting the value ofS0=λ/μ, we are left with
(47)dUdt=−λ−μS2μS+μ+δμ+γ+αμ1−l+δR0−1.At the disease extinction equilibriumE0, the basic reproductive ratio R0≤1, and for all positive values of S,L,I, and R, it is clear that dU/dt≤0. Hence, using LaSalle’s Invariance Principle [36], it is concluded that the model (1) is globally asymptotically stable.Lemma 6.
The infectious extinction equilibriumE0 of the model (1) is globally asymptotically stable when R0≤1, and the disease is naturally dispelled from the populations.


## 6.2. Global Stability at Infectious Persistence Equilibrium

Since none of the state variables are zero at the infectious persistence equilibriumE∗=S∗,L∗,I∗,R∗, thus a Lyapunov function is assumed as
(48)Ut=S−S∗−S∗lnSS∗+B1L−L∗−L∗lnLL∗+B2I−I∗−I∗lnII∗+B3R−R∗−R∗lnRR∗,where B1,B2, and B3 are all nonnegative constants to be obtained. This kind of Lyapunov function has been studied in [37–40].The infectious persistence equilibriumE∗=S∗,L∗,I∗,R∗ satisfies the following equations:
(49)λ=βS∗I∗+μS∗,(50)μ+δL∗=lβS∗I∗,(51)μ+γ+α=1−lβS∗I∗+δL∗,(52)μR∗=γI∗.Now, differentiateU with respect to time t,
(53)U′=1−S∗SS′+B11−L∗LL′+B21−I∗II′+B31−R∗RR′=1−S∗SβS∗I∗+μS∗−βSI−μS+B11−L∗LlβSI−B1μ+δL+B1μ+δL∗+B21−I∗I1−lβSI+δL−B2μ+γ+αI+B2μ+γ+αI∗+B31−R∗RγI−B3μR+B3μR∗,which can be further simplified to
(54)U′=−μS−S∗2S+βS∗I∗1−S∗S+SI−β+B1lβ+B21−lβ+I−B2μ+γ+α+B3γ+βS∗+L−B1μ+δ+B2δ+R−B3μ−B1lβSIL∗L+B1lβS∗I∗−B21−lβSI∗−B2δLI∗I+B21−lβS∗I∗+B2δL∗−B3γIR∗R+B3γI∗.For the positive constantsB1,B2, and B3, the coefficients of SI,I,L, and R must be zero, that is,
(55)−β+B1lβ+B21−lβ=0,(56)−B2μ+γ+α+B3γ+βS∗=0,(57)−B1μ+δ+B2δ=0,(58)−B3μ=0.By solving, the above equation (55) yields
(59)B1=δμ+δB2,B2=μ+δ1−lμ+δ,B3=0.For advantage, we set up new variablesx=S/S∗,y=L/L∗,z=I/I∗, and u=R/R∗ to seek S,L,I, and R and setting the expressions of B1,B2, and B3 in Equation (54), we have
(60)U′=−μS−S∗2S+B21−lβS∗I∗2−1x−x+B1lβS∗I∗2−1x−xzy+B2δL∗1−yz+B3γI∗1−zu.Multiplying byB1 to the 2nd equation of (49) and the 3rd equation of (55) by L∗ yields
(61)B1μ+δL∗=B1lβS∗I∗,B1μ+δL∗=B2δL∗.Hence, it follows that(62)−B1lβS∗I∗+B2δL∗=0.Multiplying byF1Xto the last equation, whereF1Xis considered as a general function that will be determined later andX=x,y,z,u, yields
(63)−B1lβS∗I∗F1X+B2δL∗F1X=0Multiplying the4th equation of (49) by B3 and the 4th equation of (55) by R∗ yields
(64)B3μR∗=B3γI∗,B3μR∗=0.Hence, it follows that(65)B3γI∗=0.Multiplying byF2Xto the last equation, whereF2Xis considered as a general function that will be determined later andX=x,y,z,u, yields
(66)B3γI∗F2X=0.From (54) using (63) and (66) yields
(67)U′=−μS−S∗2S+B21−lβS∗I∗2−1x−x+B1lβS∗I∗2−1x−xzy−F1X+B2δL∗1−yz+F1X+B3γI∗1−zu+F2X.Now, the functionsF1X and F2X are taken so that the coefficients of L∗ and I∗ are zero. For these cases, we get
(68)F1X=yz−1,and
(69)F2X=zu−1.Then, Equation (67) becomes
(70)U′=−μS−S∗2S+B21−lβS∗I∗2−1x−x+B1lβS∗I∗2−1x−xzy−yz+1=−μS−S∗2S+B21−lβS∗I∗2−x−1x+B1lβS∗I∗3−1x−yz−xzy.By the arithmetic mean-geometric mean inequality, for equality, if and only ifS=S∗andy=z=u, the last expression must be less than or equal to zero. Thus, we have U′≤0 with equality if and only if S=S∗ and L/L∗=I/I∗=R/R∗. By LaSalle’s Invariance Principle [36], for each solution, the omega-limit set remains in an invariant set that is contained in Ω=S,L,I,R:S=S∗,L/L∗=I/I∗=R/R∗. Since S must be in S∗,S′ turns zero, which implies that I=I∗, L=L∗, and R=R∗. Thus, there is only invariant set in Ω which is singleton E1. For each solution that intersects, ℝ+04L=I=R=0 limits to E1, which concludes that the disease persistence equilibrium E∗ of (1) is globally asymptotically stable in ℝ+04L=I=R=0 [24].Lemma 7.
The infectious persistence equilibriumE∗ of the model (1) is globally asymptotically stable when R0>1, and the disease persists in the populations for a long time.


## 7. Parameter Estimation

In this section, we obtain the value of the unknown parameters for the model (1). To estimate parameter values, we have assumed the initial condition of the state variables as S0,L0,I0,R0=10000,0,1,0. There are seven parameters in our model which are to be obtained. Among these parameters, natural mortality rateμis estimated as 0.002; the recruitment of susceptible classλ=μS0=20; the rate which leavesLtforIt, i.e., incubation periodδ=1/7; and the fraction ofStwhich moves toLt;l=0.95; and disease-induced mortality rate is estimated asα=0.0054. These are derived from the data in the literature [30]. And the rest of the parameters are disease transmission rate β and the recovered rate γ which have to be fitted; therefore, θ=β,γ. Consider the initial value of the parameters to be ω0=λ,μ,α,l,β,δ,γ=20,0.002,0.0054,0.95,0.0000065,1/7,0.005, and the initial condition of the state variables is S0,L0,I0,R0=10000,0,1,0. Using the initial value of the parameters and the initial conditions of the state variables, the value of the unknown parameters is fitted to the model (1) with the help of the nonlinear least square (NLS) method. Table 1 contains the description and estimated or the best fitted values of the parameters. Here, we have simulated the cumulative value of the daily case data, which are illustrated in Figures 2 and 3 that also represent the population dynamic of the susceptible, latent, and infected population St,Lt, and It, respectively. From these figures, it is observed that the infected population (I-class) increases significantly upon the infection and arrives at the peak at the 36th day I43=3.126×103; after that, it is decaying.Table 1
Description and value of the parameters of the diphtheria model.
ParameterDescriptionValueSourceλThe recruitment of susceptible class200personsday−1EstimatedμNatural mortality rate0.002day−1EstimatedαDisease induced mortality rate0.0054day−1EstimatedβDisease transmission rate0.000097persons−1day−1FittedlThe fraction ofSt which moves to Lt0.95EstimatedγRecovered rate0.156day−1FittedδThe rate which leavesL(t) for It0.143day−1EstimatedFigure 2
The diphtheria model (1) simulation in log scale.
(a)(b)(c)Figure 3
The fitted diphtheria model (1).
(a)(b)(c)


## 8. Numerical Results

To further investigate the behaviour of the model (1), we conducted various numerical investigations applying the estimations that are gained and given in Table 1. For this intention, we consider two parameter sets resembling the cases of stability of the infectious persistence equilibrium, where R0>1, and disease extinction steady state, where R0<1. The outcomes obtained for both equilibria with stability analysis are also numerically demonstrated using MATLAB R2018a.Using the parameter values from Table1, the basic reproductive ratio becomes R0=5.86>1 thereby signifying the asymptotic stability of the infected steady state. For this reason, different initial conditions of S0,L0,I0,R0 are chosen as IC1=10000,0,1,0,IC2=8000,0,2,0, and IC3=12000,0,3,0.Figure4 illustrates the system dynamics of the susceptible, latent, and infected population for the three initial conditions within two years, i.e., 730 days. In Figure 4(a), the susceptible population decays very sharply and reaches the nadir at 189, 283, and 134 for IC1, IC2, and IC3, respectively. As time increases, they are again increasing together and reaching a peak point of approximately 2782. Again, it is decreasing and reaches another nadir at 1525. Further, it is increasing and asymptotically stable at 1706 within two years; i.e., susceptible population would be constant. In Figure 4(b), the latent population increases sharply and reaches the first peak points 3136, 2200, and 4161 for IC1, IC2, and IC3, respectively; then, they are decreasing sharply and reach a nadir at 5 together within 3.67 months and stable about three months. As time increases, they are again increasing and reach the second peak at 305 within the next 3 months. Again, they are decaying and reach another nadir at 51 within the next 3.67 months. Further, they increase and reach the third peak point of 158 within the next 4 months. They are decaying further and reach another nadir at 85 within the next 4 months. As time increases, they are increasing further and asymptotically stable at 109 within 2 years. Moreover, in Figure 4(c), the infected population increases very sharply and reaches the first peak points at 2383, 1739, and 3058 for the same initial conditions; then, they are decaying as they are increased and reach a nadir at 5 together within 4 months and stable about three months. As time increases, they are again increasing and reaches the second peak at 280 within the next 3.33 months. Again, they are decaying and reach another nadir at 47 within the next 3.67 months. Further, they are increasing and reach the third peak at 145 within the next 4 months. They are decaying further and reach another nadir at 80 within the next 3.67 months. As time increases, they are increasing further and asymptotically stable at 100 within 2 years.Figure 4
Population dynamics interaction betweenSt, Lt, and It of diphtheria model (1) when R0=5.86>1 for different initial conditions.
(a)(b)(c)Figure5 illustrates the system’s phase portrait for different initial conditions. It represents the relative change of the susceptible St, latent Lt, and infected populations It to one another over time by a single trajectory. It also characterises the stability of the system. For different initial conditions, the trajectories are approaching a single point which specifies the disease persistence equilibrium point E∗=1706,109,100,7814 when the basic reproduction number R0=5.86>1. In this case, the trajectories approach the long-term steady state, and the disease persists in the populations for t⟶∞.Figure 5
System’s phase portrait of diphtheria model (1) in 2D and 3D when R0=5.86>1 for different initial conditions.
(a)(b)(c)For disease-extinction equilibrium, we assume the value of infection transmission rateβ different from Table 1 as β=0.000005. Therefore, the basic reproductive ratio is evaluated as R0=0.302<1. For this case, the disease dynamics are illustrated in Figure 6 for the same initial conditions. The latent and infected populations are converged to 0 within 4 months that are illustrated in Figures 6(b) and 6(c), respectively, which indicates that the disease will be extincted from the populations by itself within 120 days. However, in Figure 6(a), in the susceptible population, only positive values remain for the different initial conditions that indicate the infection-free steady state. Moreover, Figure 7 illustrates the infection-free steady states and the interaction between the populations by three trajectories for three different initial conditions. The trajectories are approaching a single point defined as infection-free steady-state E0=10000,0,0,0 and remain at this point for t⟶∞.Figure 6
Population dynamics interaction betweenSt, Lt, and It of diphtheria model (1) when R0=0.302<1 for different initial conditions.
(a)(b)(c)Figure 7
System’s phase of diphtheria model (1) in 3D when R0=0.302<1 for different initial conditions.


## 9. Conclusion

We have proposed a diphtheria epidemic model and found two steady-state equilibrium points: one is disease extinction equilibrium pointE0 (37), and another is infectious persistence equilibrium point E∗ (38). We have formulated the basic reproductive number in terms of parameters. We have also shown analytically that the infectious extinction equilibrium (37) is globally asymptotically stable when the basic reproductive ratio R0 does not exceed unity; the infection is dispelled by itself from the populations. The infectious persistence equilibrium (38) is also globally asymptotically stable when the basic reproductive ratio R0 is more than unity; the disease persists in the populations at a certain level. We have fitted the daily case data from November 8, 2017, to December 27, 2017, given in [30] to our model and have evaluated the parameter’s value. Both equilibria have been analyzed numerically and found a lot that matches the real scenario. We have found the first peak at 38 days for IC1 that matches with the real data, and the second and third peaks have been found at 310 days and 1.5 years, respectively, which also a lot matches with the real scenario, like the highest infection found after one month, given in [30, 41–44]. The enumerated infectious persistence equilibrium E∗ is 1706,109,100,7814 and infection-free steady-state E0 is 10000,0,0,0. A statistical model was used to calculate the numeric value of the basic reproductive ratio R0 in [30] and deduced a range of estimates ranging from 4.7 to 14.8 with a median estimate of 7.2. But in this study, it has calculated R0=5.86 by involving a mathematical model, which indicates that the infection rate is very high. This study suggests applying treatments to control the diphtheria epidemic. Lastly, we hope that this study will be focused on the assumption of control strategies by constituents and policymakers.



---

*Source: 6545179-2022-01-27.xml*