Tuberculosis (TB), caused by the bacterium Mycobacterium tuberculosis, remains one of the most persistent and devastating infectious diseases worldwide. Despite decades of global control efforts, TB continues to represent a major public health burden, particularly in low- and middle-income countries where healthcare infrastructure is limited and diagnostic delay is common. The airborne nature of TB transmission, mediated through infectious droplets expelled by symptomatic individuals, facilitates rapid dissemination in densely populated environments and contributes to sustained endemicity in many regions [1,2].

From a clinical perspective, TB primarily targets the pulmonary system; however, its extrapulmonary manifestations-including involvement of the central nervous system, lymphatic system, skeletal structures, and other organs-introduce significant diagnostic complexity and therapeutic challenges [3]. These extrapulmonary forms complicate the classical Susceptible–Infected–Recovered (SIR) modeling framework, as they may involve different transmission pathways, latency structures, and heterogeneous progression rates that are not easily captured within a simple three-compartment structure. Nevertheless, the present study adopted the SIR framework due to its analytical tractability and its ability to capture the essential transmission dynamics, which was sufficient for investigating stability properties and constructing dynamically consistent numerical schemes.

The standard treatment regimen requires a prolonged course of combination antibiotics lasting at least six months, which often results in poor patient adherence and incomplete treatment. This issue has contributed to the emergence and global spread of multidrug-resistant (MDR-TB) and extensively drug-resistant (XDR-TB) strains, further complicating disease management and significantly increasing mortality risk [4,5].

According to the WHO, TB remains among the top ten causes of death globally, with approximately 10.8 million new cases and 1.25 million deaths reported in 2023 alone, the majority occurring in resource-limited settings [6]. These figures highlight not only the persistent nature of TB transmission but also the urgent need for advanced analytical frameworks capable of capturing its complex epidemiological behavior and supporting effective control strategies.

Mathematical modeling has become a fundamental tool in infectious disease epidemiology, providing a structured and quantitative framework for understanding transmission dynamics, evaluating intervention strategies, and predicting outbreak behavior. Classical compartmental models, particularly the SIR framework, have played a central role in capturing disease progression and estimating key epidemiological thresholds such as the basic reproduction number [7,8].

Despite its wide applicability, the SIR model has several limitations when applied to TB dynamics. In particular, it does not explicitly account for latent infection stages, reinfection mechanisms, or time delays associated with disease progression and treatment response. These limitations may restrict its ability to fully represent the complex natural history of TB.

In addition, stability analysis of equilibrium states offers critical insight into the long-term behavior of disease systems, including conditions for persistence or eradication.

However, the nonlinear nature of epidemic models often renders analytical solutions intractable, necessitating the use of numerical approximation techniques. Standard numerical methods such as the forward Euler and classical Runge-Kutta schemes are widely applied due to their simplicity and computational efficiency. Nevertheless, these methods may fail to preserve essential qualitative properties of epidemiological systems, such as positivity, boundedness, and dynamical consistency, particularly when larger time step sizes are used [9]. Such deficiencies can lead to nonphysical solutions and unreliable long-term predictions.

To address these limitations, the Nonstandard Finite Difference (NSFD) methodology, introduced by Mickens, has emerged as a powerful alternative numerical framework. The NSFD approach is specifically designed to preserve the essential dynamical features of the continuous system at the discrete level, ensuring that key biological properties - such as positivity of populations, stability of equilibria, and invariance of solution domains-are maintained for all admissible step sizes [10, 11]. This property makes NSFD schemes particularly suitable for epidemiological modeling, where biological realism is essential.

Recent literature has demonstrated the superiority of NSFD methods in modeling infectious diseases. Applications include typhoid fever dynamics, tuberculosis transmission modeling, and comparative studies of emerging infectious diseases such as monkeypox [12-14]. These studies consistently show that NSFD schemes outperform classical numerical methods by preserving qualitative dynamics and producing stable, biologically meaningful simulations even under coarse discretization.

