
Tuberculosis (TB) remains a major global health concern due to its infectious nature and complex treatment process. In this study, we developed a mathematical model incorporating TB progression, vaccination, latency delays, and disability outcomes. The compartmental model includes seven stages: Susceptible, vaccinated, latent, infectious, quarantined, recovered, and disabled, with time-delay terms capturing disease progression dynamics. The stability analysis of the equilibria was performed, and the sensitivity analysis was conducted using the direct differentiation method. The basic reproduction number R0 was derived to assess TB spread under different interventions. Model parameters were estimated using Ordinary Least Squares (OLS) based on Saudi Arabia's TB data (2000–2023). Numerical simulations, solved via the Adams-Bashforth-Moulton method, highlight the impact of delayed latency and quarantine on TB control, emphasizing the need for timely interventions.
Citation: Kamel Guedri, Rahat Zarin, and Mowffaq Oreijah. Evaluating the impact of vaccination and progression delays on tuberculosis dynamics with disability outcomes: A case study in Saudi Arabia[J]. AIMS Mathematics, 2025, 10(4): 7970-8001. doi: 10.3934/math.2025366
[1] | Chang Hou, Qiubao Wang . The influence of an appropriate reporting time and publicity intensity on the spread of infectious diseases. AIMS Mathematics, 2023, 8(10): 23578-23602. doi: 10.3934/math.20231199 |
[2] | Mahmoud A. Ibrahim . Threshold dynamics in a periodic epidemic model with imperfect quarantine, isolation and vaccination. AIMS Mathematics, 2024, 9(8): 21972-22001. doi: 10.3934/math.20241068 |
[3] | Moh. Mashum Mujur Ihsanjaya, Nanang Susyanto . A mathematical model for policy of vaccinating recovered people in controlling the spread of COVID-19 outbreak. AIMS Mathematics, 2023, 8(6): 14508-14521. doi: 10.3934/math.2023741 |
[4] | Ahmed M. Elaiw, Amani S. Alsulami, Aatef D. Hobiny . Global properties of delayed models for SARS-CoV-2 infection mediated by ACE2 receptor with humoral immunity. AIMS Mathematics, 2024, 9(1): 1046-1087. doi: 10.3934/math.2024052 |
[5] | Muhammad Rifqy Adha Nurdiansyah, Kasbawati, Syamsuddin Toaha . Stability analysis and numerical simulation of rabies spread model with delay effects. AIMS Mathematics, 2024, 9(2): 3399-3425. doi: 10.3934/math.2024167 |
[6] | Salma M. Al-Tuwairqi, Asma A. Badrah . Modeling the dynamics of innate and adaptive immune response to Parkinson's disease with immunotherapy. AIMS Mathematics, 2023, 8(1): 1800-1832. doi: 10.3934/math.2023093 |
[7] | Iqbal M. Batiha, Abeer A. Al-Nana, Ramzi B. Albadarneh, Adel Ouannas, Ahmad Al-Khasawneh, Shaher Momani . Fractional-order coronavirus models with vaccination strategies impacted on Saudi Arabia's infections. AIMS Mathematics, 2022, 7(7): 12842-12858. doi: 10.3934/math.2022711 |
[8] | Sayed Saber, Azza M. Alghamdi, Ghada A. Ahmed, Khulud M. Alshehri . Mathematical Modelling and optimal control of pneumonia disease in sheep and goats in Al-Baha region with cost-effective strategies. AIMS Mathematics, 2022, 7(7): 12011-12049. doi: 10.3934/math.2022669 |
[9] | Fahad Al Basir, Konstantin B. Blyuss, Ezio Venturino . Stability and bifurcation analysis of a multi-delay model for mosaic disease transmission. AIMS Mathematics, 2023, 8(10): 24545-24567. doi: 10.3934/math.20231252 |
[10] | Feliz Minhós, Ali Raza, Umar Shafique . An efficient computational analysis for stochastic fractional heroin model with artificial decay term. AIMS Mathematics, 2025, 10(3): 6102-6127. doi: 10.3934/math.2025278 |
Tuberculosis (TB) remains a major global health concern due to its infectious nature and complex treatment process. In this study, we developed a mathematical model incorporating TB progression, vaccination, latency delays, and disability outcomes. The compartmental model includes seven stages: Susceptible, vaccinated, latent, infectious, quarantined, recovered, and disabled, with time-delay terms capturing disease progression dynamics. The stability analysis of the equilibria was performed, and the sensitivity analysis was conducted using the direct differentiation method. The basic reproduction number R0 was derived to assess TB spread under different interventions. Model parameters were estimated using Ordinary Least Squares (OLS) based on Saudi Arabia's TB data (2000–2023). Numerical simulations, solved via the Adams-Bashforth-Moulton method, highlight the impact of delayed latency and quarantine on TB control, emphasizing the need for timely interventions.
Tuberculosis is a chronic infectious disease that is a major cause of ill health and one of the leading causes of mortality worldwide. The causative agent of TB is the bacterium Mycobacterium tuberculosis (M.tb), which typically attacks the lungs but can also affect other parts of the body, such as the kidneys, spine, and brain. TB is spread through the air when actively infected people cough, sneeze, or spit. According to the Global Tuberculosis Report 2022 [1], about a quarter of the global population is estimated to have been infected with TB as latent tuberculosis infection (LTBI) asymptomatic and non-infectious but only about 5–10% of infected individuals eventually develop active TB disease, which is symptomatic and infectious. Some individuals may also naturally clear the infection. Without treatment, the death rate from active TB disease is high (about 50%), although 59–95% of patients can be cured with recommended treatments, which usually involve a standard 6–24 month course of four antibiotics, including rifampicin and isoniazid [2]. The only licensed vaccine for TB prevention is the Bacille Calmette-Guérin (BCG) vaccine. However, BCG's effectiveness is not lifelong, as immunity weakens over time. Moreover, not all vaccinated individuals are fully protected, and the probability of developing active TB is higher among people with weakened immune systems, particularly those living with HIV [3].
Mathematical modeling of infectious disease transmission is widely used to provide valuable insights for public health policies aimed at preventing or reducing disease spread [4,5,6,7,8,9]. Researchers have further advanced epidemic modeling by incorporating artificial neural networks and nonlinear incidence structures to simulate spatial diffusion and real-data-driven disease dynamics [10,11]. Alqahtani et al. [12] investigated the impact of mobile teams on TB treatment outcomes in the Riyadh Region of Saudi Arabia from 2013 to 2015. Their study demonstrated that mobile teams improved treatment success rates by ensuring better patient adherence and access to healthcare services. Nave et al. [13] explored mathematical modeling in the context of cancer treatment, specifically evaluating the effectiveness of a combination therapy for breast cancer. While they focus on a different disease, the modeling techniques employed could provide useful insights for optimizing TB treatment strategies. Zhang et al. [14] examined stability in complex networks using stochastic delayed models, contributing to the broader mathematical modeling field. Similarly, Zhang et al. [15] explored synchronization in stochastic systems with time-delay effects, while Zhang et al. [16] investigated control mechanisms in semi-Markov jump systems with distributed delay. However, TB's transmission is complex, and many aspects of its natural history and dynamics remain unclear [17]. Extensive research has been conducted on mathematical modeling and analysis of TB. Waaler et al. [18] introduced the first mathematical model to investigate TB epidemiological trends. Li et al. [19] later developed a TB dynamical model, fitted it to data to estimate parameters, and discussed optimal control strategies. Das et al. [20] constructed a model incorporating media impact on TB transmission rates. Numerous factors have since been incorporated into TB models, including vaccination [21,22], treatment and incomplete treatment [23,24], fast and slow progression [24,25], relapse [25], reinfection [26], co-infection with HIV [21], and drug-resistant strains [24,26].
According to the WHO, preventive treatment can be provided for people with LTBI to stop the progression to active infection [1]. At the UN high-level meeting on TB, a goal was set to provide preventive treatment to 30 million people worldwide between 2018 and 2022 [1]. This preventive treatment uses the same drugs as active TB treatment but for a shorter duration, with options now lasting only 3 to 6 months instead of the 6-month minimum required for active TB. Compared to treating active infections, preventive treatment for LTBI is more economical and poses less risk of further development into active disease. The Centers for Disease Control and Prevention (CDC) in the United States also supports preventive treatment for LTBI, highlighting its importance in controlling TB by substantially reducing the risk of progression to active disease [3]. This growing emphasis on preventive treatment has sparked interest in further modeling its impact on TB dynamics. Few models consider preventive treatment for people with LTBI. Castillo-Chavez and Feng [23], Bhunu et al. [21], and Zhang et al. have incorporated a treatment term rE to represent the movement from exposed (E) to recovered (R) individuals. However, neither preventive treatment for LTBI patients nor routine treatment for active TB is a short-term process, as both require several months. For more realistic modeling, it is therefore appropriate to place patients undergoing treatment in a separate compartment. Additionally, some individuals receiving preventive treatment may progress to active TB due to factors like reduced immunity [1]. Developing mathematical models that account for preventive treatment is thus essential for understanding TB control.
We develop a novel TB transmission model incorporating both TB's progression and key potential intervention strategies given the TB's complex dynamics especially because of its chronic nature and potential for grave outcome. This model describes important stages of TB infection (the latent period and risk of disabilities through severe or untreated TB). However, by enabling time delays and a disability compartment, the resulting framework encompasses a multi-dimensional set of TB prevention and control effort effects including vaccination and on time treat, traversing the range of impacts of TB and how interventions reduce infection rates and long term complications. Our model divides the population into seven compartments: The susceptible (S), vaccinated (V), latent (L), infectious (I), quarantined or under treatment (Q), recovered (R), and disabled (D), which are several stages of progressive TB or intervention. When exposed to infectious persons, susceptible persons become infected, enter a latent period, and then potentially become infectious. People in the vaccinated compartment have partial immunity, especially against severe forms of TB, but may gradually lose this immunity and relocate to the susceptible compartment. Our model includes time delay terms to represent the time the latent compartment takes to incubate, and the time between symptom onset and severe complications in the disabled compartment.
The primary focus of this model is to explore how preventive measures, such as vaccination and timely quarantine or treatment, impact TB transmission dynamics, while accounting for severe cases that may lead to long-term disabilities. We analyze the flow of individuals between compartments and use sensitivity analysis to assess how intervention rates, immunity waning, and time delays influence disease outcomes. Our findings provide new insights into effective TB control strategies, underscoring the importance of vaccination, rapid treatment initiation, and disability prevention as components of comprehensive TB management. This model offers a realistic framework for evaluating TB dynamics and intervention impacts, supporting policymakers and health officials in designing targeted TB control policies.
Tuberculosis can lead to long-term complications and disabilities, particularly in cases where the disease is not diagnosed or treated in a timely manner. Severe TB infections, especially drug-resistant or extrapulmonary TB, can result in permanent lung damage, musculoskeletal deformities (e.g., spinal TB leading to kyphosis), and neurological impairments (e.g., tuberculous meningitis causing cognitive and motor function loss). These long-term consequences can significantly impact the quality of life of affected individuals and impose a substantial burden on healthcare systems. In Saudi Arabia, TB remains a public health concern, and although incidence rates have declined, certain populations such as migrant workers and individuals with comorbidities like diabetes are at higher risk of developing severe TB-related complications. Studies have highlighted that delayed treatment initiation and inadequate adherence to TB therapy increase the likelihood of disease progression, leading to disability outcomes. By incorporating a disability compartment in our model, we aim to capture these long-term health effects and provide a more realistic framework for evaluating the burden of TB, beyond just infection and recovery dynamics.
In this study, we present a compartmental model to simulate TB transmission dynamics, incorporating key stages of the disease and the effects of interventions, such as vaccination and the potential for disability from severe TB cases. Given the chronic nature of TB, which includes both a latent period and potential delays in treatment initiation, the model integrates time delay terms to capture these critical aspects. The population is divided into seven compartments: Susceptible, vaccinated, latent, infectious, quarantined, recovered, and disabled individuals, as shown in Figure 1.
Susceptible individuals S(t) are those at risk of infection, and new susceptible individuals enter the population at a recruitment rate Λ. Susceptible individuals may become infected through contact with infectious individuals at rate βS(t)I(t), where β represents the transmission rate. Additionally, susceptible individuals may be vaccinated at rate ρ, which moves them to the vaccinated compartment (V(t)), and they may lose immunity over time and return to the susceptible pool at rate ε. Vaccinated individuals experience partial immunity, especially against severe TB forms, but this immunity wanes over time. They may re-enter the susceptible class or die naturally at rate μ. Susceptible individuals also have the possibility of reinfection from recovered individuals as immunity wanes, which is represented by the rate φ. When a susceptible individual is infected, they enter the latent compartment L(t), where they carry the TB bacterium but are not yet infectious. This latent stage is crucial for TB dynamics, as individuals can remain in this state for months or years before progressing to active disease. A time delay term τ is applied to the progression rate θ, reflecting this incubation period where the individual is non-infectious but is a carrier of the bacterium. Some latent individuals may naturally clear the infection at rate δ, while others progress to the infectious compartment I(t) after the delay τ. Infectious individuals represent those with active TB who can transmit the disease to others. Infectious individuals may enter quarantine or treatment at rate α, recover at rate η1, or experience progression to the disabled compartment D(t) due to complications from untreated or prolonged TB, at rate ω. Disease-induced deaths occur among the infectious at rate d2, along with the natural death rate μ.
The quarantine or treatment compartment Q(t) represents individuals who are diagnosed and isolated for treatment. These individuals have reduced infectiousness due to their isolation and may recover with treatment, entering the recovered compartment R(t) at rate η2. However, some quarantined individuals may still experience severe complications that result in disability, represented by rate κ, with a delay term σ capturing the time from symptom onset to the development of these complications. Both infectious and quarantined compartments experience disease-induced deaths at rate d2 and natural deaths at rate μ. Recovered individuals R(t) represent those who have cleared the infection and developed partial immunity to TB. However, immunity from TB is not always lifelong, so a fraction of recovered individuals may lose immunity over time and re-enter the susceptible compartment at rate φ. The recovered compartment also experiences natural deaths at rate μ. Finally, the disabled compartment (D(t)) captures individuals who have developed long-term disabilities due to TB complications. These disabilities can result from either direct progression from the infectious class or from complications arising during quarantine. Disabled individuals may recover over time at rate γ, returning to the susceptible class, or face natural mortality at rate μ. The equations governing this model are as follows:
{dS(t)dt=Λ−βS(t)I(t)−(μ+ρ)S(t)+φR(t)+εV(t)dV(t)dt=ρS(t)−(μ+ε)V(t)dL(t)dt=βS(t)I(t)−(μ+δ)L(t)−θL(t−τ)dI(t)dt=θL(t−τ)−(μ+d2+η1+α)I(t)−ωI(t)dQ(t)dt=αI(t)−(μ+d2+η2)Q(t)−κQ(t−σ)dR(t)dt=η1I(t)+η2Q(t)−(μ+φ)R(t)dD(t)dt=ωI(t)+κQ(t−σ)−(μ+γ)D(t)S(0),V(0),L(0),I(0),Q(0),R(0),D(0)≥0. | (2.1) |
Here, the time delay τ in the latent-to-infectious transition and σ in the quarantine-to-disability transition are key features of this model, capturing the natural course of TB and its progression. The model parameters are defined as follows: The recruitment rate (Λ), the death rate (μ), and the transmission rate (β) of TB from infectious individuals to susceptibles. The vaccination rate is denoted by (ρ), while the recovery rate for vaccinated individuals is given by (φ). The waning immunity rate for both recovered and vaccinated individuals is described by (ε). The rate at which latent infections clear is indicated by (δ), and the delay between the latent and infectious stages of TB is represented by (τ), while the rate of progression from latency to infectiousness after the delay is given by (θ). The natural clearance of latent infection is represented by (δ). Infectious individuals may move to quarantine at a rate of (α), recover at a rate of (η1), or become disabled due to severe complications at a rate of (ω). Upon leaving the quarantined class, individuals may recover at a rate of (η2) but could also migrate to the disabled class at a rate of (κ), following a delay (σ) that reflects the development of delayed complications. The mortality rate for individuals in the infected or quarantined classes is denoted by (d2), and the recovery rate from disability is given by (γ). This model provides a structured approach to analyze TB transmission dynamics, the effects of vaccination, and time delays in both latency and treatment, as well as the potential for severe complications leading to disability. The model offers a realistic framework for understanding the progression of TB and the impacts of interventions within a population.
System (2.1) models the dynamics of a human population across different health states. To guarantee that the state variables S(t),V(t),L(t),I(t),Q(t),R(t),D(t) remain non-negative for all t≥0, it is essential to establish their positivity throughout the time domain.
Lemma 3.1. Let U⊆Rn be an open set, and consider a solution x=(x1,…,xn)∈C1((0,T);U)∩C([0,T);U) to the ordinary differential equation:
{˙x(t)=f(t,x(t)),x(0)=x0, |
where the initial condition x0=(x1,0,…,xn,0)∈U satisfies xi,0>0 for all i∈{1,…,n}, and f=(f1,…,fn):(0,T)×U→Rn. Suppose that whenever yi=0 and all other yk>0 for k≠i, it holds that fi(t,y)>0 for t∈(0,T). Then, the solution x remains positive for all t∈[0,T), meaning xi(t)>0 for each i∈{1,…,n} and t∈[0,T) [27].
Proof. Assume that x is not strictly positive at all times. By continuity, there exists a time s∈(0,T) and an index j∈{1,…,n} such that xj(s)=0. Thus, the set of times at which any xi reaches zero, denoted as ∪ni=1x−1i({0}), forms a closed, non-empty subset that has a lower bound, indicating a minimum zero at ˉτ∈(0,T). Consequently, xi(ˉτ)≥0, while xi(t)>0 for t∈(0,ˉτ) for all i∈{1,…,n}.
Now, evaluating x at t=ˉτ:
● If x1(ˉτ)=0, then by assumption, f1(t,x(t))>0 when evaluated at x1(ˉτ)=0 and (x2(ˉτ),…,xn(ˉτ))∈U∩Rn−1≥0.
● Consequently, ˙x1(ˉτ)=f1(t,x(t))>0.
Since ˙x1 is continuous, there exists a small interval (ˉτ−ˉε,ˉτ+ˉε) where ˙x1(t)>0, leading to:
x1(ˉτ−ˉε)=x1(ˉτ)−∫ˉτˉτ−ˉε˙x1(t)dt<0, |
which contradicts the assumption that x1(t)≥0 for t∈(0,ˉτ). Thus, x1(ˉτ)>0.
For indices i=2,…,n, assume (x1(ˉτ),…,xi−1(ˉτ))∈U∩Ri−1>0 and (xi(ˉτ),…,xn(ˉτ))∈U∩Rn−i+1≥0. If xi(ˉτ)=0, then by a similar argument, fi(t,x(t))>0, which results in a contradiction. By induction, it follows that xi(ˉτ)>0 for all i, proving that the solution x remains positive for all t.
Theorem 3.1. The solutions S(t),V(t),L(t),I(t),Q(t),R(t),D(t) of system (2.1) exist for all time and are always positive, smooth, and unique for all t>0, subject to the initial conditions S0>0,V0>0,L0>0,I0>0,Q0>0,R0>0,D0>0.
Proof. Consider the ODE system with delay terms.
x=(x1x2x3x4x5x6x7)=(SVLIQRD),x0=(S0V0L0I0Q0R0D0) |
and
f=(f1f2f3f4f5f6f7)=(Λ−βSI−(μ+ρ)S+φR+εVρS−(μ+ε)VβSI−(μ+δ)L−θL(t−τ)θL(t−τ)−(μ+d2+η1+α)I−ωIαI−(μ+d2+η2+κ)Q(t−σ)η1I+η2Q−(μ+φ)RωI+κQ(t−σ)−(μ+γ)D). |
This defines the state vector x with initial values x0 and the functional vector f, where each component fi characterizes the rate of change of the corresponding state variable.
Consider the open set U=BN+1, where BN+1 is a ball of radius N+1, with N denoting the total human population. The system remains confined within this set due to population constraints. Since the initial state satisfies x0∈U, and the function f is bounded, Lipschitz continuous, and smooth in U, the "Picard-Lindelöf theorem" [28] guarantees the existence of a unique solution x(t)=(x1(t),…,x7(t))∈C1((0,T);U)∩C([0,T);U) over a short time interval T>0, given the specified initial conditions. Since the system adheres to the conditions outlined in Lemma 3.1, the state variables S(t),V(t),L(t),I(t),Q(t),R(t),D(t) are ensured to remain positive. Additionally, summing the compartmental variables, we obtain:
N=S(t)+V(t)+L(t)+I(t)+Q(t)+R(t)+D(t), |
which implies that x(t)∈¯BN⋐U. This guarantees that the solution remains positive, uniquely determined, and globally well-defined for all t≥0. By applying "Corollary 17.4" from [29], we establish the existence of a unique global positive solution.
Furthermore, utilizing a bootstrap argument, we demonstrate the smoothness of the solution over time. Since f is smooth, the solution satisfies the integral equation:
x(t)=x0+∫t0f(x(s))ds. |
If x(t) belongs to Cn, then differentiating under the integral sign implies that x(t) belongs to Cn+1. By induction, it follows that x(t) is infinitely differentiable, i.e., C∞, confirming the smoothness of the solution over time.
We examine system (2.1) within the biologically relevant feasible region DH.
Theorem 3.2. The feasible region
DH={(S(t),V(t),L(t),I(t),Q(t),R(t),D(t))∈R7+:S+V+L+I+Q+R+D≤Λμ} |
is positively invariant and attracting for system (2.1) for all t≥0.
Proof. Define N=S+V+L+I+Q+R+D. Taking the time derivative, we get:
dNdt=dSdt+dVdt+dLdt+dIdt+dQdt+dRdt+dDdt. |
By substituting the respective expressions from system (2.1), we obtain:
dNdt=Λ−μN−δL−d2I−d2Q−γD, |
where all terms associated with mortality and progression between compartments remain bounded. This simplifies to:
dNdt+μN≤Λ. |
Using the integrating factor If=eμt, we reformulate the equation as:
ddt[N(t)eμt]≤Λeμt. |
Integrating both sides from 0 to t, we obtain:
N(t)eμt−N(0)≤Λ(eμt−1). |
Rearranging, we derive:
N(t)≤Λμ+[N(0)−Λμ]e−μt. |
For all t≥0, it follows that N(t)≤Λμ whenever N(0)≤Λμ. If N(0)>Λμ, then N(t) asymptotically approaches DH, confirming that DH is both positively invariant and attracts all trajectories in R7+. Thus, system (2.1) remains well-posed, both mathematically and biologically, within the feasible human population region DH.
With the disease-free equilibrium (DFE) now defined, we can compute the basic reproduction number R0 using the next-generation matrix approach [30]. The DFE is given by:
(S,V,R,L,I,Q,D)=(S0,V0,0,0,0,0,0), |
where
S0=Λμ+ρ−ερμ+ε,V0=ρS0μ+ε. |
Several studies have estimated the basic reproduction number R0 for tuberculosis in different regions. For example, Li et al. [31] and Wu et al. [32] reported values ranging from 1.6 to 2.3 for TB in the United States, while Zhang et al. [33] found R0≈1.68 for TB transmission in China. Additionally, a systematic review by Ma et al. [34] concluded that R0 for TB generally lies between 1.0 and 4.3, depending on model structure and public health factors. The R0 obtained from our model lies within this expected range, supporting the validity of our parameter choices and the realism of our transmission assumptions. Using the next-generation matrix technique, we define F and V to represent the new infection and transition matrices, respectively. The infection matrix F evaluated at the DFE becomes:
F=[0βΛμ+ρ−ερμ+ε00000000000000]. |
The transition matrix V, which describes the rates at which individuals move out of the infected compartments, is:
V=[μ+δ+θ000−θμ+d2+η1+α+ω000−αμ+d2+η2+κ00−ω−κμ+γ]. |
The inverse of V is given by:
V−1=[1μ+δ+θ000−θ(μ+δ+θ)(μ+d2+η1+α+ω)1μ+d2+η1+α+ω000−α(μ+d2+η1+α+ω)(μ+d2+η2+κ)1μ+d2+η2+κ00−ω(μ+d2+η1+α+ω)(μ+γ)−κ(μ+d2+η2+κ)(μ+γ)1μ+γ]. |
Multiplying F by V−1, we obtain:
FV−1=[Λβθ(μ+δ+θ)(μ+ρ−ερμ+ε)(α+μ+d2+η1+ω)−Λβ(μ+ρ−ερμ+ε)(α+μ+d2+η1+ω)00000000000000]. |
The dominant eigenvalue (or spectral radius) of FV−1 gives the basic reproduction number R0, which is:
R0=Λβθ(μ+δ+θ)(μ+ρ−ερμ+ε)(α+μ+d2+η1+ω). |
This represents the basic reproduction number R0 for the system (2.1).
At equilibrium, we set all derivatives to zero:
dS(t)dt=dV(t)dt=dL(t)dt=dI(t)dt=dQ(t)dt=dR(t)dt=dD(t)dt=0. |
The equilibrium equations are:
0=Λ−βSI−(μ+ρ)S+φR+εV,0=ρS−(μ+ε)V,0=βSI−(μ+δ+θ)L,0=θL−(μ+d2+η1+α+ω)I,0=αI−(μ+d2+η2+κ)Q,0=η1I+η2Q−(μ+φ)R,0=ωI+κQ−(μ+γ)D. | (4.1) |
From the second equation of system (4.1):
V=ρSμ+ε. | (4.2) |
From the third equation of system (4.1):
L=βSIμ+δ+θ. | (4.3) |
From the fourth equation of system (4.1):
I=θLμ+d2+η1+α+ω. |
Substitute L from Eq (4.3):
I=θμ+d2+η1+α+ω⋅βSIμ+δ+θ. |
Rearrange to isolate I:
I[1−θβS(μ+d2+η1+α+ω)(μ+δ+θ)]=0. |
If I≠0:
S=(μ+d2+η1+α+ω)(μ+δ+θ)θβ. | (4.4) |
From the fifth equation of system (4.1):
Q=αIμ+d2+η2+κ. | (4.5) |
From the 6th equation of system (4.1):
R=η1I+η2Qμ+φ. |
Substitute Q from Eq (4.5):
R=η1I+η2⋅αIμ+d2+η2+κμ+φ. |
Simplify:
R=η1(μ+d2+η2+κ)I+η2αI(μ+φ)(μ+d2+η2+κ). |
Factor out I:
R=Iμ+φ⋅η1(μ+d2+η2+κ)+η2αμ+d2+η2+κ. | (4.6) |
From the last equation of system (4.1):
D=ωI+κQμ+γ. |
Substitute Q from Eq (4.5):
D=ωI+κ⋅αIμ+d2+η2+κμ+γ. |
Simplify:
D=Iμ+γ⋅[ω+καμ+d2+η2+κ]. | (4.7) |
Thus, the equilibrium is given by:
S∗=(μ+d2+η1+α+ω)(μ+δ+θ)θβ,V∗=ρS∗μ+ε,L∗=βS∗I∗μ+δ+θ,I∗=θL∗μ+d2+η1+α+ω,Q∗=αI∗μ+d2+η2+κ,R∗=I∗μ+φ⋅η1(μ+d2+η2+κ)+η2αμ+d2+η2+κ,D∗=I∗μ+γ⋅[ω+καμ+d2+η2+κ]. |
Notes:
(1) S∗ is determined independently from the system parameters and does not depend on I∗ or L∗.
(2) L∗ and I∗ are interdependent but consistent through substitution.
(3) All other variables (V∗,Q∗,R∗,D∗) depend on I∗ and system parameters.
Theorem 5.1. The disease-free equilibrium of the proposed model (2.1) is locally asymptotically stable (LAS) if R0<1.
Proof. The Jacobian matrix of the model at the disease-free equilibrium is given by:
J=[−(μ+ρ)ε0−βS00φ0ρ−(μ+ε)0000000−(μ+δ)−θe−ζτβS000000θe−ζτ−(μ+d2+η1+α+ω)000000α−(μ+d2+η2+κe−ζσ)00000η1η2−(μ+φ)0000ωκe−ζσ0−(μ+γ)]. |
The last eigenvalue of the above matrix is clearly negative, i.e., λ=−(μ+γ), and for the rest of the eigenvalue, we use the following reduced matrix:
J=[−(μ+ρ)ε0−βS00φρ−(μ+ε)000000−(μ+δ)−θe−ζτβS00000θe−ζτ−(μ+d2+η1+α+ω)00000α−(μ+d2+η2+κe−ζσ)0000η1η2−(μ+φ)]. |
We now systematically reduce the given 6×6 matrix to an upper triangular form using Gaussian elimination (row operations).
Step 1: Eliminate ρ in the second row, first column. The first pivot is at a11=−(μ+ρ). To eliminate ρ in the second row, first column, we perform the row operation:
R2←R2−ρμ+ρR1. |
Step 2: Eliminate θe−ζτ in the fourth row, third column. The third pivot is at a33=−(μ+δ)−θe−ζτ. To eliminate the non-zero term θe−ζτ in the fourth row, third column, we perform the row operation:
R4←R4−θe−ζτμ+δ+θe−ζτR3. |
Step 3: Eliminate α in the fifth row, fourth column. The fourth pivot is now at a44. To eliminate α in the fifth row, fourth column, we perform the row operation:
R5←R5−αa44R4. |
Step 4: Eliminate η1 in the sixth row, fourth column. To eliminate η1 in the sixth row, fourth column, we perform the row operation:
R6←R6−η1a44R4. |
Step 5: Eliminate η2 in the sixth row, fifth column. The fifth pivot is now at a55=−(μ+d2+η2+κe−ζσ). To eliminate η2 in the sixth row, fifth column, we perform the row operation:
R6←R6−η2a55R5. |
Thus, after performing all the row operations, the final upper triangular matrix is:
Jupper=[−(μ+ρ)ε0−βS00φ0−(μ+ε)(μ+ρ)+ρεμ+ρ0−ρβS0μ+ρ0ρφμ+ρ00−(μ+δ)−θe−ζτβS000000D44000000−(μ+d2+η2+κe−ζσ)000000−(μ+φ)] |
where
D44=−(μ+d2+η1+α+ω)−θe−ζτ⋅βS0−(μ+δ)−θe−ζτ. |
Since the determinant of an upper triangular matrix is simply the product of its diagonal elements, the eigenvalues of the system are:
λ1=−(μ+ρ),λ2=−(μ+ε)(μ+ρ)+ρεμ+ρ,λ3=−(μ+δ)−θe−ζτ,λ4=−(μ+d2+η1+α+ω)−θe−ζτ⋅βS0−(μ+δ)−θe−ζτ,λ5=−(μ+d2+η2+κe−ζσ),λ6=−(μ+ϕ),λ7=−(μ+γ). |
The following theorem establishes the stability of the endemic equilibrium.
Theorem 5.2. The endemic equilibrium of the proposed model (2.1) is locally asymptotically stable (LAS) if R0>1.
Proof. The Jacobian matrix of the system (2.1) evaluated at the endemic equilibrium is given by:
J=[−βI∗−(μ+ρ)ε0−βS∗0φ0ρ−(μ+ε)00000βI∗0−(μ+δ)−θe−ζτβS∗00000θe−ζτ−(μ+d2+η1+α+ω)000000α−(μ+d2+η2+κe−ζσ)00000η1η2−(μ+φ)0000ωκe−ζσ0−(μ+γ)]. |
The last eigenvalue of the above matrix is clearly negative, i.e., λ=−(μ+γ), and for the rest of the eigenvalue, we use the following reduced matrix:
J=[−βI∗−(μ+ρ)ε0−βS∗0φρ−(μ+ε)0000βI∗0−(μ+δ)−θe−ζτβS∗0000θe−ζτ−(μ+d2+η1+α+ω)00000α−(μ+d2+η2+κe−ζσ)0000η1η2−(μ+φ)]. |
Using elementary row operations,
Step 1: Eliminate ρ in row 2, column 1: R2←R2−ρβI∗+μ+ρR1.
Step 2: Eliminate βI∗ in row 3, column 1: R3←R3−βI∗βI∗+μ+ρR1.
Step 3: Eliminate θe−ζτ in row 4, column 3: R4←R4−θe−ζτμ+δ+θe−ζτR3.
Step 4: Eliminate α in row 5, column 4: R5←R5−αa44R4.
Step 5: Eliminate η1 in row 6, column 4: R6←R6−η1a44R4.
Step 6: Eliminate η2 in row 6, column 5: R6←R6−η2a55R5.
Step 7: Eliminate the remaining non-zero entry in row 3, column 2: R3←R3−R3[2]R2[2]R2.
We obtained the following upper triangular matrix:
Jupper=[−(βI∗+μ+ρ)ε0−βS∗0φ0a220R2[4]0R2[6]00a33R3[4]0R3[6]000a440R4[6]0000a55R5[6]00000R6[6]] |
where
a22=−(μ+ε+ρεβI∗+μ+ρ),a33=−(μ+δ)−θe−ζτR2[4]=ρβS∗βI∗+μ+ρ,R2[6]=−ρφβI∗+μ+ρR3[4]=βS∗+β2I∗S∗βI∗+μ+ρ+βI∗εa22(βI∗+μ+ρ)⋅R2[4]R3[6]=−βI∗φβI∗+μ+ρ+βI∗εa22(βI∗+μ+ρ)⋅R2[6]a44=−(μ+d2+η1+α+ω)+θe−ζτ⋅R3[4]μ+δ+θe−ζτR4[6]=θe−ζτ⋅R3[6]μ+δ+θe−ζτ,a55=−(μ+d2+η2+κe−ζσ)R5[6]=−α⋅R4[6]a44,R6[6]=−(μ+φ)−η1⋅R4[6]a44+η2⋅R5[6]a55. |
The diagonal entries of the matrix Jupper are the corresponding eigenvalues.
Precise parameter estimation is fundamental in epidemiological modeling, as it underpins the model's ability to replicate real-world data with high accuracy. Among various methods, the Least Squares Curve Fitting technique is widely utilized to refine model parameters, ensuring the alignment of model outputs with empirical data. The text is clear and well-structured but could be refined slightly for clarity and conciseness. The computational efficiency and reliability of this technique rely on the assumptions that residuals are normally distributed with constant variance, making it well-suited for obtaining accurate parameter estimates under these conditions. We analyze data on confirmed cases in Saudi Arabia from 2000 to 2023 to study tuberculosis using the Least Squares Curve Fitting approach. Specifically, the Ordinary Least Squares (OLS) method is applied to minimize discrepancies in the daily reported cases. The model's accuracy is assessed using a relative error metric, as shown below:
min(∑nι=1(Iι−ˆIι)2∑nι=1I2ι). |
In this expression, Iι denotes the model-derived cumulative count of infected individuals, based on daily transitions from infection to recovery. Moreover, ˆIι represents the observed cumulative infections reported. This distinction between model-simulated data and real-world data enables a more nuanced analysis of disease dynamics, enhancing predictions of disease spread and informing more effective disease control strategies.
In the following section, we describe the application of the least squares curve-fitting method in our epidemiological model to ensure its compatibility with empirical data. Population normalization across compartments relative to the total population, N(t), is applied for this analysis. By employing the least squares method, we accurately estimate biological parameters, as shown in Figure 2. This process minimizes the mean absolute relative error, aligning the model's forecasts with the actual incidence rates of tuberculosis. Real case data are represented by red solid circles, while the optimally fitted model curve is depicted in blue. Table 1 provides a detailed account of these optimized biological parameter estimates obtained through the least squares approach. This methodological rigor enhances the predictive precision of the model, offering a robust basis for understanding tuberculosis epidemic dynamics.
Parameters | Description | Values | Source |
Λ | Recruitment rate of susceptible individuals | 1000/year | Assumed |
μ | Natural death rate | 0.05/year | Literature |
β | Transmission rate from infectious individuals to susceptibles | 1.4137×10−04year | Estimated |
ρ | Vaccination rate of susceptible individuals | 0.65/year | Fitted |
ε | Rate of waning immunity for vaccinated individuals | 0.55/year | Estimated |
φ | Rate of waning immunity for recovered individuals | 0.05/year | Assumed |
δ | Rate of latent infection clearance | 0.001/year | Estimated |
θ | Progression rate from latent to infectious TB (with delay τ) | 0.45/year | Estimated |
α | Rate at which infectious individuals enter quarantine or treatment | 0.1/year | Estimated |
η1 | Recovery rate for infectious individuals | 0.35/year | Literature |
ω | Rate of progression from infectious to disabled due to complications | 0.1/year | Estimated |
d2 | Disease-induced death rate for infectious and quarantined individuals | 0.035/year | Estimated |
η2 | Recovery rate for quarantined individuals | 0.07/year | Literature |
κ | Rate of progression from quarantine to disability (with delay σ) | 0.01/year | Estimated |
γ | Recovery rate for disabled individuals | 0.1/year | Estimated |
τ | Delay in progression from latent to infectious TB | [0−15]/day | Fitted |
σ | Delay in progression from quarantine to disability | [0−12]/day | Fitted |
Identifying key parameters that influence the transmission dynamics of infectious diseases is crucial for devising effective control strategies. This is achieved through sensitivity analysis, which quantifies how variations in model parameters impact disease spread. Among various approaches, forward sensitivity analysis plays a pivotal role in epidemiological modeling. However, its computation can become increasingly complex for intricate biological systems. The sensitivity analysis of the basic reproduction number R0 has garnered significant attention from both ecologists and epidemiologists. To assess the impact of various parameters on R0, we conducted a sensitivity analysis using the direct differentiation method. The sensitivity index ΥR0Ω of a parameter Ω is given by:
ΥR0Ω=∂R0∂Ω⋅ΩR0. |
The influence of these parameters on the basic reproduction number R0 is illustrated graphically in Figure 3.
In this section, we apply the Adams-Bashforth-Moulton (ABM) predictor-corrector method to solve the system of delay differential equations (DDEs) [36,37,38,39,40]. The process consists of two steps: First, predicting the values at the next time step using an explicit Adams-Bashforth (AB) predictor, and second, correcting these predictions with the implicit Adams-Moulton (AM) corrector. While we employ the ABM predictor-corrector method for solving the system of DDEs, alternative numerical methods are also available. These include the Runge-Kutta methods for DDEs, Backward Differentiation Formulas (BDF), and spectral collocation methods, which are commonly used for solving both stiff and non-stiff DDEs. The choice of the ABM method in this study is based on its balance between computational efficiency and accuracy. It provides an explicit predictor step followed by an implicit corrector step, which enhances stability while maintaining a relatively low computational cost. This makes it particularly suitable for our TB model, where delay terms play a crucial role in capturing disease progression dynamics. Future work may explore more advanced adaptive-step solvers for improved precision.
Given the model equations:
dX(t)dt=f(X,t), | (8.1) |
where
X(t)=[S(t),V(t),L(t),I(t),Q(t),R(t),D(t)]T, |
and each derivative is represented by the functions on the right-hand side of each corresponding equation. Here is the step-by-step ABM scheme in mathematical form for each equation.
Step 1: Adams-Bashforth Predictor (Explicit)
(1) Predict: The values Xp(tn+1) at the next time step tn+1=tn+h using the Adams-Bashforth formula. For simplicity, we'll use a two-step AB formula:
Xp(tn+1)=X(tn)+h2(3f(X(tn),tn)−f(X(tn−1),tn−1)). | (8.2) |
This gives a predicted value Xp(tn+1) for each component in X(t). We substitute each function fi corresponding to the derivatives given in the model equations.
Step 2: Adams-Moulton Corrector (Implicit)
(2) Correct: The predicted values Xp(tn+1) using the Adams-Moulton corrector formula:
X(tn+1)=X(tn)+h12(5f(Xp(tn+1),tn+1)+8f(X(tn),tn)−f(X(tn−1),tn−1)). | (8.3) |
This corrects each component of X(tn+1) using the predicted value Xp(tn+1) and the previous function evaluations at tn and tn−1.
Let's apply this ABM scheme to each equation. For simplicity, denote:
● Sn≈S(tn), Vn≈V(tn), and so on.
● fS(tn)=Λ−βSnIn−(μ+ρ)Sn+φRn+εVn (for dS(t)dt), etc.
Thus, we have:
fS(tn)=Λ−βSnIn−(μ+ρ)Sn+φRn+εVn,fV(tn)=ρSn−(μ+ε)Vn,fL(tn)=βSnIn−(μ+δ)Ln−θL(tn−τ),fI(tn)=θL(tn−τ)−(μ+d2+η1+α)In−ωIn,fQ(tn)=αIn−(μ+d2+η2+κ)Q(tn−σ),fR(tn)=η1In+η2Qn−(μ+φ)Rn,fD(tn)=ωIn+κQ(tn−σ)−(μ+γ)Dn. |
Predictor Step
Sp,n+1=Sn+h2(3fS(tn)−fS(tn−1)), |
Vp,n+1=Vn+h2(3fV(tn)−fV(tn−1)), |
Lp,n+1=Ln+h2(3fL(tn)−fL(tn−1)), |
Ip,n+1=In+h2(3fI(tn)−fI(tn−1)), |
Qp,n+1=Qn+h2(3fQ(tn)−fQ(tn−1)), |
Rp,n+1=Rn+h2(3fR(tn)−fR(tn−1)), |
Dp,n+1=Dn+h2(3fD(tn)−fD(tn−1)). |
Corrector Step
Sn+1=Sn+h12(5fS(tn+1)+8fS(tn)−fS(tn−1)), |
Vn+1=Vn+h12(5fV(tn+1)+8fV(tn)−fV(tn−1)), |
Ln+1=Ln+h12(5fL(tn+1)+8fL(tn)−fL(tn−1)), |
In+1=In+h12(5fI(tn+1)+8fI(tn)−fI(tn−1)), |
Qn+1=Qn+h12(5fQ(tn+1)+8fQ(tn)−fQ(tn−1)), |
Rn+1=Rn+h12(5fR(tn+1)+8fR(tn)−fR(tn−1)), |
Dn+1=Dn+h12(5fD(tn+1)+8fD(tn)−fD(tn−1)). |
For delayed terms like L(t−τ) and Q(t−σ):
● Interpolation: If t−τ or t−σ lies between two time points, use interpolation (e.g., linear interpolation) to approximate the delayed values.
● Direct Use: If t−τ or t−σ coincides with a previous time step, use the stored values directly.
Figure 4 presents simulation results of the model, illustrating the time evolution of each state variable over a span of 1200 days. The dynamics exhibit oscillatory behavior, with each state variable converging to equilibrium levels. Specifically, the plots depict the temporal progression of different groups of the population as seen in the figure. These results were generated using parameter values τ=13.75 and σ=10, which reveal periodic fluctuations before reaching stability. Each trajectory highlights the transient dynamics and steady-state approach of the model. Similarly, Figure 5 displays the simulation results for the temporal dynamics of each state variable over 1200 days with parameter settings τ=14.5 and σ=10. Compared to the previous configuration, these results show prolonged oscillations before stabilizing. The figure captures the system's response and adaptive transition to steady-state levels across the population groups. This comparative analysis underscores the model's sensitivity to parameter changes in achieving equilibrium.
Figure 6 is also the dynamics of each compartment as a function of varying delay τ with fixed delay σ=10 in a similar fashion. The time evolution of a specific compartment in the tuberculosis (TB) model is shown in each subplot over 1200 days for the delay parameter values of τ = 12.25, 12.5, 12.75, and 13.0 while a constant delay value of σ=10 is maintained. The trajectories show oscillatory behavior due to the model's transient dynamics and the oscillations grow in amplitude and frequency as τ increases before they reach equilibrium. In our results, we show that the time evolution of the stabilization within each compartment is sensitive to variations in the delay of the latent period. As in Figure 7, dynamics of all compartments are shown by varying delay σ with a fixed delay τ=13.5. The time evolution of each compartment in the TB model over 1200 days is shown in each subplot using the same latent period delay τ=13.5 but with varying lapse of time before moxifloxacin delivery (σ=10,12,14,16). Patterns of oscillations in the trajectory of each compartment's oscillatory pattern are noted, demonstrating the transient dynamics dependent on changes in σ, where oscillation amplitude and frequency alternate and settle toward an equilibrium. These results stress how changes to the quarantine-to-disability delay parameter affect the stabilization of each compartment over time.
In Figure 8, the three-dimensional dynamics of the TB compartments over time are depicted, highlighting the effect of varying the delay parameter σ. This parameter represents the time required for quarantined individuals to potentially develop disabilities. The figure illustrates how changes in σ influence the progression within each compartment under a fixed delay τ. This figure keeps the latent period delay τ constant, enabling us to isolate the effect of σ on population compartments. Figure 8a shows that as σ increases, there is a slower progression of individuals from the quarantined (Q) to the disabled (D) compartment. This change preserves more individuals in the infectious compartment longer, since fewer individuals transition into quarantine. As a result, the susceptible population declines more rapidly due to the prolonged infectious contact in the population. The vaccinated group in Figure 8b shows minimal fluctuation, as vaccination primarily affects the susceptible compartment. However, with increased σ, a minor reduction in the vaccinated group suggests that the burden on health resources may slightly increase, as more individuals remain infectious or require quarantine. The latent population in the figure 8c remains relatively unaffected by variations in σ, given that this delay primarily influences the infectious-to-quarantine transition rather than the latent-to-infectious progression. Similarly, the dynomic of infected compartment in Figure 8d shows the higher values of σ lead to an increase in the infectious population as individuals delay their progression to quarantine, thus maintaining the infection potential in the population. This dynamic illustrates the physical significance of the quarantine delay, as slower movement out of the infectious state prolongs the period individuals can contribute to the infection rate. As σ increases, the quarantine compartment's 8e growth slows, resulting in fewer individuals moving to the disabled state (D). This trend highlights how adjusting the delay in quarantine impacts both the infectious duration and the accumulation of disabilities from severe TB. The recovered compartments 8f show a delay in recovery and disability onset, respectively, with higher σ values, which implies that greater quarantine delays allow for prolonged infectious stages, indirectly impacting the number of individuals progressing to recovery or disability. In subplot Figure 10e, the impact of varying the delay σ on the disabled population is illustrated. A longer delay σ results in a slower transition of individuals into the disabled state, which subsequently reduces the peak and overall number of disabled individuals. This suggests that prolonging the delay σ could potentially lessen the burden of disability by preventing rapid progression to a disabled state. The results emphasize the importance of controlling σ (representing quarantine or isolation delays) to mitigate long-term disability impacts.
Figure 9 illustrates the impact of varying the latent period delay τ on the TB compartments, with a fixed quarantine-to-disability delay σ. This figure illustrates the influence of the incubation period, showing how delayed progression from latent to infectious states affects the overall dynamics of TB transmission. As τ increases, the susceptible population in the Figure 9a decreases more gradually. A longer latent period implies that fewer individuals transition quickly to an infectious state, reducing the rate of infection spread. Similarly the vaccinated compartment 9b sees limited changes, indicating that the primary dynamics of interest are among the susceptible, latent, and infectious groups. A longer latent period τ significantly raises the latent population 9c, as individuals spend more time incubating the TB bacterium. This change physically signifies a delay in active TB emergence, reflecting the incubation period's buffering effect on the immediate spread of infection. The infectious population in the Figure 9d decreases as τ increases, due to the slower transition from the latent to infectious stages. This effect highlights how prolonged incubation periods can serve as a natural intervention, reducing the immediate infectious burden within a population. The Quarantined and Recovered in Figures 9e and 9f with an extended incubation period τ show that the delayed onset of infectious cases reduces the influx into quarantine and subsequently impacts the recovered and disabled compartments. This effect is critical in understanding TB dynamics, as a prolonged latent phase could ease the strain on quarantine resources and reduce disability rates. In subplot 10f, the delay τ is varied while keeping σ constant. Similar to the effects observed with σ, increasing the delay τ also reduces the accumulation of disabled individuals over time. Since τ relates to the progression of infection, a longer delay implies a slower transition from infection to disability. This observation underscores the importance of intervention measures that manage the infection-to-disability progression timeline, as longer delays in disease progression correlate with a lower ultimate disability burden. Together, subplots 10e and 10f highlight how both σ and τ significantly influence the disabled population's dynamics, providing insights into effective delay-based interventions to reduce the long-term disability impact of TB.
In this study, we present a comprehensive mathematical model for TB transmission in Saudi Arabia, incorporating key epidemiological factors such as vaccination, latency delays, and disability progression. The findings offer significant insights into optimizing TB control strategies. The results highlight that an extended latency period (τ) reduces the infectious population, effectively slowing the transmission rate. This suggests that early detection and treatment of latent TB cases should be a priority to prevent disease progression. Additionally, the inclusion of a disability compartment in the model underscores the long-term impact of TB on public health, emphasizing the need for timely quarantine and treatment interventions. From a public health perspective, the model provides actionable recommendations for TB control programs in Saudi Arabia. First, optimizing quarantine measures can mitigate not only the spread of TB but also the risk of long-term disability among affected individuals. Implementing timely treatment protocols, particularly for individuals transitioning from quarantine to disability, can significantly reduce the socioeconomic burden associated with TB-related disabilities. Second, the findings advocate for enhanced TB surveillance systems that integrate real-time data collection and predictive modeling to improve decision-making for intervention strategies. Third, given that vaccination plays a role in controlling TB incidence, policies should focus on increasing vaccination coverage, particularly in high-risk populations, to achieve long-term disease reduction. This research provides policymakers with a quantitative framework for assessing intervention strategies and improving TB control efforts in Saudi Arabia. By translating mathematical insights into actionable public health measures, this model contributes to the development of targeted policies that enhance disease management, reduce transmission rates, and minimize TB-related disabilities. Researchers should integrate additional socio-economic and demographic factors, along with a bifurcation study and optimal control analysis, to further refine intervention strategies and support the country's efforts toward TB eradication.
Kamel Guedri: Writing-original draft, Formal analysis, Data curation, Supervision; Rahat Zarin: Writing-original draft, Software, Investigation, Conceptualization; Mowffaq Oreijah: Funding acquisition, Project administration, Methodology, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The authors extend their appreciation to the King Salman center For Disability Research for funding this work through Research Group no KSRG-2024-308.
The authors extend their appreciation to the King Salman center For Disability Research for funding this work through Research Group no KSRG-2024-308.
All authors declare no conflicts of interest in this paper.
[1] | Global Tuberculosis Report 2022, World Health Organization, 2022. Available from: https://www.who.int/teams/global-tuberculosis-programme/tb-reports |
[2] | Updated TB Treatment Guidelines, World Health Organization, 2023. Available from: https://www.who.int/teams/global-tuberculosis-programme/tb-guidelines |
[3] | Guidelines for the Treatment of Latent Tuberculosis Infection, Centers for Disease Control and Prevention of America, 2016. Available from: https://www.cdc.gov/tb/topic/treatment/ltbi.htm |
[4] |
S. Jitsinchayakul, R. Zarin, A. Khan, A. Yusuf, G. Zaman, U. W. Humphries, et al., Fractional modeling of COVID-19 epidemic model with harmonic mean type incidence rate, Open Phys., 19 (2021), 693–709. https://doi.org/10.1515/phys-2021-0062 doi: 10.1515/phys-2021-0062
![]() |
[5] |
R. Zarin, Numerical study of a nonlinear COVID-19 pandemic model by finite difference and meshless methods, Partial Differ. Equ. Appl. Math., 6 (2022), 100460. https://doi.org/10.1016/j.padiff.2022.100460 doi: 10.1016/j.padiff.2022.100460
![]() |
[6] |
Z. Raizah, R. Zarin, Advancing COVID-19 understanding: Simulating omicron variant spread using fractional-order models and haar wavelet collocation, Mathematics, 11 (2023), 1925. https://doi.org/10.3390/math11081925 doi: 10.3390/math11081925
![]() |
[7] |
R. Zarin, U. W. Humphries, T. Saleewong, Advanced mathematical modeling of hepatitis B transmission dynamics with and without diffusion effect using real data from Thailand, Eur. Phys. J. Plus, 139 (2024), 385. https://doi.org/10.1140/epjp/s13360-024-05154-7 doi: 10.1140/epjp/s13360-024-05154-7
![]() |
[8] |
K. Guedri, R. Zarin, A. Zeb, B. M. Makhdoum, H. A. Niyazi, A. Khan, A numerical study of HIV/AIDS transmission dynamics and the onset of long-term disability in chronic infection, Eur. Phys. J. Plus, 139 (2024), 1127. https://doi.org/10.1140/epjp/s13360-024-05881-x doi: 10.1140/epjp/s13360-024-05881-x
![]() |
[9] |
K. Guedri, R. Zarin, M. Oreijah, S. K. Alharbi, H. A. E. W. Khalifa, Artificial neural network-driven modeling of Ebola transmission dynamics with delays and disability outcomes, Comput. Biol. Chem., 115 (2025), 108350. https://doi.org/10.1016/j.compbiolchem.2025.108350 doi: 10.1016/j.compbiolchem.2025.108350
![]() |
[10] |
R. Zarin, Artificial neural network-based approach for simulating influenza dynamics: A nonlinear SVEIR model with spatial diffusion, Eng. Anal. Bound. Elem., 176 (2025), 106230. https://doi.org/10.1016/j.enganabound.2025.106230 doi: 10.1016/j.enganabound.2025.106230
![]() |
[11] |
M. Alqhtani, K. M. Saad, R. Zarin, A. Khan, W. M. Hamanah, Qualitative behavior of a highly non-linear Cutaneous Leishmania epidemic model under convex incidence rate with real data, Math. Biosci. Eng., 21 (2024), 2084–2120. http://doi.org/10.3934/mbe.2024092 doi: 10.3934/mbe.2024092
![]() |
[12] |
S. Alqahtani, A. Kashkary, A. Asiri, H. Kamal, J. Binongo, K. Castro, et al., Impact of mobile teams on tuberculosis treatment outcomes, Riyadh Region, Kingdom of Saudi Arabia, 2013–2015, J. Epidemiol. Glob. Health, 7 (2017), S29–S33. https://doi.org/10.1016/j.jegh.2017.09.005 doi: 10.1016/j.jegh.2017.09.005
![]() |
[13] |
O. Nave, Y. Shor, R. Bar, E. E. Segal, M. Sigron, A new treatment for breast cancer using a combination of two drugs: AZD9496 and palbociclib, Sci. Rep., 14 (2024), 1307. https://doi.org/10.1038/s41598-023-48305-z doi: 10.1038/s41598-023-48305-z
![]() |
[14] |
N. Zhang, X. Wang, W. Li, Stability for multi-linked stochastic delayed complex networks with stochastic hybrid impulses by Dupire Itô's formula, Nonlinear Anal. Hybrid Syst., 45 (2022), 101200. https://doi.org/10.1016/j.nahs.2022.101200 doi: 10.1016/j.nahs.2022.101200
![]() |
[15] |
N. Zhang, Z. Wang, J. H. Park, W. Li, Semi-global synchronization of stochastic mixed time-delay systems with Lévy Noise under aperiodic intermittent delayed sampled-data control, Automatica, 171 (2025), 111963. https://doi.org/10.1016/j.automatica.2024.111963 doi: 10.1016/j.automatica.2024.111963
![]() |
[16] |
N. Zhang, S. Huang, L. Ning, W. Li, Semi-global sampling control for semi-Markov jump systems with distributed delay, IEEE Trans. Autom. Sci. Eng., 21 (2024), 3603–3614. https://doi.org/10.1109/TASE.2023.3282053 doi: 10.1109/TASE.2023.3282053
![]() |
[17] |
R. G. White, G. P. Garnett, Mathematical modelling of the epidemiology of tuberculosis, Infect. Dis. Clin., 24 (2010), 825–841. https://doi.org/10.1007/978-1-4419-6064-1 doi: 10.1007/978-1-4419-6064-1
![]() |
[18] | H. T. Waaler, A. Geser, S. Andersen, The use of mathematical models in the study of the epidemiology of tuberculosis, Am. J. Public Health, 52 (1962), 1002–1013. |
[19] | Y. Li, J. Wang, Y. Zhao, A dynamical model of tuberculosis with media impact, J. Epidemiol., 17 (2022), 235–243. |
[20] | A. Das, S. Rana, A. Kumar, Media impact on the spread of tuberculosis: A mathematical model, Infect. Dis. Model., 5 (2020), 327–338. |
[21] | C. Bhunu, W. Garira, Z. Mukandavire, G. Magombedze, Modeling the effects of vaccination on the transmission dynamics of tuberculosis, Infect. Dis. Model., 3 (2008), 190–200. |
[22] | D. Okuonghae, P. L. Omosigho, Mathematical modeling of tuberculosis: vaccination and drug resistance, Math. Biosci., 220 (2011), 21–36. |
[23] | C. Castillo-Chavez, Z. Feng, Mathematical models for the transmission dynamics of tuberculosis, Math. Biosci., 151 (1997), 135–154. |
[24] |
C. Dye, B. G. Williams, Criteria for the control of drug-resistant tuberculosis, Proc. Natl. Acad. Sci. USA, 97 (2000), 8180–8185. https://doi.org/10.1073/pnas.140102797 doi: 10.1073/pnas.140102797
![]() |
[25] |
C. Ozcaglar, A. Shabbeer, S. L. Vandenberg, B. Yener, K. P. Bennett, Epidemiological models of Mycobacterium tuberculosis complex infections, Math. Biosci., 236 (2012), 77–96. https://doi.org/10.1016/j.mbs.2012.02.003 doi: 10.1016/j.mbs.2012.02.003
![]() |
[26] | Y. Wang, B. Xu, Z. Zhang, Epidemiological impact and mathematical analysis of drug-resistant tuberculosis, J. Infect. Dis., 41 (2023), 112–124. |
[27] |
P. Duve, S. Charles, J. Munyakazi, R. Lühken, P. Witbooi, A mathematical model for malaria disease dynamics with vaccination and infected immigrants, Math. Biosci. Eng., 21 (2024), 1082–1109. http://doi.org/10.3934/mbe.2024045 doi: 10.3934/mbe.2024045
![]() |
[28] | W. Walter, Ordinary differential equations, Berlin: Springer Science & Business Media, 2013. |
[29] | M. W. Hirsch, S. Smale, R. L. Devaney, Discrete dynamical systems, In: Differential equations, eynamical systems, and an introduction to chaos, New York: Academic Press, 2013,329–359. https://doi.org/10.1016/b978-0-12-382010-5.00015-4 |
[30] |
P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6 doi: 10.1016/S0025-5564(02)00108-6
![]() |
[31] |
Y. Li, X. Liu, Y. Yuan, J. Li, L. Wang, Global analysis of tuberculosis dynamical model and optimal control strategies based on case data in the United States, Appl. Math. Comput., 422 (2022), 126983. https://doi.org/10.1016/j.amc.2022.126983 doi: 10.1016/j.amc.2022.126983
![]() |
[32] |
Y. Wu, M. Huang, X. Wang, Y. Li, L. Jiang, Y. Yuan, The prevention and control of tuberculosis: an analysis based on a tuberculosis dynamic model derived from the cases of Americans, BMC Public Health, 20 (2020), 1173. https://doi.org/10.1186/s12889-020-09260-w doi: 10.1186/s12889-020-09260-w
![]() |
[33] |
J. Zhang, Y. Li, X. Zhang, Mathematical modeling of tuberculosis data of China, J. Theor. Biol., 365 (2015), 159–163. https://doi.org/10.1016/j.jtbi.2014.10.019 doi: 10.1016/j.jtbi.2014.10.019
![]() |
[34] |
Y. Ma, C. R. Horsburgh Jr, L. F. White, H. E. Jenkins, Quantifying TB transmission: a systematic review of reproduction number and serial interval estimates for tuberculosis, Epidemiol. Infect., 146 (2018), 1478–1494. https://doi.org/10.1017/S0950268818001760 doi: 10.1017/S0950268818001760
![]() |
[35] | Incidence of tuberculosis in Saudi Arabia, World Bank, 2023. Available from: https://data.worldbank.org/indicator/SH.TBS.INCD?locations = SA |
[36] | J. C. Butcher, Numerical methods for ordinary differential equations, Hoboken: Wiley, 2016. |
[37] | E. Hairer, S. P. Nørsett, G. Wanner, Solving ordinary differential equations I: nonstiff problems, 2 Eds., Berlin: Springer, 1993. |
[38] |
L. F. Shampine, S. Thompson, Solving DDEs in MATLAB, Appl. Numer. Math., 37 (2001), 441–458. https://doi.org/10.1016/S0168-9274(00)00055-6 doi: 10.1016/S0168-9274(00)00055-6
![]() |
[39] | U. M. Ascher, L. R. Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, Philadelphia: SIAM, 1998. |
[40] | A. Bellen, M. Zennaro, Numerical methods for delay differential equations, Oxford: Oxford University Press, 2003. |
Parameters | Description | Values | Source |
Λ | Recruitment rate of susceptible individuals | 1000/year | Assumed |
μ | Natural death rate | 0.05/year | Literature |
β | Transmission rate from infectious individuals to susceptibles | 1.4137×10−04year | Estimated |
ρ | Vaccination rate of susceptible individuals | 0.65/year | Fitted |
ε | Rate of waning immunity for vaccinated individuals | 0.55/year | Estimated |
φ | Rate of waning immunity for recovered individuals | 0.05/year | Assumed |
δ | Rate of latent infection clearance | 0.001/year | Estimated |
θ | Progression rate from latent to infectious TB (with delay τ) | 0.45/year | Estimated |
α | Rate at which infectious individuals enter quarantine or treatment | 0.1/year | Estimated |
η1 | Recovery rate for infectious individuals | 0.35/year | Literature |
ω | Rate of progression from infectious to disabled due to complications | 0.1/year | Estimated |
d2 | Disease-induced death rate for infectious and quarantined individuals | 0.035/year | Estimated |
η2 | Recovery rate for quarantined individuals | 0.07/year | Literature |
κ | Rate of progression from quarantine to disability (with delay σ) | 0.01/year | Estimated |
γ | Recovery rate for disabled individuals | 0.1/year | Estimated |
τ | Delay in progression from latent to infectious TB | [0−15]/day | Fitted |
σ | Delay in progression from quarantine to disability | [0−12]/day | Fitted |
Parameters | Description | Values | Source |
Λ | Recruitment rate of susceptible individuals | 1000/year | Assumed |
μ | Natural death rate | 0.05/year | Literature |
β | Transmission rate from infectious individuals to susceptibles | 1.4137×10−04year | Estimated |
ρ | Vaccination rate of susceptible individuals | 0.65/year | Fitted |
ε | Rate of waning immunity for vaccinated individuals | 0.55/year | Estimated |
φ | Rate of waning immunity for recovered individuals | 0.05/year | Assumed |
δ | Rate of latent infection clearance | 0.001/year | Estimated |
θ | Progression rate from latent to infectious TB (with delay τ) | 0.45/year | Estimated |
α | Rate at which infectious individuals enter quarantine or treatment | 0.1/year | Estimated |
η1 | Recovery rate for infectious individuals | 0.35/year | Literature |
ω | Rate of progression from infectious to disabled due to complications | 0.1/year | Estimated |
d2 | Disease-induced death rate for infectious and quarantined individuals | 0.035/year | Estimated |
η2 | Recovery rate for quarantined individuals | 0.07/year | Literature |
κ | Rate of progression from quarantine to disability (with delay σ) | 0.01/year | Estimated |
γ | Recovery rate for disabled individuals | 0.1/year | Estimated |
τ | Delay in progression from latent to infectious TB | [0−15]/day | Fitted |
σ | Delay in progression from quarantine to disability | [0−12]/day | Fitted |