Abstract

Tuberculosis (TB) is a serious global health threat that is caused by the bacterium Mycobacterium tuberculosis, is extremely infectious, and has a high mortality rate. In this paper, we developed a model of TB infection to predict the impact of saturated recovery. The existence of equilibrium and its stability has been investigated based on the effective reproduction number . The model displays interesting dynamics, including backward bifurcation and Hopf bifurcation, which further results in the stable disease-free and stable endemic equilibria to be coexisting. Bifurcation analysis demonstrates that the saturation parameter is accountable for the phenomenon of backward bifurcation. We derive a condition that guarantees that the model is globally asymptotically stable using the Lyapunov function approach to global stability. The numerical simulation also reveals that the extent of saturation of TB infection is the mechanism that is fuelling TB disease in the population.

1. Introduction

Modeling and simulation are significant decision tools for controlling human diseases [1, 2]. Although every disease demonstrates its own series of biological characteristics, the models should be adjusted to every case in order to tackle real-world scenarios [3, 4].

Tuberculosis (TB) kills more people every year than any other infectious disease, including HIV and malaria, making it one of the most serious global health challenges [5]. Despite the fact that TB is a curable and preventable disease, it is still one of the top ten causes of mortality worldwide [5, 6]. TB is usually caused by the bacterium Mycobacterium tuberculosis. In 2018, there were 10 million new cases of TB reported, with 8.6 of those living with HIV [5]. To be infected with TB, a person just has to inhale a few TB bacteria [5, 6]. The infection may also affect other areas of the body, like the kidneys, brain, skin, and spine, in some cases. The TB infection is not new to the community, yet it has a long history of presence since antiquity in China, Egypt, and India [7]. Tuberculosis is transferred from one infected individual to another by germ-infested air droplets. When a person with active tuberculosis germs spits, talks, coughs, or sneezes, the TB germs are propelled into the air [5]. TB can be contracted effectively from relatives and companions. An individual becomes contaminated effectively when they inhale a couple of TB germs. Approximately ten percent of TB-infected people develop active TB infection, while the other ninety percent remain latent [5]. TB infection is not transmitted by people who are latently sick. Individuals with TB may be identified by skin or blood testing. People with immunocompromised diseases are more susceptible to TB (HIV and diabetic patients). A cough that lasts longer than three weeks, fever, coughing up blood, chest pain, fatigue, weight loss, exhaustion, and night sweats are all symptoms of active tuberculosis [5, 8]. TB is a transferable disease, which can be treated with medication therapy [5, 8].

Treatment is a significant and viable technique that can be used within a population to control the spread of infectious diseases. In modeling techniques, the role of treatment function addresses the likelihood of treatment against the disease at a given time for each infected person. In standard epidemic models, the treatment rate is proportional to the infective’s population size; however, in general, the rate of recovery is determined by medical resources such as medications, vaccinations, hospital beds, isolation areas, and the effectiveness of the treatments. It is understood in practice that societies have limited medical services for the treatment of infections, so offering these services places considerable pressure on facilities for public health. It is therefore more realistic to consider treatment functions of the saturated form, which appear to be large as the number of infectives increases [911]. Wang and Ruan [12] proposed an SIR model with a constant treatment as follows:that modeled a treatment capacity limit. The following piecewise linear treatment function was considered by Wang and Ruan [12]:where is the infectious threshold at which the healthcare system gains capacity; that is, treatment rises linearly with until capacity is achieved, and then it takes its maximal value, . This appears to be more realistic than the standard linear function. Eckalbar and Eckalbar [13] built the following SIR model using a quadratic treatment function:

Apart from that, we are aware that delaying treatment for infective persons has a negative impact on treatment efficiency. Zhang and Liu [9] employed a saturated treatment function that was both continuous and differentiable, , where corresponds to the cure rate, and the quantity represents the delay in treatment.

