Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Evaluating the impact of vaccination and progression delays on tuberculosis dynamics with disability outcomes: A case study in Saudi Arabia

  • 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

    Related Papers:

    [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.

    Figure 1.  Flow diagram of SVLIQRD model.

    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 t0, it is essential to establish their positivity throughout the time domain.

    Lemma 3.1. Let URn 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)×URn. Suppose that whenever yi=0 and all other yk>0 for ki, 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=1x1i({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(ˉτ))URn10.

    ● 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(ˉτ),,xi1(ˉτ))URi1>0 and (xi(ˉτ),,xn(ˉτ))URni+10. 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 x0U, 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)¯BNU. This guarantees that the solution remains positive, uniquely determined, and globally well-defined for all t0. 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 t0.

    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δLd2Id2Qγ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μtN(0)Λ(eμt1).

    Rearranging, we derive:

    N(t)Λμ+[N(0)Λμ]eμt.

    For all t0, 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 R01.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:

    V1=[1μ+δ+θ000θ(μ+δ+θ)(μ+d2+η1+α+ω)1μ+d2+η1+α+ω000α(μ+d2+η1+α+ω)(μ+d2+η2+κ)1μ+d2+η2+κ00ω(μ+d2+η1+α+ω)(μ+γ)κ(μ+d2+η2+κ)(μ+γ)1μ+γ].

    Multiplying F by V1, we obtain:

    FV1=[Λβθ(μ+δ+θ)(μ+ρερμ+ε)(α+μ+d2+η1+ω)Λβ(μ+ρερμ+ε)(α+μ+d2+η1+ω)00000000000000].

    The dominant eigenvalue (or spectral radius) of FV1 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 I0:

    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=βSIμ+δ+θ,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:

    R2R2ρμ+ρ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:

    R4R4θ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:

    R5R5α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:

    R6R6η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:

    R6R6η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βS0φ0ρ(μ+ε)00000βI0(μ+δ)θeζτβS00000θ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βS0φρ(μ+ε)0000βI0(μ+δ)θeζτβS0000θeζτ(μ+d2+η1+α+ω)00000α(μ+d2+η2+κeζσ)0000η1η2(μ+φ)].

    Using elementary row operations,

    Step 1: Eliminate ρ in row 2, column 1: R2R2ρβI+μ+ρR1.

    Step 2: Eliminate βI in row 3, column 1: R3R3βIβI+μ+ρR1.

    Step 3: Eliminate θeζτ in row 4, column 3: R4R4θeζτμ+δ+θeζτR3.

    Step 4: Eliminate α in row 5, column 4: R5R5αa44R4.

    Step 5: Eliminate η1 in row 6, column 4: R6R6η1a44R4.

    Step 6: Eliminate η2 in row 6, column 5: R6R6η2a55R5.

    Step 7: Eliminate the remaining non-zero entry in row 3, column 2: R3R3R3[2]R2[2]R2.

    We obtained the following upper triangular matrix:

    Jupper=[(βI+μ+ρ)ε0βS0φ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+β2ISβ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]=(μ+φ)η1R4[6]a44+η2R5[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ι)2nι=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.

    Figure 2.  Temporal evolution of yearly tuberculosis cases and reported deaths in the Saudi Arabia from 2000 to 2023 [35], along with the best-fitted curve derived from the proposed model.
    Table 1.  Notation and values of parameters in the tuberculosis model.
    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×1004year 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 [015]/day Fitted
    σ Delay in progression from quarantine to disability [012]/day Fitted

     | Show Table
    DownLoad: CSV

    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.

    Figure 3.  3D sensitivity analysis profiles showing the impact of various parameter pairs (μ, β, ω, ϵ, α, η1, ρ) on the basic reproduction number R0 for the proposed model. Each subplot illustrates how R0 responds to changes in two parameters simultaneously.

    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(tn1),tn1)). (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(tn1),tn1)). (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 tn1.

    Let's apply this ABM scheme to each equation. For simplicity, denote:

    SnS(tn), VnV(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(tn1)),
    Vp,n+1=Vn+h2(3fV(tn)fV(tn1)),
    Lp,n+1=Ln+h2(3fL(tn)fL(tn1)),
    Ip,n+1=In+h2(3fI(tn)fI(tn1)),
    Qp,n+1=Qn+h2(3fQ(tn)fQ(tn1)),
    Rp,n+1=Rn+h2(3fR(tn)fR(tn1)),
    Dp,n+1=Dn+h2(3fD(tn)fD(tn1)).

    Corrector Step

    Sn+1=Sn+h12(5fS(tn+1)+8fS(tn)fS(tn1)),
    Vn+1=Vn+h12(5fV(tn+1)+8fV(tn)fV(tn1)),
    Ln+1=Ln+h12(5fL(tn+1)+8fL(tn)fL(tn1)),
    In+1=In+h12(5fI(tn+1)+8fI(tn)fI(tn1)),
    Qn+1=Qn+h12(5fQ(tn+1)+8fQ(tn)fQ(tn1)),
    Rn+1=Rn+h12(5fR(tn+1)+8fR(tn)fR(tn1)),
    Dn+1=Dn+h12(5fD(tn+1)+8fD(tn)fD(tn1)).

    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 4.  The plots represent the simulation results of each state variable at τ=13.75 and σ=10.
    Figure 5.  The plots represent the simulation results of each state variable at τ=14.5 and σ=10.

    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.

    Figure 6.  The plots illustrate the temporal dynamics under the influence of delay τ, with a fixed delay parameter σ.
    Figure 7.  The plots illustrate the temporal dynamics under the influence of delay σ, with a fixed delay parameter τ.

    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 8.  Three-dimensional representations of state variables over time, demonstrating the effect of varying delay σ on each compartment while holding τ constant at 13.5. Each plot highlights the temporal dynamics under the influence of delay σ, with a fixed delay parameter τ, providing insights into the time-delay effects on population compartments within the model.

    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.

    Figure 9.  Three-dimensional representations of state variables over time, demonstrating the effect of varying delay τ on each compartment while holding σ constant at 8. Each plot illustrates the temporal dynamics under the influence of delay τ, with a fixed delay parameter σ, providing insights into the time-delay effects on population compartments within the model.
    Figure 10.  Dynamics of the Disable compartment.

    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.
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(417) PDF downloads(53) Cited by(0)

Figures and Tables

Figures(10)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog