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

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 as ℛ0 = 5.86.


Introduction
Diphtheria is a rapidly spreading disease which is generated by Corynebacterium 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-7aged children who did not receive the appropriate vaccine are usually infected [3][4][5][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][8][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 massivescale 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 R 0 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.

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 NðtÞ that is partitioned into four groups: susceptible ðSðtÞÞ, latent (asymptotic) stage ðLðtÞÞ, individual affected by diphtheria in the acutely infected stage ðIðtÞÞ, and recovered individuals affected by diphtheria ðRð tÞÞ; here, we suppose that the recovered people are not fur-ther 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 stateðLÞ, at ratelβSI, and the residual portion,ð1 − lÞ, departs to the active diphtheria stateðIÞ. Here, the individuals of the latent class are assumed not to transmit infection.
Assembling all the aforenamed suppositions, the model concerning the transmission dynamics of diphtheria is presented by the subsequent system of differential equations: with following subsidiary conditions:

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. Proof of Theorem 1. By Picard-Lindelöf theorem, it is narrated that the initial value problem y ′ ðtÞ = gðyðtÞÞ, yðt 0 Þ = y 0 has a unique solution yðtÞ for locally Lipschitz and continuous function g in time t ∈ ½t 0 − ϵ, t 0 + ϵ, 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 The Jacobian matrix of g is obtained as This Jacobian is linear in ℝ 4 . Thus, ∇gðyÞ satisfies the continuity and differentiability for an interval I ∈ ℝ 4 . According to the mean value theorem, where y * ∈ I 1 . By assuming j∇gðy * Þj = M, we obtain jgðy 1 Þ − gðy 2 Þj ≤ Mjy 1 − y 2 j for y 1 , y 2 ∈ I 1 and thus, gðyÞ 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Þ = gðyðtÞÞ, yð0Þ = y 0 for t 0 > 0 has a unique solution yðtÞ.
Proof. Let Y = ðS, L, I, RÞ T ; then, model (1) will take the form where 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].
The solution of Equation (12) is for all t > 0. Hence, the R.H.S. of Equation (13) is greater than or equal to zero, i.e., SðtÞ > 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 Theorem 4 (boundedness). Suppose the model (1) satisfies S 0 > 0, L 0 > 0, I 0 > 0, and R 0 > 0 and has a unique solution on ½0, t 0 for some t 0 > 0 by Theorem 1; then, the state functions SðtÞ, LðtÞ, IðtÞ, and RðtÞ will be bounded and be positive ∀t ∈ ½0, t 0 .
Proof. Initially, suppose that the values of SðtÞ, LðtÞ, IðtÞ, and RðtÞ 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 Since all initial conditions are nonnegative and the solutions are continuous, hence, the solutions must be positive on an interval which is denoted as T * > 0. Therefore, we calculate each term on ½0, T * : instantly, the lower bounds on L, I, and R can be placed.
as the reduction expressions are linear; this achieves or Applying initial condition, we get for t ∈ ½0, T * : Again,   Computational and Mathematical Methods in Medicine as the reduction expressions are linear; this achieves Further, i.e., Similarly, by placing the upper bound on dS/dt, we get i.e., where C is an arbitrary constant which is depending on the upper bound of Sð0Þ 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 ofSðtÞ, we get where C 1 ≥ max fβC, μ, ðμ + αÞg, i.e., where the constant C 2 > 0 for t ∈ ½0, T * that only depends on Lð0Þ, Ið0Þ, Rð0Þ, and C 1 . For the positivity of LðtÞ, IðtÞ, and RðtÞ are positive, an upper bound can be placed on both L, I, and R by  Now, SðtÞ can be bounded from below using where C 3 ≥ max fβC 2 , μg, ⇒ðdS/dtÞ 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 SðtÞ, LðtÞ, IðtÞ, and RðtÞ are still positive, which contradicts with the definition of T * and specifies that SðtÞ, LðtÞ, IðtÞ, and RðtÞ 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.

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 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 populationsðL ≠ 0, I ≠ 0, R ≠ 0Þastgrows large.
We start to determine the equilibria from the nonlinear intercommunicated terms into Equations (33), (34), and (35) that give Thus, either I = 0 or S = ðμ + δÞ ðμ + γ + αÞ/ðð1 − lÞμ + δ Þβ. Using I = 0 in Equations (33), (34), and (35), we get the disease extinction equilibrium point as By setting S = ðμ + δÞðμ + γ + αÞ/ðð1 − lÞμ + δÞβ into Equations (32) and (35) yields the infectious persistence equilibrium that exists at the point In the biological sense, E 0 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.

Basic Reproductive Ratio
The basic reproductive ratio is also called basic reproductive rate or basic reproduction number and is denoted by R 0 . It   Computational and Mathematical Methods in Medicine 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 R 0 . Using the next-generation matrix method on the model (1), we get Therefore, we have, Thus, the spectral radius of the matrixFV −1 and the basic reproductive ratioR 0 are obtained [35].
Putting S 0 = λ/μ, we obtain, This expression of R 0 represents the basic reproductive ratio for the model (1).