It has been noted that saturation in recovery or treatment function allows multiple endemic equilibria to occur when changes [9, 10, 1417]. For , the presence of backward bifurcation (the existence of endemic equilibrium along with disease-free equilibrium) implies that only is insufficient to eradicate the disease from the community [9, 14, 16, 17].In this situation, an endemic equilibrium is established along with the stable disease-free equilibrium (bistability of equilibria when ). Moreover, some researchers have observed the presence of limit cycles alongside the stable disease-free equilibrium when [16, 17]. Subsequently, if , the occurrence of a stable endemic equilibrium and oscillations in the infective population indicate that maximum efforts are needed to eliminate or monitor the disease [14, 16, 17].

Moreover, when , some different scenarios, for example, Hopf bifurcation and the occurrence of several endemic equilibria, also have been observed as a result of saturated recovery or treatment [12, 14, 15, 17]. For example, Cui et al. [14] proposed and investigated an SIS model with saturated recovery rate and discovered the occurrence of oscillatory behavior in the population via Hopf bifurcation because of the treatment capacity limitations. Wang and Ruan [12] presented and analyzed an SIR model with a constant treatment rate of infective. They noticed that when , the system experiences different bifurcations including saddle-node bifurcation and Hopf bifurcation. Global stability is observed without removal rate, while oscillatory behavior for constant removal rate was discovered. The behaviour of SEIR epidemic model with saturated treatment function is investigated by [18], a backward bifurcation that results in bistability emergence, and numerical result work suggested that they should improve the medical conditions to control the epidemic [19]. To fully understand the influence of delayed treatment on transmission of the dynamics of disease, Zhang et al. [20] considered an SEIR model approach with saturated incidence and saturated treatment method. Their finding recommends that to eliminate the disease, they should raise the effectiveness and extend the capacity of treatment. In other words, our medical technology and contribution of more drugs and beds provide timely treatment to patients.

In this paper, we incorporated saturated recovery which was excluded in previous four-dimensional -like models with exogenous infection used to study the phenomenon of backward bifurcation (see [21, 22]).

2. Construction of the Model

We examined the transmission dynamics of TB by employing a nonlinear ODE system of type system, which was developed in [21, 22]. Our aim is to better understand the interplay between exogenous reinfection and saturated treatment (recovery) rate in determining the outbreaks of TB. The whole population is classified into four (9) compartments, , , , and , which, respectively, denote susceptible, exposed, infected, and treated individuals. A normal mass action known as density-dependent incidence form interaction for the TB epidemic is examined and the population is considered to also be homogeneously mixed. The process of TB spreading is shown in Figure 1.

2.1. The Model with Saturated Recovery Rate

The nonlinear model that is considered in this study consists of a system of ordinary differential equations :with the initial conditions given by the following equation:

In this epidemiological system, the susceptible compartment is increased by recruiting individuals, either by immigration or birth, into the population at the constant rate . The term is the natural death rate. The exposed compartment becomes infectious at the rate of E and progresses to actively infected state. Exogenous reinfection can enable individuals that are already infected to develop active TB. This happens when they contract new infection from infectious individuals at a constant rate . Infected individuals are decreased in number because of treatment, which is given at the rate , . respectively. The term is the saturation form of treatment (recovery) to the infected populations with the recovery rate , and the parameter determines the magnitude of the impact of delaying the treatment of infected individuals.

3. Basic Properties of the TB Model with Saturated Recovery

Following the technique in [23, 24] (and furthermore the proof in [24]), it is not hard to demonstrate that when all model parameters remain nonnegative, the state variables are all positive for all time t given they are initially positive.

Theorem 1. The regionis positively invariant and attracts all nonnegative solutions of model system (4).

Proof. As mentioned in the study of [23], considering the nonlinear system of model (4), we take the first equation and let .that is,Integrating (9) by separation of variables givesand thereforeSimilarly, it can also been shown that , , , for all .

3.1. Invariant Region

Theorem 2. The solution of TB model system (4) is enclosed in the region subset , given byfor the initial conditions (5) in .

4. Analysis of Disease-Free Equilibrium , , and Effective Basic Reproduction Number

The disease-free equilibrium state, , of system model (4) is obtained by setting the right-hand side of (4) to zero which is given by

