A novel tuberculosis model incorporating a Caputo fractional derivative and treatment effect via the homotopy perturbation method

Tuberculosis (TB), caused by Mycobacterium tuberculosis, is a contagious infectious disease that primarily targets the lungs but can also impact other critical systems such as the bones, joints, and neurological system. Despite significant efforts to combat TB, it remains a major global health concern. To address this challenge, this study aims to explore and evaluate various tuberculosis control approaches using a mathematical modeling framework. The study utilized a novel SEITR mathematical model to investigate the impact of treatment on physical limitations in tuberculosis. The model underwent qualitative analysis to validate key aspects, including positivity, existence, uniqueness, and boundedness. Disease-free and endemic equilibria were identified, and both local and global stability of the model was thoroughly examined using the derived reproduction number. To estimate the impact of each parameter on each compartment, sensitivity analysis was conducted, and numerical simulations were performed using Maple 18 software with the homotopy perturbation method. The obtained results are promising and highlight the potential of the proposed interventions to significantly reduce tuberculosis virus prevalence. The findings emphasize the significance of fractional-order analysis in understanding the effectiveness of treatment strategies for mitigating tuberculosis prevalence. The study suggests that the time fractional dynamics of TB treatment correspond to the treatment’s efficacy, as the conceptual results showed that non-local interactions between the disease and the treatment may lead to more accurate ways of eradicating tuberculosis in real-world scenarios. These insights contribute to a better understanding of effective treatment strategies and their potential impact on tuberculosis control and public health. In conclusion, scientists, researchers, and healthcare personnel are urged to take action and utilize the discoveries from this research to facilitate the eradication of the hazardous tuberculosis bacteria.


Background
Tuberculosis (TB) is a highly contagious disease caused by the Mycobacterium tuberculosis bacteria, primarily affecting densely populated areas and posing a significant public health threat.It spreads through airborne transmission when infected individuals cough or have close contact with others.Common symptoms include a persistent cough, chest pain, weight loss, and a fever.Diagnosis involves a comprehensive assessment of medical history and various examinations, such as X-rays and tests.TB is particularly prevalent in low-and middleincome countries with limited access to healthcare, impacting millions of people worldwide.An active TB patient can infect 5 to 15 susceptible individuals in their environment (WHO 2021).
To combat TB, several preventive measures have been proposed, including treatment interventions (Luju and Yan 2014;Ramli et al. 2019;Clark et al. 2019).(Luju and Yan 2014) incorporated treatment interruptions and latent periods into their approach, while (Ramli et al. 2019) emphasized the importance of multiple drug administration's within specific timeframes.Clark et al. (2019) conducted sensitivity analyses to identify effective strategies for treatment and prevention during epidemics.(Ullah et al. 2020a, b) employed the modified Adams-Bashforth technique to calculate numerical solutions for a fractional tuberculosis model, highlighting the significance of the proposed derivative in understanding TB dynamics and control.(Fatmawati et al. 2020) focused on investigating TB infection dynamics in children and adults using a novel fractional model and fractional calculus, particularly the Caputo and Atangana-Baleanu derivatives.Their results indicated a significant reduction in TB infection among individuals already affected by the disease when the fractional order was decreased.(Ullah et al. 2020a, b) explored optimal TB control techniques involving vaccines and treatment, emphasizing the importance of community-based TB control and showcasing the impact of disease control through simulations.Furthermore, (Ullah et al. 2018) proposed a fractionalorder model for simulating TB dynamics in Khyber Pakhtunkhwa, Pakistan, utilizing the Caputo derivative to demonstrate its superiority over classical models.
Recent research has focused on investigating the agespecific aspects of TB infection, with newborn vaccination being suggested as a preventive measure (Tilahun et al. 2020;Mengistu et al. 2020).The urgency to develop a more effective and longer-lasting TB vaccine has heightened due to the emergence of drug-resistant strains.However, surpassing the efficacy and durability of the current BCG vaccine has proven to be a significant challenge (Schrager et al. 2020).To address this challenge, researchers have developed the TB Vaccine Development Pathway web tool (Roordink et al. 2021), which serves as a valuable resource for guiding vaccine development from discovery to deployment, ultimately aiming to eradicate TB.
Various preventive strategies,including therapy interventions, infant immunization, and innovative vaccine development according to (Zhang et al. 2023;Kelemu Mengistu and Witbooi 2019), have been proposed.It is crucial to employ well-designed trials and models to plan for success and ensure the timely distribution of vaccines to those most in need.Different numerical methods, such as the homotopy perturbation method (Kolawole et al. 2023;Olayiwola et al. 2023) and the Laplace Adomian decomposition method (Yunus et al. 2023), have been employed to control various diseases, including COVID-19 and Lassa fever.Fractional calculus has also been utilized to develop mathematical models for epidemic diseases, including TB.The Caputo fractional derivative operator has been employed to describe TB transmission dynamics, taking into account memory behavior (Danane et al. 2020).Additionally, Zhang et al. (2021) delve further into the application of the Caputo operator for TB transmission and prevention, highlighting the crucial role of appropriate treatment in minimizing TB transmission and prevalence.Ahmad et al. (2021) consider a mathematical model for TB that incorporates the Caputo and Fabrizio derivatives with an exponential kernel.Furthermore, (Zafar et al. 2022) examine a newly discovered fractional operator (FO) with a ML non-singular kernel for different fractional orders, utilizing the fractional Caputo and Atangana-Baleanu-Caputo predictor and corrector approaches.Fractional-order derivatives in the context of epidemic disorders have been extensively studied in previous research (Liu et al. 2023;Farman et al. 2023;Nazir et al. 2020;Okyere and Ackora-Prah 2023;Ahmed et al. 2021;Hajaj and Odibat 2023) using operators such as Atangana-Baleanu-Caputo, Caputo-Fabrizio, and Caputo fractional derivatives.Several fractional-order TB models have investigated the impact of quarantine, hospitalization, and treatment, while (Ullah et al. 2023) took a unique deterministic epidemic approach, examining the influence of treatment adherence and awareness on TB dynamics.
Despite the wealth of research on fractional-order models for tuberculosis, the specific effects of treatment rate and effective contact rate on the disease have not been thoroughly explored.This research seeks to fill this gap and contribute valuable insights into how these factors influence tuberculosis dynamics.By analyzing the influence of treatment rate on the fractional order and the impact of the effective contact rate on tuberculosis transmission, this study aims to enhance our understanding of disease control strategies and pave the way for more effective and targeted interventions to combat tuberculosis.To achieve this, we modified the model presented by (Kereyu and Demie 2021) to suit the study of tuberculosis dynamics in Nigeria, specifically considering the incorporation of the treatment class, leading to a significant improvement in disease cure.
The paper's structure includes a detailed description of the model formulation and description in "Model formulation and description" section, the model analysis in "Model analysis" section, an exploration of the method in "Methods" section, the presentation of results in "Results" section, and a conclusion in "Discussion" section.