Furthermore, recent advancements in computational epidemiology have introduced hybrid modeling approaches that combine mechanistic models with machine learning techniques. In particular, dynamics-informed neural networks have been applied to model COVID-19 transmission, demonstrating improved predictive accuracy by integrating differential equation structures with data-driven learning [15]. While dynamics-informed neural networks provide a data-driven and flexible framework for capturing complex epidemic patterns, they are fundamentally different from NSFD schemes. The NSFD method is a structure-preserving numerical technique designed to maintain the qualitative properties of deterministic models, whereas neural networks focus on predictive accuracy and data assimilation. Therefore, we do not claim the superiority of NSFD over machine learning approaches, but rather emphasize that NSFD serves as a complementary tool for ensuring reliable and biologically consistent numerical simulations.

In addition, recent developments in mathematical epidemiology have highlighted the importance of delay differential models and advanced stability analysis techniques. For example, delay-induced stability transitions and Hopf bifurcation phenomena have been studied in various applied contexts, providing deeper insight into oscillatory behavior and long-term dynamics [16, 17]. These approaches offer valuable extensions to classical models and motivate future research directions beyond the scope of the present study.

Motivated by these challenges, we developed a nonlinear SIR-type mathematical model for tuberculosis transmission and proposed a dynamically consistent NSFD scheme for its numerical approximation. The proposed framework ensures the preservation of essential qualitative properties, including positivity, boundedness, and stability, at the discrete level. Moreover, we provide a reliable computational tool for analyzing TB dynamics which could be extended to other nonlinear epidemiological systems where classical numerical methods fail to maintain biological consistency.

The modeling approach adopted in this study was grounded in classical compartmental epidemiology, where the host population was partitioned into epidemiological classes based on disease status. The formulation was developed under standard assumptions of homogeneous mixing, constant recruitment, and time-invariant parameters, which allowed for tractable analytical and numerical investigation while preserving essential biological realism. Particular attention was given to ensuring that the resulting dynamical system satisfied fundamental epidemiological properties such as positivity, boundedness, and biological feasibility of solutions.

The mathematical framework was further structured to facilitate the derivation of key epidemiological quantities, including equilibrium points and the basic reproduction number , which served as threshold indicators for disease persistence or elimination. In addition, the system was formulated in a way that supported both qualitative analysis and efficient numerical implementation using structure-preserving computational methods.

To ensure clarity and reproducibility, all state variables, parameters, and modeling assumptions were explicitly defined in the subsequent subsection. The overall objective of this section was to establish a rigorous and self-consistent mathematical foundation for analyzing the transmission dynamics of tuberculosis and to provide a suitable framework for stability analysis and numerical simulation in the later sections of the paper.

Mathematical Model of Tuberculosis

To study the transmission dynamics of tuberculosis (TB), we considered a classical compartmental SIR epidemic model. The total population is partitioned into three mutually exclusive epidemiological classes: susceptible individuals , infectious individuals , and recovered individuals .

The susceptible class represents individuals who are at risk of contracting TB. The infectious class consists of individuals with active pulmonary tuberculosis who are capable of transmitting the disease. The recovered class includes individuals who have completed treatment and are assumed to have acquired partial or full immunity, depending on the biological assumption of the model.

The schematic representation of TB transmission dynamics is illustrated in Figure 1.

Figure 1: Schematic diagram of tuberculosis transmission dynamics.

Model Formulation

The transmission dynamics of tuberculosis are governed by the following system of nonlinear ordinary differential equations:

where represents the total population size.

The initial conditions are given by:

Parameter Description

The biological meaning of each parameter is defined in Table 1.

Table 1: Description of parameters in the tuberculosis transmission model

All parameters are assumed to be positive constants. The incidence term represents a standard incidence function under homogeneous mixing assumptions.

Model Reduction

Since the recovered class does not influence the dynamics of and directly, the system can be reduced to a two-dimensional subsystem for analytical purposes:

This reduced system is sufficient for studying key dynamical properties such as equilibrium points, stability analysis, and threshold behavior (basic reproduction number ), which will be investigated in the subsequent sections.

Equilibrium Points

Disease-Free Equilibrium (DFE)

The DFE represents a state where no infection is present in the population. By setting the infectious class to zero ( ) and solving the steadystate equations ( ) from System (1), the DFE, denoted by , is obtained as:

This equilibrium represents a healthy population with no infection present and serves as a baseline for disease invasion analysis.

Endemic Equilibrium (EE)