The basic reproduction number of system model (4) was calculated using the next generation operator approach (see [25, 26]), where it is obtained to bewhere is the effective reproduction number of system model (4). The following outcomes would guaranty the stability of , when the threshold quantity ascertains that the system will be free from TB infection under a threshold condition. A significant result is the following theorem.

Theorem 3. The disease-free equilibrium state, , of model system (4) is locally asymptotically stable (LAS) when and unstable if .

Proof. Consider model system (4), and the Jacobian of system (4) at is given byThen, the root characteristic equation of is , , , and . If , all the characteristic roots are negative and thus is locally asymptotically stable and unstable when .

4.1. Global Stability of State

In order to ensure that TB is independent of the initial size of the subpopulations of model system (4), it is necessary to show that is globally asymptotically stable (GAS) [2730]. Several approaches, such as the Lyapunov theorem, [27, 29, 31, 32] and global stability theorem have been used to prove the global stability of the disease free equilibrium, . Here, we shall apply the method that was introduced in [32].

Theorem 4. The disease-free equilibrium of model system (4) is global asymptotically stable in if and conditions and are satisfied as given in [32].

Proof. The two conditions and in Castillo-Chavez global stability theorem must be satisfied for to ascertain the global stability of the . Write system model (4) in the following form:Here, the components and where denotes uninfected population and denotes the infected population. The disease-free equilibrium is defined bywhereFor condition (the global asymptotical stability of ) to be fulfilled, we havewritten in the form of linear differential equations as follows:Solving linear differential equation (20) yieldsNow, clearly and as , regardless of the value of and .
Hence, is globally asymptotically stable.
Next, applying second condition of the theorem proposed in [32], that is, , , we haveThis is certainly an M- matrix (the off-diagonal elements of A are nonnegative) [32].
Sincethen,so thatSince , we have showed that condition H2 of the theorem proposed in [32, 33] is satisfied by our model. This shows that regardless of the initial population of infected, TB still can be controlled.

4.2. Endemic Equilibria and Bifurcation Analysis

To find endemic equilibrium as before, we set the right-hand side of system model (4) to zero. Solving for , , and from first, third, and last equation of (4), respectively, we havewhere and

By substituting and simplifying equations (26)–(28), we get the following equation in :where

Clearly, . When , then and thus . Furthermore, as as , . Accordingly, from the continuity of , at least one positive exists, such that , and thus there will be at least one endemic equilibrium of the model system (4). Notice that if and are positive, then model system (4) has a unique endemic equilibrium defined as from Descartes’ rule of sign. The number of possible positive real roots of the cubic polynomial (29) depends on the signs of , and . We, therefore, propose the following results.

Theorem 5. The TB model system (4) has(1)A unique endemic equilibrium when .(2)One or more than one endemic equilibrium when .(3)No endemic equilibrium when .

Now we will look at the polynomial in (29) to see how it influences the existence of endemic equilibrium.

If , the polynomial (29) reduces to a linear equation in , with a unique solutionwhich is positive if and only if . As a result, we can now say that if , then a unique endemic equilibrium exist if and there cannot be an endemic equilibrium when . To put it in another way, if the treatment of certain infected outpatients is not delayed, when the effective reproduction number is smaller than unity, an endemic steady state solution is not possible, ignoring the possibility of a backward bifurcation in this scenario.

4.3. Local Stability of Endemic Equilibrium

The Jacobian matrix of (4) at is given as

Here

We obtain the characteristic equation after some row and column operations of given bywhere

Thus, all the roots of the characteristic equation are with negative real parts if , , , and . Therefore, from Routh–Hurwitz criterion [34], the endemic equilibrium state is locally asymptotically stable if these conditions are true. The following result can be obtained.

Theorem 6. If and and are positive, then the unique endemic equilibrium of system model (4) is locally asymptotically stable provided .

4.4. Analysis of Backward Bifurcation

The existence and stability of endemic equilibrium are determined by investigating the possibility of backward or forward bifurcation due to the existence of endemic equilibrium. To investigate the possibility of backward or forward bifurcation of model system (4), we use the center manifold theory [35, 36]. Let the bifurcation parameter be .

Firstly, we obtained the bifurcation parameter at , and thusand therefore

To investigate the use of center manifold theory in [35], it is convenient to make simplification and transform the variables on model system (4). This is done by rewriting our system model (4).

Let

Moreover, by using the vector notation , model system (4) can be restated in the form of as follows:

We next show that the Jacobian matrix of (4) at the point has a simple zero eigenvalue, that is,

With , system model (39) has a simple zero real part eigenvalue, and all other eigenvalues are negative (i.e., has a hyperbolic equilibrium point). The center manifold theory [35] can therefore be used to examine the dynamics of the system (39) close to . It is practicable to derive the right eigenvectors of represented by , where

Conversely, the left eigenvectors represented as are easy to obtain, where

4.5. Computation of a and b

Hence, associated bifurcation coefficients are denoted by and , respectively; as explained in [35], when bifurcation coefficients and are both negative, then the system undergoes a backward bifurcation; otherwise, forward bifurcation will occur.

It is convenient to find both a and b defined by [35]

Taking into account model system (39) and examining a and b only, which are not equal to zero derivatives, for the termswe get

Now considering (43) and (44), we obtain

As indicated by the result shown in [35], our system experiences backward bifurcation at , only when both and are positive at . Obviously, is consistently positive. Hence, the positivity of offers the threshold circumstance for the phenomenon of backward bifurcation.

4.6. Analysis of Hopf Bifurcation

Hopf bifurcation of system (4) around endemic equilibrium point will be considered in this section. We employ the criterion that was derived in [37] to show that Hopf bifurcation exists. Let be an endemic equilibrium of system (4) when and be a bifurcation parameter for which is locally stable when and unstable when . At , loses its stability and periodic solutions emerge. The characteristic equation at endemic equilibrium equivalent to model system (4) at is given by

It was also observed that the component of the endemic equilibrium smoothly relies on , and the corresponding characteristic equation is given as

To this place, the coefficient is a smooth function of . Following the approach in [38], a Hopf bifurcation occurs when the following two conditions are satisfied:The Jacobian matrix has a pair of purely imaginary eigenvalues and the remaining eigenvalues have negative real parts..

According to [37], the coefficients of the characteristic equation have to meet the prerequisite for a pair of purely imaginary eigenvalues, which are, and ..

To complete the discussion, we need to drive the transversality condition (ii). For this reason, we let be a pair of purely imaginary eigenvalues corresponding to . Here, differentiating characteristic equation (48) with respect to , we get

Further, we get

The following relation will be used.

Here

If , then . We then summarize this result in the following equation:

Therefore, the transversality condition holds, and the system undergoes the phenomena of Hopf bifurcation at . Thus, transmission rate crosses its critical value, , and individuals start oscillating around endemic equilibrium point . This ends the proof of the theorem.

4.7. Global Stability of Endemic Equilibrium When

To prove the global stability of endemic equilibrium state of system model (4), we consider the case when exogenous reinfection parameter . We employed the Lyapunov function method and LaSalle’s invariance principles [39, 40]. System model (4) with therefore reduces towith

Furthermore, define

The following result is obtained.

Theorem 7. If , then the endemic equilibrium state given by (55) is in .

Proof. Consider the model without exogenous reinfection given by (55), and , so that the endemic equilibrium of the model (55) exists. In addition, consider the following Lyapunov function method of this type which has been utilized in the study of epidemiology and ecology, such as in [29, 4143], for model system (55) comprising the first three equations of (55), defined by the following. Define the Lyapunov functionand the derivatives of U are given byFilling in the derivatives and using (55) in , we haveAt the endemic steady state, it can be observed thatUsing expression (62) in equation (61) yieldsWe also observed from (55) at steady state thatUsing expression (64) in equation (63) givesWe can simplify the first term in the last equation of (65) as follows.Finally, it follows that since the arithmetic mean exceeds the geometric mean,Accordingly, it follows from this that when with if and only if , , and . Thus, is a Lyapunov function for system model (55) in . Thus, the largest compact invariant set where is a singleton , and therefore it follows that by [39], , , and as .