Model formulation and description
A five-compartmental model of tuberculosis is described, consisting of susceptible (S(t)), exposed E(t) , infected I(t) , treatment T (t) , and recovered R(t) compartments.Susceptible individuals are those who are vulnerable to contracting the disease if they come into contact with someone who is infected.Exposed individuals E(t) are the group of individuals who have come into contact with an infectious agent but have not yet developed symptoms of the disease.The infected class comprises confirmed cases of people who have developed the disease and are prone to transmitting it to others.Individuals undergoing treatment for any illness are at risk of developing tuberculosis, while those who have recovered are denoted as R(t) and receiving treatment are denoted as T(t).The parameter shows the recruitment rate which occurs in the susceptible class only, and susceptible individuals have the tuberculosis through contact with the infected individuals at the rate β N S(t)I(t) where β is the effective contact rate.δ represents the death due to TB, while η also rep- resents the treatment rate of infected class.τ indicates the progression rate of exposed to infected class, and σ represents natural cure rate.The re-infection rate is ϕ, and the individual leaves the treatment class at the rate θ .Based on the details mentioned earlier, a nonlinear SEITR model for tuberculosis is formulated, and this model is illustrated through the schematic diagram in Fig. 1. (1) Equation ( 1) can be rewritten as where K = β N .With the initial conditions: Table 1 shows the description of variables and parameters.

Preliminaries of fractional calculus
In this section of the paper, some preliminaries definitions are discussed based on fractional calculus.
Definition 1 A fractional order can be represented by κ ∈ (0, 1] using the Riemann-Liou- ville method of the function g (3) dx, and the following prop- erties are satisfied:

Definition 2
The Caputo derivative of a in the interval [0, T ] can be written as

Model analysis
In this section, we shall conduct the qualitative analysis of the mathematical model, and we shall do this with an integer order α = 1 to fully study the properties of the mathematical model and show its potential for real-life applications.