Remark 5. The infectious equilibrium point with the expression of basic reproduction number
6. Global Stability Analysis 6.1. Global Stability at Infectious Extinction Equilibrium. For disease extinction equilibrium E 0 = ðS 0 , L 0 , I 0 , R 0 Þ = ðλ/μ, 0, 0, 0Þ, we assume the following Lyapunov function: By differentiation, we get Substituting the values of S ′ , L ′ , and I ′ in the above equation, we have After substituting the value of S 0 = λ/μ, we are left with At the disease extinction equilibrium E 0 , the basic reproductive ratio R 0 ≤ 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 equilibrium ðE 0 Þ of the model (1) is globally asymptotically stable when R 0 ≤ 1, and the disease is naturally dispelled from the populations.

Global Stability at Infectious Persistence Equilibrium.
Since none of the state variables are zero at the infectious persistence equilibrium E * = ðS * , L * , I * , R * Þ, thus a Lyapunov function is assumed as where B 1 , B 2 , and B 3 are all nonnegative constants to be obtained. This kind of Lyapunov function has been studied in [37][38][39][40]. The infectious persistence equilibrium E * = ðS * , L * , I * , R * Þ satisfies the following equations: which can be further simplified to For the positive constants B 1 , B 2 , and B 3 , the coefficients of SI, I, L, and R must be zero, that is, By solving, the above equation (55) yields For advantage, we set up new variables x = S/S * , y = L/ L * , z = I/I * , and u = R/R * to seek S, L, I, and R and setting the expressions of B 1 , B 2 , and B 3 in Equation (54), we have Multiplying by B 1 to the 2 nd equation of (49) and the 3 rd equation of (55) by L * yields Hence, it follows that Multiplying byF 1 ðXÞto the last equation, whereF 1 ðXÞis considered as a general function that will be determined later andX = ðx, y, z, uÞ, yields Multiplying the 4 th equation of (49) by B 3 and the 4 th equation of (55) by R * yields Hence, it follows that Multiplying byF 2 ðXÞto the last equation, whereF 2 ðXÞis considered as a general function that will be determined later andX = ðx, y, z, uÞ, yields From (54) using (63) and (66) yields Now, the functions F 1 ðXÞ and F 2 ðXÞ are taken so that the coefficients of L * and I * are zero. For these cases, we get and Then, Equation (67) becomes 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 Ω = fðS, L, I, RÞ: S = S * , L/L * = I/I * = R/R * g. 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 fE 1 g. For each solution that intersects, ℝ 4 +0 fL = I = R 10 Computational and Mathematical Methods in Medicine = 0g limits to E 1 , which concludes that the disease persistence equilibrium E * of (1) is globally asymptotically stable in ℝ 4 +0 fL = I = R = 0g [24].

Lemma 7.
The infectious persistence equilibrium ðE * Þ of the model (1) is globally asymptotically stable when R 0 > 1, and the disease persists in the populations for a long time.

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 ðS 0 , L 0 , I 0 , R 0 Þ = ð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λ = μS 0 = 20; the rate which leavesLðtÞforIðtÞ, i.e., incubation periodδ = 1/7; and the fraction ofSðtÞwhich moves toLðtÞ;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 ðS 0 , L 0 , I 0 , R 0 Þ = ð 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 SðtÞ, LðtÞ, and IðtÞ, 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 ðIð43Þ = 3:126 × 10 3 Þ; after that, it is decaying.

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 R 0 > 1, and disease extinction steady state, where R 0 < 1. The outcomes obtained for both equilibria with stability analysis are also numerically demonstrated using MATLAB R2018a. Using the parameter values from Table 1, the basic reproductive ratio becomes R 0 = 5:86 > 1 thereby signifying the asymptotic stability of the infected steady state. For this reason, different initial conditions of ðS 0 , L 0 , I 0 , R 0 Þ are chosen as IC1 = ð10000, 0, 1, 0Þ, IC2 = ð8000, 0, 2, 0Þ, and IC3 = ð12000, 0, 3, 0Þ. Figure 4 illustrates the system dynamics of the susceptible, latent, and infected population for the three initial con-ditions 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 5 illustrates the system's phase portrait for different initial conditions. It represents the relative change of the susceptible SðtÞ, latent LðtÞ, and infected populations IðtÞ 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 R 0 = 5:86 > 1. In this case, the trajectories approach the long-term steady state, and the disease persists in the populations for t ⟶ ∞.
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 R 0 = 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 11 Computational and Mathematical Methods in Medicine different initial conditions. The trajectories are approaching a single point defined as infection-free steady-state E 0 = ð 10000, 0, 0, 0Þ and remain at this point for t ⟶ ∞.

Conclusion
We have proposed a diphtheria epidemic model and found two steady-state equilibrium points: one is disease extinction equilibrium point E 0 (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 R 0 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 R 0 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][42][43][44]. The enumerated infectious persistence equilibrium E * is ð1706,109,100,7814Þ and infection-free steady-state E 0 is ð10000, 0, 0, 0Þ. A statistical model was used to calculate the numeric value of the basic reproductive ratio R 0 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 R 0 = 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.

Data Availability
The data supporting the study are accessed from the study "R.

Conflicts of Interest
The authors declare that they have no conflicts of interest.