5. Numerical Simulation and Discussion

To confirm the feasibility of theoretical findings regarding stability and bifurcation of model system (4), we employed numerical continuation package XPPAUT [44] to perform one parameter bifurcation analysis and Matlab to perform numerical experiment. Motivated by some epidemiological literature of TB model [21, 45], model system (4) was validated by choosing a set of parameter values as illustrated in TB relevant literature, and the estimated parameters used are captured in each respective figure.

Figure 2 demonstrates the population of exposed individuals as the transmission rate is varied from 0 to 2. The y-axis defines the exposed individuals and x-axis defines the transmission rate beta . The red and black lines represent stable and unstable steady states, respectively. The green (blue, respectively) cycles illustrated the occurrence of stable and unstable limit cycles in this epidemiological system. In general, this epidemiological system demonstrates rich bifurcation dynamics which are epidemiologically important. Several dynamical behaviors of the system are observed as changes. As demonstrated in Figure 2, the system is identified by , which corresponded to Hopf bifurcation, transcritical bifurcation, saddle-node bifurcation, and saddle-node bifurcation of cycles. van Voorn [46] stated that “when an unstable limit cycle is born, while the equilibrium becomes stable, we have a subcritical Hopf bifurcation.” Hence, the bifurcation is subcritical. Indeed, Figure 2 reveals that for the bifurcation point, the stable equilibrium (red line) becomes unstable (black line). In addition, stable and unstable limit cycles also emerge. As the bifurcation parameter increases, the unstable limit cycles collide with stable limit cycles via a saddle-node bifurcation of cycles .

Figure 3 depicts a codimension-two bifurcation graph as and are changed. To examine the interaction of subcritical Hopf, saddle-node, and transcritical bifurcations, we conducted a codimension-two bifurcation analysis of TB model (4) as we varied the transmission rate and recovery parameter , and the result can be seen in Figure 3. There is indeed a codimension-2 point, black curve with red point, referring to cycle saddle-node bifurcation. This kind of bifurcation frequently occurs as a result of a collision between stable limit cycles (green curve) and unstable limit cycles (blue curve). In general, the codimension-2 bifurcation point serves as an organizing center, separating the parameter space into various curves with distinct qualitative findings: red curve corresponds to saddle-node bifurcation, blue curve corresponds to Hopf bifurcation, and cyan curve corresponds to transcritical bifurcation as demonstrated in previous result (that is, Figure 2, which illustrates the same qualitative mechanisms as Figure 3).

In Figures 4 and 5, we set the parameter values as , , , , , , and , and is varied around the critical value . The threshold quantity , and the eigenvalues are (−0.2163698605 + 1.349717265I, −3.652541848, −0.2163698605–1.349717265I) and the model system endemic equilibrium is given as (2.53865, 3.9028, 0.820965). Here, is found to be unstable. The endemic equilibrium is stable for and unstable for . We discover that the condition for Hopf bifurcation is fulfilled resulting in the presence of periodic oscillations around . The graphs for the results obtained are given in Figures 4 and 5.

Furthermore, we consider another set of representative parameters as , , , , , , , and . The threshold quantity , and the equilibrium is (2.11674, 0.997134, 1.0057) while the eigenvalues are given as (−0.231843, −0.231843 5.08514I). We can conclude that is stable since all of the eigenvalues are negative or negative real parts (see Figure 6).

Figure 7 depicts the phase diagram of and planes with different values of saturated recovery and with , , , , b = 0.1, and . When , the endemic equilibrium is unstable, and as we increase the value of to 0.1, the endemic equilibrium becomes stable; furthermore, we observed the appearance of two (5) stable limit cycles. We observed that the stability endemic equilibrium state changes from stable to unstable, and this kind of bifurcation is called subcritical Hopf bifurcation (see [46, 47]).

