# Analysis of the Epidemic Biological Model of Tuberculosis (TB) via Numerical Schemes

Authors: S. Kanwal; M.K. Siddiqui; E. Bonyah; K. Sarwar; T.S. Shaikh; N. Ahmed
Journal: Complexity
 (2022)

---

## Abstract

Tuberculosis (TB) is caused by bacillusMycobacterium tuberculosis (MTB). In this study, a mathematical model of tuberculosis (TB) is analyzed. The numerical behaviour of the considered model is analyzed including basic reproduction number and stability. We applied three numerical techniques to this model, i.e., nonstandard finite difference (NSFD) scheme, Runge–Kutta method of order 4(RK-4), and forward Euler (FD) scheme. NSFD scheme preserves all the essential properties of the model. Acquired results corroborate that NSFD scheme converges for each step size. While the other two schemes failed to preserve some properties of the model such as positivity and convergence. A graphical comparison presented in this study confirms the numerical stability of the NSFD technique shown here is maintained over a large area.


---

## Body

## 1. Introduction

Tuberculosis (TB) is a chronic, bacterial infectious disease. TB is one of the oldest diseases. In 1882, Robert Koch discovered the tubercle bacillus andMycobacterium tuberculosis which is the contributing agent of tuberculosis. It attacks the lungs and other organs of body as well. Tuberculosis is an airborne infection. When TB active individuals (infectious) cough, sneeze, speak, or spit, the tubercle bacilli spread in the air. When susceptible people inhale that air, they might become infected. The people who are regularly in contact with TB active are at higher risk. The bacillus mycobacterium sets up in lungs and transmits to other organs of the body if immune system of that individual is not strong enough to suppress it.The people who are infected have 10 percent risk to be infectious (active TB). Most of the infected individuals stay latently infected (noninfectious) for whole life, and this latent period may be from months to centuries. However, this period may depend on the co-infectious diseases that individual contains. With co-infectious diseases immunity decrease, the risk of infected to be infectious increases. The presence of HIV enhances faster the risk towards the active TB stage.According to world health organization (WHO) report 2019, about 10.0 million fell ill with tuberculosis in 2018 [1]. Globally, an estimated 1.2 million people expired due to tuberculosis among HIV negative and about 251,000 among HIV positive, in 2018. Pakistan reports about 518,000 new TB cases per year, and among 30 high burden countries, Pakistan ranks 6th [2]. Annually, an estimated 44,000 people expire due to tuberculosis, in Pakistan. Four lakh sixty two thousand nine hundred and twenty cases were reported in Khyber Pakhtunkhwa, province of Pakistan, in the duration from 2002 to 2017, according to NTP [3]. Newly registered case statics will be analyzed in this research work. Now, there is no direct way to discover if MTB has been removed or not. In this situation, mathematical models help us to estimate the future trend.To discover the parameters for infectious diseases, mathematical models use all collected data and statics and execute mathematical operations on them. These parameters are used to investigate the effects of various control strategies. Mathematical modeling determines which controls should be implemented and which should not. Modeling also aids in predicting the spread or contraction of infectious diseases in the future. Scientific experts have developed a large number of mathematical models for many communicable diseases throughout history. For a wide range of infectious diseases and epidemics, precise mathematical models are developed and implemented. The system of differential equations is often used in mathematical models. Different approaches can be used to control different contagious diseases. Vaccination and other control methods, as well as treatment, are used to combat some infectious diseases. To address the models’ flaws, several numerical techniques are employed, and results are derived. Simulations corroborate these findings.Waaler et al. [4] contributed significantly in the epidemiology of tuberculosis. They separated the population into three classes and proposed a model according to the characteristics of tuberculosis communication. The system of nonlinear ordinary differential equation was studied and a model for tuberculosis was constructed [5]. Carlos Castillo-Chavez and Feng [6, 7] worked a lot in the field of tuberculosis epidemiology. They analyzed the dynamics of tuberculosis and presented different models for detailed observations and provided results. Castillo and Song (2004) [8] presented a thorough analysis of the work on tuberculosis dynamics. They gathered different dynamical models of tuberculosis and provided a theoretical structure.Heesterbeek et al. [9] worked to explore the basic reproduction number and the estimated number of resulting cases generated by a typical primary infected individual during its period of infection. In [10], analysis of a SEIR model was presented for infectious disease which includes exponential normal birth and death rate and deaths due to disease, so the size of population might be changed with time. In 1994, Mickens suggested that numerical methods are the only way to find out exact solution of differential equations.Anguelov and Lubuma [11] presented their contributions to the construction of the nonstandard finite difference method and its application. Garba et al. [12] considered the problem of creating qualitative consistent finite difference schemes; they approximate with the innovative continuous-time model. To accomplish their goal, they considered a deterministic continuous-time model for the dynamics of transmission for any disease, for the presence of deficient vaccine. The model was thoroughly analyzed to investigate the dynamic features of the model.Gurski [13] presented a general construction of simple NSFD schemes for simple nonlinear systems and for construction used approximation techniques of standard differential equations, for example, predictor-corrector method and artificial viscosity. Memarbashi et al. [14] presented NSFD technique with the help of Mickens discretization scheme. Authors proved that the proposed NSFD scheme conserves the equilibrium points of the related continuous system. They analyzed the qualitative properties of the proposed system, and it preserves stability of equilibrium points, positivity as well as Neimark–Sacker bifurcations. The results proved the dynamical consistency of the considered discretized SEI epidemic model through the continuous system. Ahmed et al. considered an SEIQV model of saturated incidence rate and presented its numerical modeling. They developed a NSFD technique to solve the continuous model [15]. Authors analyzed the convergence of developed NSFD technique as well as showed the unconditional stability of the NSFD scheme. They compared the results of the proposed scheme to the RK-4 finite difference scheme. Comparison showed that RK-4 scheme produces nonphysical oscillations, fails to conserve the positivity property, and diverges even for small step size, as well as it converges to false stable position. Results were verified of the NSFD technique by simulations. Fatima et al. [16] presented the mathematical modeling of computer virus dynamics (SLBQRS). Ahmed et al. presented two operator splitting nonstandard finite difference schemes (NSFD) for the reaction diffusion SEIR epidemic model to solve numerically [17]. In 2019, Ahmed et al. proposed a Brusselator reaction diffusion model and checked the stability of the model by Neumann criteria of stability. They examined the solution by presenting finite difference scheme (FD), forward Euler explicit FD scheme, and semi-implicit Crank–Nicolson FD scheme [18].Ahmed et al. considered a model of autocatalytic glycolysis and investigated the numerical solution. Unknown variables showed the concentration of chemical elements, so the considered model established the positive results in [19]. Authors presented three different numerical techniques: the nonstandard finite (NSFD) method, the Runge–Kutta method (RK-4) method, and the forward Euler (FD) finite difference method.Kim et al., in 2018, presented a fitted model of tuberculosis to analyze the dynamics of TB transmission. They fitted data of Philippines and examined various control strategies such as hidden case discoveries, presence of case, active cases discoveries, and distancing. The analysis suggests that most favorable control strategies can reduce the number of infectious cases with minimum cost of implementation [20]. For tuberculosis (TB) transmission, authors presented a fractional-order delay differential model that includes the impacts of endogenous reactivation and external reinfections. All through the local stability of the steady states and bifurcation studies, they analyzed the qualitative characteristics of the model. The framework incorporates a discrete time delay to account for the time it takes to evaluate the ailment [21]. In [22], authors presented a stochastic epidemic model with time delays for COVID-19 dynamics (SIAQR). In the model, there is only one global favorable response with expecting value one. They derived a generalized stochastic threshold as a requirement of permanence and availability of an ergodic stationary population.In [1], to solve a SIR epidemic model, the authors provided two nonstandard finite difference (NFSD) approaches. The proposed approaches have key qualities including positivity and boundedness, as well as conservation law preservation. Numerical comparisons show that our method outperforms other existing standard methods including the second-order Runge–Kutta (RK2) method, the Euler method, and several readymade MATLAB scripts in terms of accuracy. In [2], authors constructed a nonstandard finite difference scheme to solve numerically a mathematical model for obesity population dynamics. Numerical comparisons between the proposed nonstandard numerical scheme and Euler’s method reveal that the suggested nonstandard numerical scheme is more effective. The nonstandard difference scheme methodology is a suitable alternative for solving numerically varied mathematical models where important properties of the populations must be achieved in order to represent the real world, as demonstrated by numerical examples. The creation of an NSFD scheme that is compatible with the aspects of a continuous dynamic model, such as positivity and population conservation, is thus the novel aspect of this work. In contrast to the smaller regions of other common numerical approaches, all of the numerical simulations with varying parameter values imply that the NSFD scheme described here keeps numerical stability in wide areas.In this work, the mathematical model of tuberculosis is numerically analyzed. Three numerical techniques, nonstandard finite difference (NFSD), forward Euler (FD), and Runge–Kutta of 4th-order (RK-4) schemes, are used to analyze the tuberculosis model. Numerical schemes analyze the mathematical presentation of proposed model resolution, permanency property, stability analysis, and threshold criteria. As a result, this study is divided into the following sections. We begin with a generalized tuberculosis model and description of the suggested model in Section2. Followed by a mathematical analysis, we compute the basic reproduction number and demonstrate the solution’s positivity, as well as the derivation of equilibrium points (DFE and EE) and equilibrium point stability, existence, and uniqueness in Section 3. Before discussing the results, some simulations are shown in Section 4. We constructed the different numerical schemes for the mathematical model of tuberculosis and provided the convergence analysis of the NSFD scheme for the proposed model. In Section 5, we present the graphical comparison of these numerical schemes. The last remarks bring this study article to a conclusion.The goal of this research is to build and analyze numerical schemes whose equilibrium points correspond with the continuous system’s equilibrium points while preserving all of the continuous system’s key attributes. The motivation of this study is to create consistent positivity-preserving numerical algorithms for the continuous model of tuberculosis (TB).