The EE occurs when the disease persists in the population . By setting and from System (2) and assuming , expressions for , and can be derived. Due to algebraic complexity, the detailed derivation is omitted here but is included in the full analysis. The existence of a unique endemic equilibrium is guaranteed when the basic reproduction number exceeds unity ( ).

Basic Reproduction Number

The basic reproduction number is a fundamental epidemiological threshold that determines whether the infection will die out or become endemic. As discussed by Ochwach and Okongo [18], the next-generation matrix method is a robust approach for deriving in complex epidemic models. It reflects the average number of secondary infections produced by a single infectious individual in a fully susceptible population.

We calculate using the next-generation matrix method. At the DFE, the matrices for new infection and transition are:

The basic reproduction number, which is the spectral radius of the next-generation matrix , is derived as:

If , the disease-free equilibrium is stable and the disease will die out. If , the infection will spread, potentially becoming endemic.

Forward Euler Finite Difference Scheme

First order differential equations can be solved using the Euler scheme, which is purely numerical. In order to solve the given SIR model, Euler's method is more appropriate [19]. The description of Euler method is given below.

Finally, we get the Euler scheme:

As shown in Table 2, the fixed parameter values used in System (1) were selected to investigate the numerical behavior and stability properties of the tuberculosis transmission model under different discretization schemes.

Table 2: Fixed values of parameters used in system (1)

It is evident that Euler's method, when applied to the DFE, yields positive numerical results only at very small step sizes. Figure 2 illustrates that increasing the step size leads to numerical instability and non-physical negative solutions. This failure occurs because the Euler scheme treats all terms on the right-hand side explicitly. Consequently, the discrete solution does not inherently enforce the invariant region (positivity and boundedness) of the continuous system. For large , the explicit linear approximations overshoot, causing the solution to leave the feasible domain and eventually diverge. Thus, we conclude that the stability and positivity of the Euler scheme cannot be maintained at higher step sizes.

Figure 2: Mathematical results of SIR epidemic system (1) gained by Euler pattern through (a) , (b) , (c) and (d) . All parameter values are given in Table 1.

Runge-Kutta Order 4 (RK-4)

Due to its precision and ease of use, the Runge-Kutta method is a widely utilized approach for solving ordinary differential equations (ODEs) [20]. Not all derivatives of functions are easy to compute, which is why the Taylor series method relies on calculating them [21]. Let , and , then the RK-4 scheme is represented as follows.

Step-1:

Step-2:

Step-3:

Step-4:

Final step:

Although the RK-4 method achieves higher accuracy than Euler for small step sizes due to its multi-stage implicit averaging, it remains fundamentally an explicit scheme. Figure 3 demonstrates that when the step size exceeds a certain threshold, the RK-4 method also becomes unstable and diverges. The reason is similar to Euler: the scheme does not enforce the model's invariant region (positivity and boundedness) and its stability region in the complex plane is limited. For stiff epidemiological systems like TB transmission, explicit RK-4 cannot guarantee biologically meaningful solutions for large .

Figure 3: Mathematical results of SIR epidemic system (2) gained by RK-4 structure through (a) , (b) , (c) and (d) . Table 1 contains the values for every other parameter

NSFD

Construction of the NSFD Scheme

Our primary goal in this section is to propose an NSFD scheme for this model. Mickens presents a novel finite difference method. Their dynamical analysis demonstrates that, regardless of the step size, the NSFD scheme maintains the stability properties of all prevailing equilibrium points [22]. The key difference between NSFD and classical methods lies in the implicit (or partially implicit) treatment of nonlinear terms and the use of nonlocal approximations. Specifically, the incidence term is approximated in a way that preserves the conservation laws and positivity of the continuous system. It is possible to construct the NSFD scheme for the continuous dynamical System (2) as shown below:

Notice the strategic placement of the index: the outflow terms from (the incidence and death) are treated implicitly, while the inflow is kept explicit. Similarly, the equation treats its loss terms implicitly. This implicit formulation prevents the denominator from becoming zero and guarantees that for any positive step size , the updated populations , and remain nonnegative. Moreover, the discrete system preserves the same equilibrium points as the continuous model.

The initial conditions and of the discrete NSFD SIR system (3) are likewise expected to be nonnegative. To obtain this in explicit form, the discrete NSFD system (3) can be reorganized:

In this explicit form, all right-hand side terms are positive when the previous iterates are positive. The denominators are always greater than 1 , which acts as a contraction mapping that prevents blow-up. This is fundamentally different from Euler or RK-4, where the new values are computed by adding potentially large multiples of derivatives, leading to overshoot and instability when is large. The NSFD denominator provides a natural "saturation" effect, keeping the solution within the feasible region for any step size .

Since all parameters in (3) are positive, it is proven that and for all and for any . This qualitative behavior of the NSFD scheme, including the preservation of positivity and stability for different step sizes, is illustrated in Figure 4.

Figure 4: Mathematical results of SIR epidemic system (2) obtained by NSFD structure through (a) , (b) , (c) and (d) . All remaining parameter values are given in Table 1

Local Stability Analysis for the Discrete NSFD Scheme

We consider the following functions to demonstrate that DFE and DEE points are locally asymptotically stable (LAS):

Theorem 1 If , then the DFE point of the NSFD structure (3) is LAS for all .

Proof Let us consider the Jacobian matrix:

We compute the partial derivatives:

Substituting these values into equation :

The eigenvalues are clearly shown by the aforementioned matrix:

Since both eigenvalues are strictly less than 1 in magnitude for any , the DFE is locally asymptotically stable for the discrete NSFD system whenever , regardless of the step size. This is in sharp contrast to explicit methods, whose stability depends critically on being sufficiently small.

This verifies that if , then the DFE point is locally asymptotically stable; conversely, if , it is not.

In the next section, we discuss the application of the Schur-Cohn criterion [23, 24], which plays a significant role in analyzing the local stability of the disease-endemic equilibrium point .

Lemma 2 (Schur-Cohn Criterion) The solutions of the equation satisfy for if and only if the conditions listed below are met:

,

,

, where and represent the determinant and trace of the Jacobian matrix, respectively.

Theorem 3 If , then the endemic equilibrium point of the NSFD structure (3) is LAS for all .

Proof Let us consider the Jacobian matrix at :

By substituting and , and then putting the DEE point into (5), we obtain:

where , and . Here .

The eigenvalues of are determined by the characteristic equation:

where

The following outcomes are derived:

Due to , we have:

It is clear that .

Through straightforward calculation, it is shown that if , then:

Thus, all three Schur-Cohn conditions are satisfied for any whenever . This proves that the NSFD scheme preserves the local asymptotic stability of the endemic equilibrium without any step-size restriction. Classical explicit methods cannot guarantee this property.

If for any step size , then the endemic equilibrium point is locally asymptotically stable according to the Schur-Cohn criterion presented in Lemma 1.

3.5.3 Global Stability Analysis of the Discrete NSFD System

Now, using the appropriate Lyapunov function, we can examine the global stability of disease-free and disease-endemic equilibria.

Theorem 4 If , then the DFE point of the discrete NSFD scheme (3) is globally asymptotically stable.

Proof Let us consider the following discrete Lyapunov function:

where . Using the first and second equations of System (3), we can derive (6). Now:

Using the Enatsu et al. criterion [25,26], , we obtain:

Certainly, if , then for all . This shows that is a monotonic decreasing sequence. Since , we have and ( . Consequently, we obtain and . It is clear that if , then . Thus, we conclude that is globally asymptotically stable.

Theorem 5 If , then the endemic equilibrium point of the NSFD system (3) is globally asymptotically stable for any step size .

Proof We apply the Enatsu et al. [26] criterion to establish necessary conditions for the global stability of . For this, the discrete Lyapunov function is defined as follows:

where and . Equation (7) is created by substituting the first and second equations into the Lyapunov function.

The global stability behavior of the NSFD scheme for different step sizes is illustrated in Figure 5.

Figure 5: Mathematical results of SIR system (2) achieved from NSFD structure through (a) , (b) , (c) , (d)

Now:

Let us denote , and . Then (8) can be expressed as:

Using the definition of , we obtain:

Thus, is a monotonic decreasing sequence. We can demonstrate that by using the same methods as in Theorem 3. As a result, for all , we obtain and . The proof is thus conclusive.

3.6 Comparative Numerical Performance

To further demonstrate the effectiveness of the NSFD scheme, we conducted a comparative analysis with the classical Euler and RK-4 methods. Table 3 presents the relative error values for all three methods at different time steps. It is evident that the NSFD scheme maintains accuracy and positivity even at larger step sizes, unlike the Euler and RK-4 schemes which diverge or yield biologically invalid (negative) values. The fundamental reason, as established above, is that the NSFD scheme employs a nonlocal implicit discretization of the nonlinear incidence term and treats loss terms implicitly. This approach enforces the invariant region of the continuous model at the discrete level, guaranteeing that the numerical solution remains bounded and positive for any step size . In contrast, Euler and RK-4 are purely explicit methods that approximate derivatives linearly; their stability regions are finite, so when exceeds the stability limit (approximately where is the largest eigenvalue of the system's Jacobian), the numerical solution becomes unstable and blows up.

As shown in Table 3, when increases from 0.1 to 1, both Euler and RK-4 become unstable or diverge, while the NSFD method remains stable and accurate. This confirms that the NSFD scheme provides robust numerical solutions suitable for long-term simulations of TB dynamics, whereas the other methods fail to retain essential properties under coarse discretization due to their explicit nature and lack of structure-preserving design.

Table 3: Comparison of Euler, RK-4, and NSFD methods in terms of error, stability, and positivity at different step sizes.

The findings of this study demonstrate that the NSFD scheme significantly outperforms classical numerical methods (Euler and RK-4) in preserving the essential qualitative properties of tuberculosis transmission dynamics. While Euler and RK-4 methods are computationally simple and widely used, they fail to maintain positivity and stability when larger time steps are employed a critical limitation for long-term epidemiological simulations where computational efficiency is important. The threshold parameter derived in this study provides a practical metric for public health intervention planning. However, the real-world applicability of the NSFD scheme's stability results must be interpreted in light of the model's underlying assumptions and limitations.

Impact of Modeling Limitations on NSFD Stability Results

The theoretical stability results proved for the NSFD scheme-namely, the preservation of positivity, boundedness, and asymptotic stability of equilibria for any step size -are mathematically rigorous under the assumptions of the proposed SIR model. However, each acknowledged limitation of the model affects how these stability results translate to real-world tuberculosis transmission. Below, we analyze each limitation separately.

Permanent immunity assumption: The model assumes that recovered individuals acquire permanent immunity and cannot be reinfected. In reality, TB reinfection is well-documented, particularly in high-burden settings. From the perspective of NSFD stability, this assumption has two opposing effects. First, it overestimates the effectiveness of recovery in reducing the susceptible population, which may lead to an underestimation of the true and an overestimation of the stability region of the disease-free equilibrium. Consequently, the NSFD scheme would predict disease elimination under conditions where, in reality, reinfection would sustain transmission. Second, however, the NSFD scheme's ability to preserve stability for any step size is a structural property of the discrete scheme itself, not contingent on the biological accuracy of the immunity assumption. Therefore, while the scheme remains dynamically consistent with the model, the model's deviation from biological reality means that the NSFD-predicted stability thresholds may not match real-world epidemiological thresholds. This limitation does not reduce the numerical reliability of NSFD but rather highlights the need to extend the model to include reinfection and waning immunity.

Homogeneous mixing assumption: The model assumes that all individuals mix uniformly and have equal contact rates. In real populations, contact patterns are heterogeneous, age-structured, and influenced by social networks, geography, and behavior. The NSFD stability analysis relies on the standard incidence term , which implicitly assumes homogeneous mixing. Under heterogeneous mixing, the effective reproduction number and the stability properties of the system can differ significantly, with certain subgroups (e.g., healthcare workers, densely populated households) acting as high-transmission cores that sustain endemicity even when the population-average is below unity. The NSFD scheme, as constructed, would still preserve positivity and stability for any step size when applied to a homogeneous model, but it cannot capture the stabilizing or destabilizing effects of heterogeneity. Therefore, the NSFD results are most applicable to well-mixed populations (e.g., small communities, schools, or institutional settings) and should be extended using age-structured or network-based models for heterogeneous populations.

Parameter values not representing all settings (e.g., Afghanistan): The parameter values used in numerical simulations (Table 1) were obtained from the literature and may not accurately reflect the epidemiological conditions of specific regions, particularly Afghanistan, where factors such as malnutrition, delayed diagnosis, limited healthcare access, high population density, and co-morbidities (e.g., HIV, diabetes) significantly alter transmission dynamics. From the perspective of NSFD stability, the scheme's mathematical properties (positivity, boundedness, step-size independence) are parameter-agnostic-they hold for any positive parameter values. However, the biological interpretation of the stability results is parameter-dependent. For example, if the true TB-induced death rate in Afghanistan is lower than the literature value used, or if the natural recovery rate is overestimated due to poor treatment adherence, the true may be substantially higher than our computed value. Consequently, the NSFD scheme would still produce stable, positive solutions, but the predicted equilibrium (disease-free vs. endemic) might be incorrect. The robustness of NSFD lies in its ability to reliably simulate the dynamics given the parameters, but accurate parameter estimation for specific regions remains a critical prerequisite for practical applicability.

In summary, the NSFD scheme preserves the mathematical stability properties of the model regardless of step size. However, the model's simplifications (permanent immunity, homogeneous mixing) and the use of generic parameter values mean that the NSFD-predicted stability thresholds and long-term dynamics may deviate from real-world TB transmission. These limitations do not diminish the NSFD scheme's numerical superiority over Euler and RK-4, but they do define its scope of applicability: the NSFD method is most reliable for comparative studies, sensitivity analyses, and simulations where the model structure is accepted as a reasonable approximation. For high-fidelity predictions in specific regions like Afghanistan, the NSFD framework should be extended with more realistic model components (see Future Directions) and calibrated with local data.

Future Directions

Building on the limitations discussed above, future work should extend the NSFD framework to address each modeling shortcoming while preserving dynamical consistency. First, to address the permanent immunity limitation, an SIRS model (incorporating waning immunity and reinfection) or a model with an explicit latent compartment (SEIR) should be developed, and dynamically consistent NSFD schemes for these extended models should be constructed. Second, to address homogeneous mixing, network-based or age-structured NSFD formulations should be explored, potentially using operator-splitting techniques to preserve stability. Third, to improve parameter relevance for specific settings like Afghanistan, future studies should undertake systematic parameter estimation using local epidemiological data, including contact surveys, treatment outcome studies, and mortality surveillance. Additional promising directions include:

Fractional-order derivatives for memory effects in TB incidence

Spatio-temporal dynamics for geographic spread of drug-resistant strains

Stochastic NSFD schemes for demographic variability in low-incidence settings

Optimal control strategies integrated with NSFD for intervention optimization

Delay differential equation models with NSFD to capture treatment and progression delays

The impact of treatment on TB patients is profound, particularly at the level of recovery; the greater the emphasis on therapeutic parameters, the more likely it is for TB patients to achieve successful recovery. In this study, we derived the basic reproduction number, which is essential for examining both the local and global stability of disease-free (DFE) and endemic equilibrium (DEE) points using discrete Lyapunov and monotonic criteria.

Although Euler and RK-4 methods provide basic numerical approximations, they exhibit limited convergence and fail to preserve positivity at larger step sizes, resulting in instability. To address these limitations, we developed a nonstandard finite difference (NSFD) scheme that is not only mathematically consistent but also biologically meaningful, as it preserves the essential features of the continuous TB model across all step sizes.

Unlike the Euler and RK-4 methods, which yield divergence and biologically invalid (negative) outcomes at coarse time steps, the NSFD scheme guarantees positivity, stability, and convergence regardless of step size. Furthermore, the NSFD method maintains the dynamical integrity of the TB model, ensuring reliable long-term simulations-a feature critical for accurate disease prediction.

Theoretical results confirm that the DFE and DEE points are locally and globally asymptotically stable under specific criteria, supported by discrete Lyapunov analysis and the Schur-Cohn condition. This advancement strengthens the model's applicability for real-world epidemiological planning and makes it suitable for policy simulations and long-term disease forecasting, especially in resource-limited settings. Finally, the proposed framework provides a solid foundation for future research, including extensions to fractional-order, fuzzy logic-based, and spatio-temporal epidemic models, where preserving biological realism and mathematical robustness remains vital.

The authors would like to express their sincere gratitude to the Department of Mathematics, Faculty of Science, Herat University, for providing the necessary academic environment and support to carry out this research.

No external funding was received for this study.

The authors declare that they have no competing interests.

Code availability

MATLAB codes used for numerical simulations are available from the corresponding author upon reasonable request.