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.

1. 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-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 ℛ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.

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 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 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 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.

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: with following subsidiary conditions:

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 that S0, 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, R)(0) = (S0, L0, I0, R0).

Proof of Theorem 1.

By Picard-Lindelöf theorem, it is narrated that the initial value problem y′(t) = g(y(t)), y(t0) = y0 has a unique solution y(t) 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

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∗ ∈ I1. By assuming |∇g(y∗)| = M, we obtain |g(y1) − g(y2)| ≤ M|y1 − y2| for y1, y2 ∈ I1 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) = y0 for t0 > 0 has a unique solution y(t).

Theorem 2 .

The proposed model (1) is invariant in the nonnegative orthant ℝ+4.

Proof

Let Y = (S, L, I, R); then, model (1) will take the form where and

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 .

For t > 0, any solution (S, L, I, R) of the model (1) with condition (2) is positive.

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

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 respectively. Those solutions show that all L(t), I(t), and R(t) 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 S(t), L(t), I(t), and R(t) will be bounded and be positive ∀t ∈ [0, t0].

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 𝒯∗ 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 𝒯∗ > 0. Therefore, we calculate each term on [0, 𝒯∗]: instantly, the lower bounds on L, I, and R can be placed. as the reduction expressions are linear; this achieves or or

Applying initial condition, we get for t ∈ [0, 𝒯∗].

Again, 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 C1 ≥ max{βC, μ, (μ + α)}, i.e., where the constant C2 > 0 for t ∈ [0, 𝒯∗] that only depends on L(0), I(0), R(0), and C1. 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 C3 ≥ max{βC2, μ}, ⇒(dS/dt) + C3(1 + e)S ≥ 0, i.e.,

Therefore, S, L, I, and R remain rigorously positive ∀t ∈ [0, 𝒯∗]. Hence, there exists a 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 𝒯∗ 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].

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

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, 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 by ℛ0. 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 ℛ0.

Using the next-generation matrix method on the model (1), we get

Therefore, we have,

Thus, the spectral radius of the matrixℱ𝒱−1and the basic reproductive ratioℛ0are obtained [35].

Putting S0 = λ/μ, we obtain,

This expression of ℛ0 represents the basic reproductive ratio for the model (1).

Remark 1 .

The infectious equilibrium point with the expression of basic reproduction number ℛ0

6. Global Stability Analysis

6.1. Global Stability at Infectious Extinction Equilibrium

For disease extinction equilibrium E0 = (S0, L0, I0, R0) = (λ/μ, 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 S0 = λ/μ, we are left with

At the disease extinction equilibrium E0, 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 1 .

The infectious extinction equilibrium (E0) of the model (1) is globally asymptotically stable when ℛ0 ≤ 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 equilibrium E∗ = (S∗, L∗, I∗, R∗), thus a Lyapunov function is assumed as 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 equilibrium E∗ = (S∗, L∗, I∗, R∗) satisfies the following equations:

Now, differentiate U with respect to time t, which can be further simplified to

For the positive constants B1, B2, and B3, 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 B1, B2, and B3 in Equation (54), we have

Multiplying by B1 to the 2nd equation of (49) and the 3rd equation of (55) by L∗ yields

Hence, it follows that

Multiplying byF1(X)to the last equation, whereF1(X)is considered as a general function that will be determined later andX = (x, y, z, u), yields

Multiplying the 4th equation of (49) by B3 and the 4th equation of (55) by R∗ yields

Hence, it follows that

Multiplying byF2(X)to the last equation, whereF2(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 F1(X) and F2(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 Ω = {(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, ℝ+04{L = I = R = 0} limits to E1, which concludes that the disease persistence equilibrium E∗ of (1) is globally asymptotically stable in ℝ+04{L = I = R = 0} [24].

Lemma 2 .

The infectious persistence equilibrium (E∗) of the model (1) is globally asymptotically stable when ℛ0 > 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 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 (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 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 × 103); after that, it is decaying.

Table 1

Description and value of the parameters of the diphtheria model.

Parameter	Description	Value	Source
λ	The recruitment of susceptible class	200 persons day−1	Estimated
μ	Natural mortality rate	0.002 day−1	Estimated
α	Disease induced mortality rate	0.0054 day−1	Estimated
β	Disease transmission rate	0.000097 persons−1 day−1	Fitted
l	The fraction of S(t) which moves to L(t)	0.95	Estimated
γ	Recovered rate	0.156 day−1	Fitted
δ	The rate which leaves L(t) for I(t)	0.143 day−1	Estimated

Figure 2

The diphtheria model (1) simulation in log scale.

Figure 3

The fitted diphtheria model (1).

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 ℛ0 > 1, and disease extinction steady state, where ℛ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 ℛ0 = 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).

Figure 4 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 between S(t), L(t), and I(t) of diphtheria model (1) when ℛ0 = 5.86 > 1 for different initial conditions.

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 ℛ0 = 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 ℛ0 = 5.86 > 1 for different initial conditions.

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 ℛ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 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 between S(t), L(t), and I(t) of diphtheria model (1) when ℛ0 = 0.302 < 1 for different initial conditions.

Figure 7

System's phase of diphtheria model (1) in 3D when ℛ0 = 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 point E0 (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 ℛ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 ℛ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–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 ℛ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 ℛ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.