## 2. Materials and Methods

In this section, we present the generalized mathematical model of tuberculosis and describe its parameters.


### 2.1. Mathematical Model

The mathematical model for tuberculosis (TB) is given by [23]. We have a system of equations (1)–(5):(1)dSdt=Δ−θSIN−ωS,(2)dLdt=θSIN−ω+ξL+1−ΨφT,(3)dIdt=ξL+ΨφT−ω+Υ+ψ1I,(4)dTdt=ΥI−ω+φ+ψ2+χT,(5)dRdt=χT−ωR,where all the coefficients are positive and Nt=St+Lt+It+Tt+Rt. Here, we have the variables.St: suceptible classLt: exposed classIt: infected classTt: under treatment classRt: recovered classAt timet, initially, we have(6)S0>0,L0≥0,I0≥0,T0≥0,R0≥0.


### 2.2. Description of Parameters

The model parameters are estimated and described below:Δ: the recruitment rate of susceptible populationω: the coefficient of the natural death rate of all epidemiological human classesθ: the coefficient of transmission of TB infection from susceptible to infectedξ: the progression rate from L to IΥ: the rate at which infected individuals are treatedφ: the rate at which treated individuals leave the class TΨ: the parameter of treatment failure1−ΨφT: when Ψ=0, it means that all the infected class under treatment become latent, and when Ψ=1, it means the treatment failure and treatment failure and treated class will remain as a latent class due to the remainder of Mycobacterium tuberculosis or infected class I due to the failure of treatment at the rate in the infected class φψ1: the disease mortality rate in the infected class Iψ2: mortality ratio of treated peopleχ: the coefficient of recovery of treated individuals