Positivity and boundedness
This section focuses on investigating the positivity and boundedness of solutions to the model in Eq. ( 2) with the initial condition described in Eq. ( 3).
Theorem 1 For t ≥ 0 , the system of Eq. (2) and the starting condition (3) have a positive solution.
Proof Assume that Eq. ( 2) has a solution S(t), E(t), I(t), T(t), and R(t) with the initial condition (3).So, we have i.
Therefore, solution is positive for all t ≥ 0 and hence bounded.

Dynamic behaviors of the model
The dynamic behavior of the model was investigated by obtaining the disease-free and endemic equilibrium of the tuberculosis model.

Disease-free equilibrium
A disease-free equilibrium (DFE) is a population that is entirely free of infectious diseases.Thus, when there is no infection I = 0 , the equilibrium points result in:E 1 = S 0 , E 0 , I 0 , T 0 , R 0 = � µ , 0, 0, 0, 0 .

Endemic equilibrium points
Endemic equilibrium is a steady state of an infectious disease in a population in which the prevalence of the disease remains stable over time.Hence, E = I = 0 and the follow- ing thresholds are obtained: .

Basic reproduction number
The reproductive number of the SEITR model associated with the reproductive power of the disease is defined by R 0 = ρ(G) where ρ is the spectral radius of the next-gener- ation matrix G = FV −1 .So the R 0 of the model Eq. ( 2) at diseases-free equilib- rium is obtained as: So, the next-generation matrix is resulted in Therefore, R 0 = β�τ N (µ+σ +τ )(η+µ+δ) which represent the dominant eigenvalue.

Local stability for disease-free equilibrium
At the equilibrium state, the equations of the model (2) can be written as follows where K = β N .At the disease-free equilibrium point, the Jacobian matrix of Eq. ( 9) is obtained as Also, the characteristics equation can be obtained using det (J (E 1 ) − I) to obtain the eigenvalues .
The characteristics equation derived from Eq. ( 11) is Using the Routh-Hurwitz criterion stability (Gantmacher 1959) of order four, a 3 a 2 a 1 > a 2 1 + a 2 3 a 0 where a i > 0.
From the above equations, it can be observed that a 3 > 0, a 2 > 0, a 1 > 0 and a 0 > 0.
So, the Routh-Hurwitz criterion is satisfied.Hence, the disease-free equilibrium of the system is asymptotically stable. (10)

Global stability of endemic equilibrium
Theorem 2 If R 0 > 1 , the endemic equilibrium point of the model Eq. ( 2) is globally asymptotically stable.
Proof The global stability is constructed using Lyapunov function as follows: Obtaining the derivative of Lyapunov of the solution of Eq. ( 2) results in: where we have Inserting ( 18) into (17), we arrive at Then, we obtain ( 16) , where P and Q can be extracted from Eq. ( 20) and then give the formula as dV dt = P − Q which is then written as: (21) and Therefore, if P < Q , the dV dt will be negative, also if Hence, the largest invariant set in (S(t), E(t), I(t), T (t), R(t)) : dV dt = 0 is the single point set E * .Using the LaSalle's invariance principle adopted by (LaSalle 1976), the endemic equilibrium is globally asymptotically stable if P < Q.

Homotopy perturbation method
Our study focuses on utilizing the homotopy perturbation method as a numerical approach to solve the mathematical model.The homotopy perturbation method is presented below, and we explore its effectiveness in accomplishing this task.
Firstly, consider the classical order procedure for solving coupled system of differential equation Subject to the boundary condition where indicates the operator of differential, represents the boundary condition operator, m(t) denotes the analytic function, and υ n is the normal vector derivative.The function �(υ) is divided into two which are: L H (υ) and N H (υ) denotes the linear and nonlinear operator in the equation.

So, (25) implies
The homotopy for Eq. ( 26) can be obtained as where p is an encoding parameter that can change from [0, 1] by a deformation process.Equation ( 27) is reduced further to yield: as p → 0, Equation (28) gives: And when p → 1, Equation ( 30) can be written in the form of power series and gives As a result of combining ( 31) and ( 28), and comparing coefficients of equal powers of p, we have the solution: (23) (32) The application of He's homotopy perturbation method (He 1999) is subsequently demonstrated for a system of ordinary differential equations of fractional order using the following algorithm: The fractional-order derivative of Eq. ( 2) can be written as: Using the HPM method as described in Olayiwola et al. (2023) and Abdulaziz et al. (2008), the theory of HPM can be applied to Eq. ( 33).Consider a system of ordinary differential equations of Caputo fractional order α defined as: Constructing a homotopy for Eq. ( 34) is written as follows: Simplifying Eq. ( 35) results in where p ∈ (0, 1) is an embedding parameter.At p = 0, Eq. ( 36) becomes linear such that the following equation is obtained: At p = 1, the original equation in ( 37) is obtained.Let's assume a solution series embedding the parameter p for (36) such that Substituting Eq. ( 38) into Eq.( 36) and comparing the coefficients of equal powers of p , the following series of equations is obtained: (33) (38) S i (t) = s i0 + ps i1 + p 2 s i2 + p 3 s i3 , . . .p n s in , and so on.Also, these systems of equation in ( 39) can be solved by applying the Riemann-Liouville fractional integral operator I α 1 to obtain the values of s i1 (t), s i2 (t), s i3 (t) . . . .Thus, the solution of ( 34) is obtained as:

Numerical experiment
This section detailed how the homotopy perturbation technique was used to run the numerical simulation that produced the approximate solution to the SEITR epidemic model.
The SEITR epidemic model's Caputo time fractionalorder derivative in (2) is provided as: dt α i for i = 1, 2, 3, 4, 5. Constructing a homotopy for (39), (39) The approximate solution of (41) can be assumed as: substituting ( 43) for ( 42) and trying to compare coefficients with identical p powers Solving Eq. ( 42) yields Similarly comparing the coefficients of p 1 , Applying the property of Riemann-Liouville integral (Liouville 1832) in Definition 1, and using solution (45), results in the system of equations, ( 43) (46) The coefficients of p 2 equally yield: By solving these equations, the second approximation is achieved as follows: . .

Results
Numerical simulations and analysis of ℜ 0 .

Experiment I Sensitivity analysis on R 0
We analyze the effect of some model parameter on R 0 using the formula presented below and parameter in Table 2 T 2 (t) =µ(−ηi 0 + (θ + µ)v 0 )  The sensitivity of parameters η changed in the infected group.This suggests that raising the treatment rate η leads to a fall in the proportion of infected persons seeking treatment, which decreases the R 0 .To reduce the virus transmission, we must lower the rate at which exposed individuals become infectious, denoted by σ .Additionally, a decrease in sigma also negatively impacts the basic reproduction number R 0 .Hence, the rate should be decreased.Lastly, prioritizing treatment plays a vital role in reducing TB virus transmission effectively.

Experiment II Numerical result analysis
To eradicate tuberculosis infection, we present the numerical results graphically by applying the homotopy perturbation method to the tuberculosis model.Additionally, we present and assess the result of the effective contact rate (β) for the susceptible, exposed, infected, treated, and recovered groups in Fig. 2.

Discussion
In this research, two experiments were conducted to investigate the dynamic behavior of the tuberculosis model.Experiment II is presented in Figs. 3 and 4, illustrating the impact of the treatment rate η on the suscepti- ble, exposed, infected, treated, and recovered population growth.The graphs revealed a decrease in the susceptible, exposed, and infected populations as the fractional order increased.Interestingly, the treated class showed an increase with higher treatment values, while the recovered class indicated a rise in population, suggesting (49 where F = (β, σ , η, µ, δ, τ ). that many individuals would recover from the disease.Notably, a higher treatment rate resulted in a decrease in the population.Furthermore, the study analyzed the effect of the contact rate, which revealed that a higher rate led to fewer susceptible individuals becoming infected.Experiment I provided sensitivity index results in Table 3, and Fig. 5 displays the graph of each parameter affecting the basic reproductive number R 0 .Through this analysis, the study identified key parameters significantly influencing the growth of the basic reproductive number.The results and graph demonstrated that increasing the treatment rate η correlated with a decline in the proportion of infected individuals seeking treatment, consequently reducing the reproductive number.Additionally, the rate at which exposed individuals become infectious τ positively impacts R 0 , highlighting the need to reduce the value of τ to control virus spread.On the other hand, σ had a nega- tive effect on the reproductive number, indicating that this rate should be decreased as well.
Stability analysis for both local and global equilibrium points was conducted.The disease-free equilibrium was found to be locally asymptotically stable when R 0 using the Routh-Hurwitz criterion for stability, implying that the tuberculosis disease would be eradicated from the population in this case.Furthermore, the endemic equilibrium was confirmed to be globally stable using the Lyapunov function, and no other equilibrium points were identified.
In conclusion, the study highlights that efficient treatment is a crucial control measure for eradicating tuberculosis.By understanding the impact of various parameters on the disease dynamics, we can develop more effective strategies to combat tuberculosis and minimize its spread in the population.

Fig. 2
Fig. 2 a-d Display the behavior of the all population at various β levels

Table 2
Parameters of the model and their respective values Estimated Olayiwola and Adedokun Bulletin of the National Research Centre (2023) 47:121