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 .

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 [36]. 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 [79]. 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 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 is indicated by that is partitioned into four groups: susceptible , latent (asymptotic) stage , individual affected by diphtheria in the acutely infected stage , and recovered individuals affected by diphtheria ; 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 from state at a rate of . Infectious people are successfully treated with a fixed rate , listing to the recovered state. Susceptible people acquire diphtheria infection among active diphtheria at rate , where signifies the infection transmission coefficient. Moreover,signifies a fraction of susceptible people that earn diphtheria infection and migrate to the latent diphtheria state, at rate, and the residual portion,, departs to the active diphtheria state. 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:

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 . Then, there exists continuous differentiable functions for positive time such that the 4-tuple covers (1) and .

Proof of Theorem 1. By Picard-Lindelöf theorem, it is narrated that the initial value problem has a unique solution for locally Lipschitz and continuous function in time , where . As the system (1) is autonomous, it is enough to prove that the function is locally Lipschitz in . Here, is defined as The Jacobian matrix of is obtained as This Jacobian is linear in . Thus, satisfies the continuity and differentiability for an interval . According to the mean value theorem, where . By assuming , we obtain for and thus, is bounded locally for each . Therefore, for all compact subset of , the derivative of is continuous and bounded and thus, is locally Lipschitz. Hence, according to the Picard-Lindelöf theorem, the initial value problem for has a unique solution .

Theorem 2. The proposed model (1) is invariant in the nonnegative orthant .

Proof. Let ; then, model (1) will take the form where and Here, and in matrix , all off-diagonal elements are greater than or equal zero. Hence, is a Metzler matrix and the system (1) is positive invariant in [32].

Theorem 3. For , any solution 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 The solution of Equation (12) is for all . Hence, the R.H.S. of Equation (13) is greater than or equal to zero, i.e., for all . 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 , and are greater than or equal zero .

Theorem 4 (boundedness). Suppose the model (1) satisfies and and has a unique solution on for some by Theorem 1; then, the state functions , and will be bounded and be positive .

Proof. Initially, suppose that the values of , and are positive. From Theorem 1, for , there exists a solution on . 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 . Therefore, we calculate each term on : instantly, the lower bounds on and can be placed. as the reduction expressions are linear; this achieves or or Applying initial condition, we get for
Again, as the reduction expressions are linear; this achieves Further, i.e., Similarly, by placing the upper bound on , we get i.e., where is an arbitrary constant which is depending on the upper bound of and . Now, by adding the equations forandand placing the bounds on this sum and by the positivity of these functions, for the upper bound of, we get where i.e., where the constant for that only depends on and . For the positivity of , and are positive, an upper bound can be placed on both , and by Now, can be bounded from below using where , i.e., Therefore, , and remain rigorously positive . Hence, there exists a for the continuity, at which the state variables and are still positive, which contradicts with the definition of and specifies that , and are rigorously positive on the whole interval . Moreover, all functions remain bounded with this interval; thus, the existing interval can be further extended. Actually, the bounds on and obtained earlier exist on each compact time interval. For the extension of the time interval to at which the solution endures and of the above discussion, the solutions continue both positive and bounded on .

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.,, or insists in the populationsasgrows large.

We start to determine the equilibria from the nonlinear intercommunicated terms into Equations (33), (34), and (35) that give

Thus, either or . Using in Equations (33), (34), and (35), we get the disease extinction equilibrium point as

By setting into Equations (32) and (35) yields the infectious persistence equilibrium that exists at the point

In the biological sense, 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 , 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 . 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 .

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

Therefore, we have,

Thus, the spectral radius of the matrixand the basic reproductive ratioare obtained [35].

Putting , we obtain,

This expression of 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 , we assume the following Lyapunov function:

By differentiation, we get

Substituting the values of , and in the above equation, we have

After substituting the value of , we are left with

At the disease extinction equilibrium , the basic reproductive ratio , and for all positive values of and , it is clear that . Hence, using LaSalle’s Invariance Principle [36], it is concluded that the model (1) is globally asymptotically stable.

Lemma 6. The infectious extinction equilibrium of the model (1) is globally asymptotically stable when , 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 , thus a Lyapunov function is assumed as where , and are all nonnegative constants to be obtained. This kind of Lyapunov function has been studied in [3740].

The infectious persistence equilibrium satisfies the following equations:

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

For the positive constants , and , the coefficients of , and must be zero, that is,

By solving, the above equation (55) yields

For advantage, we set up new variables , and to seek , and and setting the expressions of , and in Equation (54), we have

Multiplying by to the equation of (49) and the equation of (55) by yields

Hence, it follows that

Multiplying byto the last equation, whereis considered as a general function that will be determined later and, yields

Multiplying the equation of (49) by and the equation of (55) by yields

Hence, it follows that

Multiplying byto the last equation, whereis considered as a general function that will be determined later and, yields

From (54) using (63) and (66) yields

Now, the functions and are taken so that the coefficients of and are zero. For these cases, we get and

Then, Equation (67) becomes

By the arithmetic mean-geometric mean inequality, for equality, if and only ifand, the last expression must be less than or equal to zero. Thus, we have with equality if and only if and . By LaSalle’s Invariance Principle [36], for each solution, the omega-limit set remains in an invariant set that is contained in . Since must be in turns zero, which implies that , , and . Thus, there is only invariant set in which is singleton . For each solution that intersects, limits to , which concludes that the disease persistence equilibrium of (1) is globally asymptotically stable in [24].

Lemma 7. The infectious persistence equilibrium of the model (1) is globally asymptotically stable when , 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 . There are seven parameters in our model which are to be obtained. Among these parameters, natural mortality rateis estimated as 0.002; the recruitment of susceptible class; the rate which leavesfor, i.e., incubation period; and the fraction ofwhich moves to;; and disease-induced mortality rate is estimated as. 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 , and the initial condition of the state variables is . 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 , and , respectively. From these figures, it is observed that the infected population (-class) increases significantly upon the infection and arrives at the peak at the 36th day ; after that, it is decaying.

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 , and disease extinction steady state, where . 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 thereby signifying the asymptotic stability of the infected steady state. For this reason, different initial conditions of are chosen as , and .

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 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 , latent , and infected populations 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 when the basic reproduction number . In this case, the trajectories approach the long-term steady state, and the disease persists in the populations for .

For disease-extinction equilibrium, we assume the value of infection transmission rate different from Table 1 as . Therefore, the basic reproductive ratio is evaluated as . 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 and remain at this point for .

9. Conclusion

We have proposed a diphtheria epidemic model and found two steady-state equilibrium points: one is disease extinction equilibrium point (37), and another is infectious persistence equilibrium point (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 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 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, 4144]. The enumerated infectious persistence equilibrium is and infection-free steady-state is . A statistical model was used to calculate the numeric value of the basic reproductive ratio 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 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. Matsuyama, A. R. Akhmetzhanov, A. Endo, H. Lee, T. Yamaguchi, S. Tsuzuki, H. Nishiura, Uncertainty and sensitivity analysis of the basic reproductive number of the diphtheria: a case study of a Rohingya refugee camp in Bangladesh, November-December 2017, PeerJ, 6 (2018) 4582. doi:10.7717/peerj.4583.”


This work has been partially presented in my M.Sc. thesis work.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors like to express their gratitude to the Bangladesh University of Engineering and Technology (BUET) for providing financial support under the Basic Research Grant No. 1111202109017.