## 2.1. Mathematical Model

The mathematical model for tuberculosis (TB) is given by [23]. We have a system of equations (1)–(5):(1)dSdt=Δ−θSIN−ωS,(2)dLdt=θSIN−ω+ξL+1−ΨφT,(3)dIdt=ξL+ΨφT−ω+Υ+ψ1I,(4)dTdt=ΥI−ω+φ+ψ2+χT,(5)dRdt=χT−ωR,where all the coefficients are positive and Nt=St+Lt+It+Tt+Rt. Here, we have the variables.St: suceptible classLt: exposed classIt: infected classTt: under treatment classRt: recovered classAt timet, initially, we have(6)S0>0,L0≥0,I0≥0,T0≥0,R0≥0.


## 2.2. Description of Parameters

The model parameters are estimated and described below:Δ: the recruitment rate of susceptible populationω: the coefficient of the natural death rate of all epidemiological human classesθ: the coefficient of transmission of TB infection from susceptible to infectedξ: the progression rate from L to IΥ: the rate at which infected individuals are treatedφ: the rate at which treated individuals leave the class TΨ: the parameter of treatment failure1−ΨφT: when Ψ=0, it means that all the infected class under treatment become latent, and when Ψ=1, it means the treatment failure and treatment failure and treated class will remain as a latent class due to the remainder of Mycobacterium tuberculosis or infected class I due to the failure of treatment at the rate in the infected class φψ1: the disease mortality rate in the infected class Iψ2: mortality ratio of treated peopleχ: the coefficient of recovery of treated individuals


## 3. Analysis of the Model

The present section contains the mathematical analysis of basic properties of the epidemic tuberculosis model: threshold quantity, equilibria, and stability.


### 3.1. Basic Reproduction Number

Basic reproduction number is of central importance and usually denoted byℜ0. It represents the number of resulting infections affected by primary when main infectious individuals are presented to host inhabitants where every individual is susceptible. It is the threshold quantity that predicts whether the infection will expire or not in a population. If ℜo>1, then infection will remain in population, and if ℜ0<1, then infection will vanish from the population. We calculated the basic reproduction number using the next generation method. Let x=S,L,I,T,Rτ; then, we can write the system of equations (1)–(5) as(7)dxdt=Fx−νx,where(8)Fx=θSIN00,and(9)νx=d1L−1−ΨφTd2I−ξL−ΨφTd3T−ΥI.Here,Fx and νx are the Jacobian matrices at DFE. Hence, we have(10)ℜ0=θξd3d1d2d3−d1ΥΨφ−1−ΨφξΥ.This is the spectral radius of next generation matrixFν−1.