Figure 8 illustrates a typical bifurcation graphs of system model (4). In this diagram, the recruitment rate is changed, whereas other parameters are kept constant. The set of representative parameter values which result in Figure 8 are , , , , , , , and . Figure 8(a) demonstrates the scenario of forward bifurcation that when , the TB-free equilibrium is globally asymptotically stable, while if , the TB infection can persist. Conversely, as we observed from Figure 8(b), increasing the value of the parameter from 0.1 to 0.25, the TB infection can persist once defined for the range of values below unity which suggests the occurrence of the phenomenon of backward bifurcation. This suggests that diminishing below one will not be fundamentally enough for the eradication of TB epidemic from the population. When is adequately reduced such that , the positive equilibrium does not exist anymore and TB epidemic will cease to develop and will gradually fall from its generally high endemic to the TB-free equilibrium. From Figure 8(b), we observed that when , there are stable endemic equilibrium, an unstable endemic equilibrium, and a stable TB-free equilibrium, whereas when , there is just one steady endemic equilibrium.

Figure 9 illustrates the backward bifurcation diagram of system model (4) as the effective basic reproduction number varies against the infected population . Figures 9(a)9(d) demonstrate that increasing the value of contributes to the extension of the region of bistability, while reducing the value of brings about contraction of the bistability region. The TB eradication quantity also known as critical reproduction number shifts from left to the right when decreases and vice versa. Here, high value of means insufficient treatment for a huge population of TB infections; therefore, we prefer a scenario where there will be consistent TB epidemics within the population despite the fact that .

6. Conclusions

Dynamical models of TB infection have been analyzed by numerous authors [21, 22, 4850]. However, they all proposed linear recovery function to investigate transmission dynamic behavior of their model. However, we discovered that recovery from TB infection depends upon numerous factors such as antibiotic treatment and the number of hospital beds. This leads to nonlinearity in the number of recoveries. The nonlinear recovery function rate principle was then introduced to reveal some insight into the eradication of TB disease by examining the effect of antibiotic treatment and the number of hospital beds.

The main goal of this study is to investigate the qualitative dynamics of the TB infection incorporating saturated recovery (treatment) of the form into the model proposed in [21, 22] which gives a more realistic model. On the other hand, complicated dynamic behaviors, including oscillations, can be caused by saturated recovery [9, 10, 1417].

Typically in epidemiology, the criterion for the persistence and extinction of an infection is vital. The effective reproduction number is computed and the conditions for the local stability of equilibria and the existence of backward bifurcation and Hopf bifurcation are determined. Under the condition of , the infection disappears completely from Theorem 3. But if , then according to Theorem 3, there will be endemic disease. In reality, the epidemic cannot be eliminated except that the basic reproduction number is reduced under a lower level to such an extent that (see Figure 8(b)). The mathematical analysis of the TB model system (1) has shown two important mathematical phenomena: bistability and periodicity (oscillatory), which have been observed in some other infectious diseases. In the bistability phenomenon, which usually comes with a conversation of backward bifurcation (see [21] and the references therein), the system exhibits multiple endemic equilibria despite . In such situation, a stable endemic equilibrium has been shown to compete with a stable disease-free equilibrium. From these results, it has been confirmed that decreasing cannot control the spread of TB infection in the population. The second phenomenon undergoes oscillatory behavior under some certain conditions. The occurrence of limit cycles supports the behavior stated in several studies on the dynamics of some contagious diseases such as whooping cough, measles, rubella, and so on [5154]. Numerical findings indicate that which is the saturation parameter is accountable for backward bifurcation. Inability to intervene before TB infection has accumulated in the population will prompt a circumstance where a TB epidemic exists despite the fact that the . Improving existing medical technologies and channeling adequate resources in medicines can essentially encourage early intervention by ensuring that individuals infected with TB get treatment immediately.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

We would like to acknowledge the funder from the Ministry of Higher Education Fundamental Research Grant Scheme (FRGS) – Grant code (FRGS/2021/STG06/USM/02/9) and Article Processing Charge (APC) Fund from Research Creativity and Management Office (RCMO), Universiti Sains Malaysia. We would also like to thank the School of Mathematical Sciences for providing research and computing facilities.