### 3.2. Existence of Equilibria

Equilibrium points of the considered epidemic model exist and are as follows.Disease-free equilibrium: disease-free equilibrium is a state when disease dies out of the systemF0=S0,0,0,0,0, where S0=Δ/ω is a number of susceptible at DFE.Endemic equilibrium: endemic equilibrium is a state when the disease remains in the systemF1=S∗,L∗,I∗,T∗,R∗, where(11)S∗=N∗ℜ0,L∗=ω+ψ1+Υω+ψ2+χ+Υφ1−ΨN∗ℜ0−1ℜ0ω+ψ1+Υω+ψ2+χ+Υφ1−Ψ+ξΥ+d3,I∗=ξd3N∗ℜ0−1ℜ0ω+ψ1+Υω+ψ2+χ+Υφ1−Ψ+ξΥ+d3,T∗=ξΥN∗ℜ0−1ℜ0ω+ψ1+Υω+ψ2+χ+Υφ1−Ψ+ξΥ+d3,R∗=χT∗ω.There exist a unique endemic equilibrium, forℜ0>1.


### 3.3. Stability of the Model

#### 3.3.1. Local Stability at DFE

For disease-free equilibrium,F0=S0,0,0,0,0, the following is the variational matrix of the system of equations (1)–(5):(12)Z=z11z12z13z14z15z21z22z23z24z25z31z32z33z34z35z41z42z43z44z45z51z52z53z54z55,where(13)z11=−ω,z12=0,z13=−θ,z14=0,z15=0,z21=0,z22=−d1,z23=θ,z24=1−Ψφ,z25=0,z31=0,z32=ξ,z33=−d2,z34=Ψφ,z35=0,z41=0,z42=0,z43=Υ,z44=−d3,z45=0,z51=0,z52=0,z53=0,z54=χ,z55=−ω.Thus, we have(14)Z=−ω0−θ000−d1θ1−Ψφ00ξ−d2Ψφ000Υ−d30000χ−ω.By finding the eigenvalues of the Jacobian matrix, disease-free equilibrium is locally asymptotically stable, forℜ0<1. When ℜ0>1, the system becomes unstable at S0.


#### 3.3.2. Local Stability at EE

(15)E=−θI∗N∗−ω0−θS∗N∗00θI∗N∗−d1θS∗N∗1−Ψφ00ξ−d2Ψφ000Υ−d30000χ−ω.We have the characteristic equation of the Jacobian matrix of the form(16)λ4+λ3c1+λ2c2+λc3+c4=0,where(17)c1=θI∗N∗+ω+d1+d2+d3,c2=θI∗N∗d1+d2+d3+ωd1+d2+d3+d1d2+d2d3+d1d3−ΥΨφ−θS∗N∗ξ,c3=θI∗N∗(d1d2+d2d3+d1d3+ωd1d2+d2d3+d1d3+d1d2d3−θI∗N∗ΥΨφ−ωΥΨφ+ωθS∗N∗⋅ξ−d1ΥΨφ+θS∗N∗ξd3+1−ΨφξΥ,c4=θI∗N∗d1d2d3+ωd1d2d3−θI∗N∗d1ΥΨφ+θI∗N∗1−ΨφξΥ−ωd1ΥΨφ+ωθS∗N∗ξk3+1−ΨωφξΥ.According to Routh–Hurwitz criteria, all the eigenvalues of the characteristic polynomial are negative and conditions are satisfied, i.e,ci>0, for i=1,2,3,4 and c1c2c3>c32+c12c4. Endemic equilibrium (EE) of the system of equations (1)–(5) is locally asymptotically stable for ℜ0>1.


## 3.1. Basic Reproduction Number

Basic reproduction number is of central importance and usually denoted byℜ0. It represents the number of resulting infections affected by primary when main infectious individuals are presented to host inhabitants where every individual is susceptible. It is the threshold quantity that predicts whether the infection will expire or not in a population. If ℜo>1, then infection will remain in population, and if ℜ0<1, then infection will vanish from the population. We calculated the basic reproduction number using the next generation method. Let x=S,L,I,T,Rτ; then, we can write the system of equations (1)–(5) as(7)dxdt=Fx−νx,where(8)Fx=θSIN00,and(9)νx=d1L−1−ΨφTd2I−ξL−ΨφTd3T−ΥI.Here,Fx and νx are the Jacobian matrices at DFE. Hence, we have(10)ℜ0=θξd3d1d2d3−d1ΥΨφ−1−ΨφξΥ.This is the spectral radius of next generation matrixFν−1.


## 3.2. Existence of Equilibria

Equilibrium points of the considered epidemic model exist and are as follows.Disease-free equilibrium: disease-free equilibrium is a state when disease dies out of the systemF0=S0,0,0,0,0, where S0=Δ/ω is a number of susceptible at DFE.Endemic equilibrium: endemic equilibrium is a state when the disease remains in the systemF1=S∗,L∗,I∗,T∗,R∗, where(11)S∗=N∗ℜ0,L∗=ω+ψ1+Υω+ψ2+χ+Υφ1−ΨN∗ℜ0−1ℜ0ω+ψ1+Υω+ψ2+χ+Υφ1−Ψ+ξΥ+d3,I∗=ξd3N∗ℜ0−1ℜ0ω+ψ1+Υω+ψ2+χ+Υφ1−Ψ+ξΥ+d3,T∗=ξΥN∗ℜ0−1ℜ0ω+ψ1+Υω+ψ2+χ+Υφ1−Ψ+ξΥ+d3,R∗=χT∗ω.There exist a unique endemic equilibrium, forℜ0>1.


## 3.3. Stability of the Model

### 3.3.1. Local Stability at DFE

For disease-free equilibrium,F0=S0,0,0,0,0, the following is the variational matrix of the system of equations (1)–(5):(12)Z=z11z12z13z14z15z21z22z23z24z25z31z32z33z34z35z41z42z43z44z45z51z52z53z54z55,where(13)z11=−ω,z12=0,z13=−θ,z14=0,z15=0,z21=0,z22=−d1,z23=θ,z24=1−Ψφ,z25=0,z31=0,z32=ξ,z33=−d2,z34=Ψφ,z35=0,z41=0,z42=0,z43=Υ,z44=−d3,z45=0,z51=0,z52=0,z53=0,z54=χ,z55=−ω.Thus, we have(14)Z=−ω0−θ000−d1θ1−Ψφ00ξ−d2Ψφ000Υ−d30000χ−ω.By finding the eigenvalues of the Jacobian matrix, disease-free equilibrium is locally asymptotically stable, forℜ0<1. When ℜ0>1, the system becomes unstable at S0.


### 3.3.2. Local Stability at EE

(15)E=−θI∗N∗−ω0−θS∗N∗00θI∗N∗−d1θS∗N∗1−Ψφ00ξ−d2Ψφ000Υ−d30000χ−ω.We have the characteristic equation of the Jacobian matrix of the form(16)λ4+λ3c1+λ2c2+λc3+c4=0,where(17)c1=θI∗N∗+ω+d1+d2+d3,c2=θI∗N∗d1+d2+d3+ωd1+d2+d3+d1d2+d2d3+d1d3−ΥΨφ−θS∗N∗ξ,c3=θI∗N∗(d1d2+d2d3+d1d3+ωd1d2+d2d3+d1d3+d1d2d3−θI∗N∗ΥΨφ−ωΥΨφ+ωθS∗N∗⋅ξ−d1ΥΨφ+θS∗N∗ξd3+1−ΨφξΥ,c4=θI∗N∗d1d2d3+ωd1d2d3−θI∗N∗d1ΥΨφ+θI∗N∗1−ΨφξΥ−ωd1ΥΨφ+ωθS∗N∗ξk3+1−ΨωφξΥ.According to Routh–Hurwitz criteria, all the eigenvalues of the characteristic polynomial are negative and conditions are satisfied, i.e,ci>0, for i=1,2,3,4 and c1c2c3>c32+c12c4. Endemic equilibrium (EE) of the system of equations (1)–(5) is locally asymptotically stable for ℜ0>1.


## 3.3.1. Local Stability at DFE

For disease-free equilibrium,F0=S0,0,0,0,0, the following is the variational matrix of the system of equations (1)–(5):(12)Z=z11z12z13z14z15z21z22z23z24z25z31z32z33z34z35z41z42z43z44z45z51z52z53z54z55,where(13)z11=−ω,z12=0,z13=−θ,z14=0,z15=0,z21=0,z22=−d1,z23=θ,z24=1−Ψφ,z25=0,z31=0,z32=ξ,z33=−d2,z34=Ψφ,z35=0,z41=0,z42=0,z43=Υ,z44=−d3,z45=0,z51=0,z52=0,z53=0,z54=χ,z55=−ω.Thus, we have(14)Z=−ω0−θ000−d1θ1−Ψφ00ξ−d2Ψφ000Υ−d30000χ−ω.By finding the eigenvalues of the Jacobian matrix, disease-free equilibrium is locally asymptotically stable, forℜ0<1. When ℜ0>1, the system becomes unstable at S0.


## 3.3.2. Local Stability at EE

(15)E=−θI∗N∗−ω0−θS∗N∗00θI∗N∗−d1θS∗N∗1−Ψφ00ξ−d2Ψφ000Υ−d30000χ−ω.We have the characteristic equation of the Jacobian matrix of the form(16)λ4+λ3c1+λ2c2+λc3+c4=0,where(17)c1=θI∗N∗+ω+d1+d2+d3,c2=θI∗N∗d1+d2+d3+ωd1+d2+d3+d1d2+d2d3+d1d3−ΥΨφ−θS∗N∗ξ,c3=θI∗N∗(d1d2+d2d3+d1d3+ωd1d2+d2d3+d1d3+d1d2d3−θI∗N∗ΥΨφ−ωΥΨφ+ωθS∗N∗⋅ξ−d1ΥΨφ+θS∗N∗ξd3+1−ΨφξΥ,c4=θI∗N∗d1d2d3+ωd1d2d3−θI∗N∗d1ΥΨφ+θI∗N∗1−ΨφξΥ−ωd1ΥΨφ+ωθS∗N∗ξk3+1−ΨωφξΥ.According to Routh–Hurwitz criteria, all the eigenvalues of the characteristic polynomial are negative and conditions are satisfied, i.e,ci>0, for i=1,2,3,4 and c1c2c3>c32+c12c4. Endemic equilibrium (EE) of the system of equations (1)–(5) is locally asymptotically stable for ℜ0>1.


## 4. Numerical Schemes

In this section, three numerical schemes are developed according to systems (1)–(5), i.e., NSFD scheme, RK-4, and forward Euler (FD) scheme. NSFD scheme is used to solve systems (1)–(5) numerically and results will be compared with the other two schemes.


### 4.1. Runge–Kutta (RK-4) Scheme of Fourth Order

RK-4 is developed for the system of equations (1)–(5) as follows:Step 1:(18)u1=h−θSnInN−ωSn,v1=hθSnInN−ω+ξLn+1−ΨφTn,w1=hξLn+ΨφTn−ω+Υ+ψ1In,x1=hΥIn−ω+φ+ψ2+χTn,y1=hχTn−ωRn.Step 2:(19)u2=h−θNSn+u12In+w12−ωSn+u12,v2=hθNSn+u12In+w12−ω+ξLn+v12+1−ΨφTn+x12,w2=hξLn+v12+ΨφTn+x12−ω+Υ+ψ1In+w12,x2=hΥIn+w12−ω+φ+ψ2+χTn+x12,y2=hχTn+x12−ωRn+y12.Step 3:(20)u3=h−θNSn+d22In+w22−ωSn+u22,v3=hθNSn+u22In+w22−ω+ξLn+v22+1−ΨφTn+x22,w3=hξLn+v22+ΨφTn+x22−ω+Υ+ψ2In+w22,x3=hΥIn+w22−ω+φ+ψ2+χTn+x22,y3=hχTn+x22−ωRn+y22.Step 4:(21)u4=h−θNSn+u32In+w32−ωSn+u32,v4=hθNSn+u32In+w32−ω+ξLn+v32+1−ΨφTn+x32,w4=hξLn+v32+ΨφTn+x32−ω+Υ+ψ3In+w32,x4=hΥIn+w32−ω+φ+ψ3+χTn+x32,y4=hχUn+x32−ωRn+y32.Step 5:(22)Sn+1=Sn+16u1+2u2+2u3+u4,Ln+1=Ln+16v1+2v2+2v3+v4,In+1=In+16w1+2w2+2w3+w4,Tn+1=Tn+16x1+2x2+2x3+x4,Rn+1=Rn+16y1+2y2+2y3+y4.


### 4.2. Forward Euler (FD) Scheme

Forward Euler finite difference scheme developed for system of equations (1)–(5) is as follows:(23)Sn+1=Sn+h−θSnInN−ωSn,Ln+1=Ln+hθSnInN−ω+ξLn+1−ΨφTn,In+1=In+hξLn+ΨφTn−ω+Υ+ψ1In,Tn+1=Tn+hΥIn−ω+φ+ψ2+χTn,Rn+1=Rn+hχTn−ωRn.


### 4.3. Nonstandard Finite Difference Scheme

Nonstandard finite difference scheme was proposed by Mickens in 1989. NFSD scheme maintains dynamic consistency as well as numerical stability in terms of initial restriction with irregular step length. The system of equations (1)–(5) is developed according to the laws given by Mickens:(24)Sn+1−Snh=−θSn+1InN−ωSn+1,Sn+1=Sn+h−hθSn+1InN−hωSn+1,Sn+11+hθInN+hω=Sn+h,Sn+1=Sn+h1+hθIn/N+hω.Similarly, we have(25)Ln+1=Ln+hθSnIn/N+h1−ΨφTn1+hω+ξ,In+1=In+hξLn+hΨφUn(1+hω+Υ+ψ1,Tn+1=Tn+hΥIn1+hω+φ+ψ2+χ,Rn+1=Rn+hχTn1+hω.See more details in Tables1 and 2.Table 1
Disease-free equilibrium,R0<1.
Δ450,862.200 862 6θ0.243 3χ0.396 8Υ0.287 3ω0.014 78ψ10.220 2ψ20.055 0φ1.199 6Ψ0.150 0ξ0.200 7Table 2
Endemic equilibrium,R0>1.
Δ450,862.200 862 6θ0.543 3χ0.396 8Υ0.287 3ω0.014 78ψ10.220 2ψ20.055 0φ1.199 6Ψ0.150 0ξ0.200 7


### 4.4. Convergence Analysis of NSFD Scheme

Nonstandard finite difference scheme (NSFD) converges unconditionally at disease-free equilibrium (DFE) and endemic equilibrium (EE), see Appendix.


### 4.5. Positivity of NSFD Scheme

We have assumed that all the initial conditions are nonnegative.i.e., (26)S0≥0,L0≥0,I0≥0,T0≥0,R0≥0.These variables have approximate values which also are nonnegative due to the following supposition:Sn≥0,Ln≥0,In≥0,Tn≥0,andRn≥0.Solution of NSFD scheme’s formulas suggests the positivity of NSFD scheme,i.e.,Sn+1≥0,Ln+1≥0,In+1≥0,Tn+1≥0,andRn+1≥0.


## 4.1. Runge–Kutta (RK-4) Scheme of Fourth Order

RK-4 is developed for the system of equations (1)–(5) as follows:Step 1:(18)u1=h−θSnInN−ωSn,v1=hθSnInN−ω+ξLn+1−ΨφTn,w1=hξLn+ΨφTn−ω+Υ+ψ1In,x1=hΥIn−ω+φ+ψ2+χTn,y1=hχTn−ωRn.Step 2:(19)u2=h−θNSn+u12In+w12−ωSn+u12,v2=hθNSn+u12In+w12−ω+ξLn+v12+1−ΨφTn+x12,w2=hξLn+v12+ΨφTn+x12−ω+Υ+ψ1In+w12,x2=hΥIn+w12−ω+φ+ψ2+χTn+x12,y2=hχTn+x12−ωRn+y12.Step 3:(20)u3=h−θNSn+d22In+w22−ωSn+u22,v3=hθNSn+u22In+w22−ω+ξLn+v22+1−ΨφTn+x22,w3=hξLn+v22+ΨφTn+x22−ω+Υ+ψ2In+w22,x3=hΥIn+w22−ω+φ+ψ2+χTn+x22,y3=hχTn+x22−ωRn+y22.Step 4:(21)u4=h−θNSn+u32In+w32−ωSn+u32,v4=hθNSn+u32In+w32−ω+ξLn+v32+1−ΨφTn+x32,w4=hξLn+v32+ΨφTn+x32−ω+Υ+ψ3In+w32,x4=hΥIn+w32−ω+φ+ψ3+χTn+x32,y4=hχUn+x32−ωRn+y32.Step 5:(22)Sn+1=Sn+16u1+2u2+2u3+u4,Ln+1=Ln+16v1+2v2+2v3+v4,In+1=In+16w1+2w2+2w3+w4,Tn+1=Tn+16x1+2x2+2x3+x4,Rn+1=Rn+16y1+2y2+2y3+y4.


## 4.2. Forward Euler (FD) Scheme

Forward Euler finite difference scheme developed for system of equations (1)–(5) is as follows:(23)Sn+1=Sn+h−θSnInN−ωSn,Ln+1=Ln+hθSnInN−ω+ξLn+1−ΨφTn,In+1=In+hξLn+ΨφTn−ω+Υ+ψ1In,Tn+1=Tn+hΥIn−ω+φ+ψ2+χTn,Rn+1=Rn+hχTn−ωRn.


## 4.3. Nonstandard Finite Difference Scheme

Nonstandard finite difference scheme was proposed by Mickens in 1989. NFSD scheme maintains dynamic consistency as well as numerical stability in terms of initial restriction with irregular step length. The system of equations (1)–(5) is developed according to the laws given by Mickens:(24)Sn+1−Snh=−θSn+1InN−ωSn+1,Sn+1=Sn+h−hθSn+1InN−hωSn+1,Sn+11+hθInN+hω=Sn+h,Sn+1=Sn+h1+hθIn/N+hω.Similarly, we have(25)Ln+1=Ln+hθSnIn/N+h1−ΨφTn1+hω+ξ,In+1=In+hξLn+hΨφUn(1+hω+Υ+ψ1,Tn+1=Tn+hΥIn1+hω+φ+ψ2+χ,Rn+1=Rn+hχTn1+hω.See more details in Tables1 and 2.Table 1
Disease-free equilibrium,R0<1.
Δ450,862.200 862 6θ0.243 3χ0.396 8Υ0.287 3ω0.014 78ψ10.220 2ψ20.055 0φ1.199 6Ψ0.150 0ξ0.200 7Table 2
Endemic equilibrium,R0>1.
Δ450,862.200 862 6θ0.543 3χ0.396 8Υ0.287 3ω0.014 78ψ10.220 2ψ20.055 0φ1.199 6Ψ0.150 0ξ0.200 7


## 4.4. Convergence Analysis of NSFD Scheme

Nonstandard finite difference scheme (NSFD) converges unconditionally at disease-free equilibrium (DFE) and endemic equilibrium (EE), see Appendix.


## 4.5. Positivity of NSFD Scheme

We have assumed that all the initial conditions are nonnegative.i.e., (26)S0≥0,L0≥0,I0≥0,T0≥0,R0≥0.These variables have approximate values which also are nonnegative due to the following supposition:Sn≥0,Ln≥0,In≥0,Tn≥0,andRn≥0.Solution of NSFD scheme’s formulas suggests the positivity of NSFD scheme,i.e.,Sn+1≥0,Ln+1≥0,In+1≥0,Tn+1≥0,andRn+1≥0.


## 5. Graphical Comparison of the Epidemic Model Using Numerical Methods

The graphical comparison of the epidemic model using numerical methods is shown in Figures1–4.Figure 1
Simulations of forward Euler at disease-free equilibrium for step sizeh=1.2 and h=2. (a) Simulations of forward Euler. (b) Simulations of forward Euler.
(a)(b)Figure 2
Simulations of forward Euler at endemic equilibrium for step sizeh=1.2 and h=2 (a) Simulations of forward Euler. (b) Simulations of forward Euler.
(a)(b)Figure 3
Simulations of forward Euler at endemic equilibrium for step sizeh=1.2 and h=2. (a) Simulations of NSFD. (b) Simulations of NSFD. (c) Simulations of NSFD. (d) Simulations of NSFD.
(a)(b)(c)(d)Figure 4
Simulations of forward Euler at endemic equilibrium forh=1.2, h=2, and h=100. (a) Simulations of NSFD. (b) Simulations of NSFD (c) Simulations of NSFD.
(a)(b)(c)Figures1(a) and 1(b) show L (exposed), I (TB active), T (under treatment), and R (recovered) population using the Euler method at disease-free equilibrium for step size h=1.2 and h=2, respectively. Figures 1(a) and 1(b) show the behaviour of state variables obtained from forward Euler’s method at disease-free equilibrium. It is shown in Figure 1(a) that the forward Euler’s method gives the negative behaviour which is always meaningless in the population dynamics model. Also, by changing the step size in Figure 1(b), the forward Euler’s scheme gives the divergence behaviour.Similarly, Figures2(a) and 2(b) depict the unusual behaviour of the state variables evaluated from forward Euler’s method at endemic equilibrium point. On the same step sizes, Figures 2(a) and 2(b) give the negativity and divergence, respectively.Figures2(a) and 2(b) show L (exposed), I (TB active), T (under treatment), and R (recovered) population using the Euler Method at endemic equilibrium for step size h=1.2 and h=2, respectively.Figures3(a)–3(d) represent all population for the NSFD scheme for disease-free equilibrium. The NSFD method converges at different step sizes h=1.2, h=2, and h=100, respectively.Figures4(a)–4(c) represent the NSFD scheme for all population at endemic equilibrium. The NFSD scheme converges for even for large step sizes here the graphical representation for h=1.2, h=2, and h=100 respectively.


## 6. Conclusion

This study is based on the numerical study of the epidemic model of tuberculosis (TB). Basic reproduction number and stability analysis for disease-free equilibrium and endemic equilibrium points are evaluated. Three numerical schemes are applied to this model. The nonstandard finite difference (NSFD) scheme is being used to find out the numerical solution of the model. The NSFD scheme preserves essential properties, e.g., positivity, unconditional convergence, and stability. Here, we analyzed the convergence of nonstandard finite difference scheme. Results prove the unconditional convergence of the NSFD scheme of the epidemic model of tuberculosis. Results of forward Euler (FD) and Runge–Kutta (RK-4) are presented and compared. The NSFD scheme converges unconditionally for each step size, while forward Euler (FD) and Runge–Kutta (RK-4) fail to converge even at small step sizeh. Results are presented graphically and confirm the numerical stability of the NSFD technique shown here which is maintained over a large area, and a comparison is presented here. [24].



---

*Source: 5147951-2022-03-11.xml*