A mathematical model for the study of latent tuberculosis under 3HP and 1HP regimens

1.   Introduction

Tuberculosis (TB) has a significant impact on public health worldwide. It is estimated that in 2019, approximately ten million people worldwide developed TB and 1.2 million died from the disease. As for treatment outcomes, in 2018, the percentage of treatment success was 85% for new cases ^[1].

Monotherapy with isoniazid for 6–9 months has been used for decades, and its efficacy in preventing progression to active TB is approximately 90%. But this overall efficacy has been affected by low adherence and completion rates due to its prolonged duration and risks of hepatoxicity ^[2,3].

Treatment of TB infection is called TB preventive therapy and helps prevent the development of the disease. If the infection is not treated, it can develop into an active infection, causing the person to become ill and transmit TB.

Rifapentine is in the class of drugs called rifamycins, which form the basis of the new short-acting TB preventive therapy drugs. When rifapentine is combined with isoniazid, which is another anti-TB drug, two possible uses are obtained: the 3HP regimen (given once a week for 12 weeks) and the 1HP regimen (given once a day for one month). The 3HP and 1HP regimens are shorter alternatives compared to the older standard of care called isoniazid preventive therapy, in which people need to take isoniazid daily for 6 up to 36 months ^[2,3,4]. Any therapy applied to cases with latent TB other than 3HP and 1HP for our work are called classical treatment.

Short-duration treatment regimens, such as 1HP and 3HP, are effective and safe and have higher completion rates than isoniazid monotherapy classical therapy. Short regimens with rifamycin use have a lower risk of hepatoxicity ^[2,3,4].

Models for TB treatment have been increasing over the last decade. For example, Delgado et al. ^[5] present a mathematical model with ordinary differential equations for the efficacy of TB treatment, taking into account the influence of HIV and diabetes, and a version of this study using fractional order derivatives can be found in ^[6]. Chong et al. presented a mathematical model to assess the impact of treatment of latent TB infection in the elderly and the use of current TB control strategies^[7]. Sulayman et al. ^[8] extended a mathematical model that studies the impact of a vaccine with a limited partial effect and other exogenous factors, such as re-infection between treated individuals, and proved that vaccination can be effective for TB control in a population. Odibat et al. ^[9] proposed a mathematical model of TB transmission dynamics incorporating first and second-line treatment.

Treatment of latent TB is essential to prevent future active cases that may spread the disease in society. Current latent TB therapies can have a high degree of toxicity and poor adherence, so 3HP and 1HP therapies would have a positive effect. This work aims to present a mathematical model for TB taking into account the impact of latent TB on disease transmission. For that, we take into account the cases that have latent TB and eliminate the bacteria naturally, the effect of classical therapy, and also the use of 3HP and 1HP therapies. By using fractional derivatives in the Caputo sense, we exploit the advantages of such a modeling technique concerning the memory effect ^[10], comparing the influence of the 3HP and 1HP latent treatments on a population with respect to the classical treatment applied therapies.

This paper is organized as follows: In Section 2, we present the theoretical background with the basic definitions used in the paper. In Section 3, we introduce the model and its basic properties, and we study the basic reproduction number. Section 4 is dedicated to numerical simulations. Section 5 closes the paper with conclusions.

2.   Theoretical background

For the sake of completeness, in this section, we present the mathematical definitions used in the paper ^[11,12,13]. In what follows, we assume that $\alpha \in \mathbb{R}_{+}$, $b > 0$, $f\in AC^{n}[a, b]$, and $n = [\alpha]$ ($AC^{n}$ is the set of functions with order derivative $n-1$ absolutely continuous, and $[\alpha]$ is the integer part of $\alpha$ (the integer part of $\alpha$ is the greatest relative integer less than or equal to $\alpha$)).

Definition 2.1. We define the left-sided and right-sided fractional integral Riemann-Louville for $f$: $\mathbb{R}_{+}\longrightarrow \mathbb{R}$ and $\alpha > 0$ as:

where $\Gamma(\cdot)$ is the Gamma function.

Remark 2.1. Let us define

Definition 2.2. The left-sided and right-sided Riemann-Liouville fractional derivatives are defined as:

Remark 2.2. Let us denote

Definition 2.3. The left-sided and right-sided fractional derivatives proposed by Caputo are given by:

Remark 2.3. Let us define

In the order-fractional derivatives, we find the memory effect, which is an important factor in epidemic modeling (see ^{[10,14,15,16,17]}).

3.   Mathematical modeling

The model is defined in terms of 6 compartments: susceptible $S$, exposed $E$, latent $L$, latent diagnosed, selected latent to use new treatments $L_{T}$, drug-sensitive $D$, and drug-resistant $R$. The latent compartment $L$ contains diagnosed and undiagnosed cases that use classical therapy elements. In order to need to enter treatment for latent TB, they need to be diagnosed and selected for treatment. $L_{T}$ represents the cases that are diagnosed and enter treatment.

The parameter $M$ represents the recruitment rate at which new susceptibles are entering into the dynamic, and $\mu$ is the natural death rate. The TB transmission rate is defined as:

where $\beta$ is the effective transmission rate, $\epsilon_{1}$ represents the modification parameter associated with the infectivity of the resistants, and $N$ is the total population

People are infected from the latent state with the transmission rate of $\lambda_{T}$ with $\epsilon_{2}$, which is the modification parameter related to infection in the latent state. The modification parameters $\epsilon_{1}$ and $\epsilon_{2}$ adapt the infectivity of the resistant ones, which is higher than in the drug-sensitive ones, and the latent ones can become infected compared to the susceptible ones. These parameters, by definition, take values greater than unity and can be obtained by performing a statistical study in the respective subpopulations. The parameter $v$ is the rate at which bacteria are acquired and lead patients to the active or latent stage. The rate $r_{L}$ represents patients with latent TB who are selected to use the 3HP and 1HP therapies. The expressions $(1-p_{R})\alpha_{R}$ and $(1-p_{E})(1-\theta)\psi D$ represent cases where the treatment was not effective and the patient becomes latent. The expression $(1-\phi)v$ represents the exposed cases that acquire the bacteria, but the bacteria are in a non-active state. Moreover, the expression $(1-p_{L})\sigma_{L}$ represents the cases where 3HP or 1HP were not successful, causing a return to the latent state. The parameter $\sigma$ represents the cases that do not develop active TB and become susceptible naturally (i.e., their immune system responds without treatment) and those who use treatment and achieve treatment success. The parameter $\omega$ represents the activation cycle of the bacteria. The parameter $p_{L}$ is associated with the latents that eliminate the bacteria.

The exposed that continue the bacterial cycle in a drug-sensitive manner are represented by the expression $(1-p_{E})\phi vE$, and those that do so in a drug-resistant manner are expressed as $p_{E}\phi vE$. Latents that activate the bacteria as sensitive to TB drugs are represented by the expression $(1-p_{E})\omega L$, and those that make them resistant by $p_{E}\omega L$. Those that enter the drug-sensitive compartment and develop resistance are expressed as $(1-\theta)p_{E}\psi D$. The modification parameter $t_{h}$ adapts TB deaths to resistant cases. Patients entering treatment in drug-sensitive and drug-resistant are represented by the $\psi$ and $\alpha_{R}$ parameters, respectively.

Figure 1 illustrates its flow diagram, and Tables 1 and 2 present the description of the variables and parameters of the model. The fractional derivative operator ${}^{c}\mathbb{D}_{t}{^\alpha}$ has time dimension ${\alpha}$, then the parameters have power dimension $\alpha$, except for the modification parameters.

This dimensional consistency can also be found in the works ^[16,17]. We assume that the modification parameters do not change dimension. The following model describes the transmission dynamics of TB, taking into account the use of 3HP and 1HP using fractional derivatives in the Caputo sense:

The initial conditions are:

3.1.   Stability analysis

Now, let us prove the existence and positivity of the solution of systems (3.2)–(3.7), and let's find the biologically feasible region.

3.1.1.   Existence and non-negativity of solutions

Let us denote

The following lemma and corollary are used in the proof of Theorem 3.1 and can be found in ^[18,19].

Lemma 3.1. (Generalized mean value theorem) Suppose that $f\in C[a, b]$ and ${}^{c}\mathbb{D}_{t}{^\alpha}f\in C[a, b]$, for $\alpha \in (0, 1]$. Then, $\forall t\in (a, b],$ with $a\leq \epsilon \leq t$, we have

Corollary 3.1. Consider that $f\in C[a, b]$ and ${}^{c}\mathbb{D}_{t}{^\alpha}f\in C[a, b]$, for $\alpha \in (0, 1]$. Then if

● ${}^{c}\mathbb{D}_{t}{^\alpha}f(t) \geq 0$, $\forall t\in (a, b)$, then $f(t)$ is non-decreasing for each $t\in[a, b]$,

● ${}^{c}\mathbb{D}_{t}{^\alpha}f(t) \leq 0$, $\forall t\in (a, b)$, then $f(t)$ is non-increasing for each $t\in[a, b]$.

Theorem 3.1. There is a unique solution

of models (3.2)–(3.7) for $t\geq 0$ and the solution remains in $\Omega$.

Proof. Using Theorem 3.1 and Remark 3.2 of ^[20], we obtain the solution of the initial value problems (3.2)–(3.7) exists and is unique at $(0, \infty)$.

Now, we prove the positivity of the solution of the systems (3.2)–(3.7). We obtain the positivity of the solution if we show that for each hyperplane bounding the nonnegative orthonormal, the vector field points to $\Omega$. From the models (3.2)–(3.7), we have:

Using the Corollary 3.1, we conclude that the solution remains in $\Omega$. □

3.1.2.   Biologically feasible region

Our next step is that we need to find the region where our models (3.2)–(3.7) make biological sense, which we call the biologically feasible region ($\Omega_{\alpha}$). In this case, we are working with compartments of humans, so our variables have to be positive or zero.

Lemma 3.2. The closed set

is positively invariant with respect to models (3.2)–(3.7).

Proof. For our total population, the fractional derivative in the Caputo sense is defined as:

and we have

The following definitions contribute to the demonstration of the result:

Definition 3.1. The Laplace transform of the Caputo fractional derivatives of the function $\phi(t)$ with order $\alpha > 0$ is defined as

Definition 3.2. The Laplace transform of the function $t^{\alpha_{1}-1}\mathbb{E}_{\alpha, \alpha_{1}}(\pm \lambda t^{\alpha})$ is defined as

where $\mathbb{E}_{\alpha, \alpha_{1}}$ is the two-parameters Mittag-Leffler function $\alpha, \alpha_{1} > 0$. Furthermore, the Mittag-Leffler function satisfies the following equation:

Applying the Laplace transform to (3.10), we have

which further gives

Using the Eqs (3.11)–(3.13) and we assumed that $(S(0), E(0), L(0), L_{T}(0), D(0), R(0))\in \mathbb{R}^{6}_{+}$, then

Using the asymptotic behavior of the Mittag-Leffler function, we can observe that

as $t \rightarrow \infty$.

The $\Omega_{\alpha}$ region is well established, and all solutions with initial values that belong to $\Omega_{\alpha}$ remain in $\Omega_{\alpha}$ for each time $t > 0$. □

3.2.   Basic reproduction number

In a population composed only of susceptible individuals, the average number of infections caused by an infected individual is defined as $\Re_{0}$. This section aims to present the calculation of the basic reproduction number ($\Re_{0}$).

The relevance of the basic reproduction number provides key information on the future behavior of an epidemic, which can contribute to decision making for epidemic control and eradication strategies. If $0 < \Re_{0} < 1$, the infection dies out in the long run, and if $\Re_{0} > 1$, the infection is able to spread within a population ^[21,22]. To find the basic reproduction number, we use the new generation matrix method presented in ^[21,22].

Due to the used methodology, we need an infection-free equilibrium point. Here we are interested in which conditions centered on $\Re_{0}$ and on the structure of the models (3.2)–(3.7) imply that this equilibrium point is globally asymptotically stable.

The infection-free equilibrium point of models (3.2)–(3.7) is

To compute the basic reproduction number, we use the next-generation matrix method presented in ^[21,22]. Using the notation in ^[21,22] on systems (3.2)–(3.7), matrices for the new infection terms, $F$, and the other terms, $V$, are given by:

and

where

Then, the basic reproduction number is

where

and

The stability of $\epsilon_{0}$ can be determined using the following theorem:

Theorem 3.2. ([23,Theorem 2]) Let

where $p, q \in \mathbb{Z}_{+}$ and

Define

then the disease-free equilibrium point $\epsilon_{0}$ of the models (3.2)–(3.7) is asymptotically stable if

for all roots $\lambda$ of the following equation

where $M_{1}$ is the matrix of the linearization of the models (3.2)–(3.7) at $\epsilon_{0}$.

This theorem can be demonstrated with analogous ideas in ^[19,24,25]. From Theorem 2 of ^[19,22], we have the following lemma characterizing the instability with the $\Re_{0}$:

Lemma 3.3. The disease-free equilibrium point $\epsilon_{0}$ of the models (3.2)–(3.7) is unstable if $\Re_{0} > 1$.

Now, we prove the global stability of the infection-free equilibrium point using the methodology presented in ^[26]. We can rewrite the models (3.2)–(3.7) as

where $S\in \mathbb{R}_{+}$ is the vector with susceptible individuals and $I \in \mathbb{R}^{5}_{+}$ have the other compartment of models (3.2)–(3.7).

The disease-free equilibrium is now denoted by

where

The conditions $(H^{\alpha}_{1})$ and $(H^{\alpha}_{2})$ below must be satisfied to guarantee the global asymptotic stability of $E_{0}^{\alpha},$

$(H^{\alpha}_{1})$: For ${}^{c}\mathbb{D}_{t}{^\alpha}S = F(S, 0),$ $S_{0}$ is globally asymptotically stable;

$(H^{\alpha}_{2})$: $G(S, I) = AI-G^{*}(S, I),$ $G^{*}(S, I)\geq 0,$ for $(S, I) \in \Omega_{\alpha},$

where $A$ is the Jacobian of G at $(S_{0}, 0_{\mathbb{R}^{5}})$ and is an $M$-matrix (the off-diagonal elements of $A$ are non-negative) and $\Omega_{\alpha}$ is the biologically feasible region. The following theorem allows us to characterize the global stability of the infection-free equilibrium point based on the behavior of $\Re_{0}$.

Theorem 3.3. The fixed point $E_{0}^{\alpha}$ is a globally asymptotically stable equilibrium of models (3.2)–(3.7) provided that $\Re_{0} < 1$ and that the conditions $(H^{\alpha}_{1})$ and $(H^{\alpha}_{2})$ are satisfied.

Proof. Let

We have the global stability of $S_{0}$ because $F(S, 0)$ is a linear equation. Let's

and

Since $S$ is always less than or equal to $N,$

Thus,

for all $(S, I) \in \Omega_{\alpha}.$ The $E_{0}^{\alpha}$ is globally asymptotically stable. □

3.3.   Incidence and mortality

Before including incidence and mortality in our model, we briefly present these epidemiologically important concepts. Incidence proportion (or attack rate or risk) is the number of new cases of disease during a specific time interval divided by the population at the start of the time interval. Mortality is related to the number of deaths caused by the epidemic ^[27,28].

In our study dedicated to TB, active cases are defined as active TB cases, including exposed cases entering the drug-sensitive and resistant compartments and latent cases entering an active state of the bacterium.

The following differential equations define the number of infected cases and the number of cases that die from the disease:

Equations (3.19) and (3.20) are incorporated into the systems (3.2)–(3.7), they define the active cases and deaths at time $t$. Then, incidence and mortality are defined as:

with $t$ as the current moment of time and $t-1$ as the previous moment of time.

In the sequel, we study how the implementation of 3HP and 1HP therapies in latents influences the incidence and mortality of TB.

4.   Numerical simulations

In this section, we perform computational illustrations of our model. We study the behavior of the $\Re_{0}$ and active TB compartments for different fractional orders

The objective is to study the basic reproduction number when implementing new treatments for latent TB and the impact of these implementations in the compartments compared to the maintenance of classical treatment schemes and to quantify this impact. For the computational simulations, we use the following values for the initial conditions and parameters:

The difference between the 3HP and 1HP is in the change of the value of the parameters $p_{T}$ and $\sigma_{T}$. For 3HP

and 1HP

We start with

to simulate the implementations of the onset of its activation in the population.

The numerical results of the Caputo derivative are obtained by the predictor-corrector PECE method of the Adams-Bashforth-Moulton type and implemented in MATLAB2021B. The methodology and theoretical results of the method can be found in ^[29,30], its application to HIV models in ^[19,25], to a TB model in ^[6], and to a model for the study of obesity in ^[31]. We use 15 years as the time horizon.

The value of the basic reproduction number for the scenario under study without the implementation of the 3HP and 1HP is approximately 1.5101, which is greater than unity, and it means that TB will spread in the population. Now, we study the behavior of $\Re_{0}$ when we implement the 3HP and 1HP concerning the selection rate for the application of these implementations and the success of this treatment. Then, the parameters $r_{L} \in [0.10, 0.50]$ and $p_{L} \in [0.80, 0.95]$. We can observe that when introducing these new therapies in latents, the basic reproduction number is less than unity. The higher the $\alpha$-value under study the higher the basic reproduction number but for any $\alpha$-value under study it is less than unity. This shows that the implementation of effective therapies in latents will reduce the impact of TB transmission since latents can develop active TB and be elements of transmission; see Figures 2–5.

We studied the impact of drug-sensitive and drug-resistant therapies on the basic reproduction number. For the parameters associated with the treatment event in drug-sensitive $\theta \in (0.65, 0.95)$ and $\psi \in (0.30, 0.98)$, we observed that the $\Re_{0}$ increases when we increase the $\alpha$-value, and in particular for $\alpha = 0.9$ and $\alpha = 1$ we found $\mathfrak{R}_0$ greater and less than unity. This is evidence that only the treatment of drug-sensitive individuals will not reduce the impact of TB in the population and that it is necessary to apply strategies in other compartments; see Figures 6–9.

For the parameters associated with the success of the treatment in the resistant $\alpha = 1.0$ and $p_{R} \in (0.4, 0.95)$, we find that as $\alpha_{R} \in (0.55, 0.92)$ increases, $\alpha$ increases, but not significantly. For the variations of these parameters and the different $\Re_{0}$-values, $\alpha$ is always less than unity. Thus, the variation of the parameters associated with the treatment of the resistors has a positive effect; see Figures 10–13.

In cases with drug-sensitive TB with only the classical treatment, the number of reported cases increases significantly over time, and the higher the $\alpha = 1.0$-values, the higher the number reported; see Figures 14–17.

The implementation of the 3HP and 1HP significantly reduces the number of reported cases of drug-sensitive TB compared to the application of the classical therapy alone. The higher the $\alpha$-values, the greater the reduction. Table 3 shows the difference between classical therapy and the implementation of the 3HP and 1HP. The 1HP reported better results than the 3HP in the reduction of drug sensitivity, and in Table 3, we can observe the difference between the numbers reported between 3HP and 1HP.

In resistant cases for $\alpha = 1.0$ with only the classical treatment, there is an increase at the beginning, followed by a decrease, and then an increase again until the end of the study. This shows that with the implementations, we avoid this volatile behavior; see Figure 18. For the other $\alpha = 0.5$-values under study, we have an analogous behavior; see Figures 19–21. At the end of the study, the number of cases reported concerning the initial one when we apply only the classical therapy is significantly exceeded, but with the implementations, the opposite situation occurs. Table 3 shows that 1HP reaches better results in the sense of reducing the number of cases.

When the $\alpha = 1.0$-values increase, the number of new active cases of TB increases, and the 3HP and 1HP implementations significantly reduce the number of active cases that transmit the bacteria; see Figures 22–25. The study without the implementations in this scenario shows the need to apply control strategies due to the growth of the number of active cases and the value of the basic reproduction number over 7.5 time.

The implementation of the 1HP therapy showed better results for the reduction of active cases in comparison with 3HP. Using Eqs (3.21) and (3.22). Table 4 shows the incidence of TB and we can observe that the higher the $\alpha$-value, the higher the incidence of TB and how the 3HP and 1HP implementations reduce the impact of TB on the population.

For deaths associated with TB, we see that all the different $\alpha = 1.0$-values grow over time, and in particular when we do not apply implementations on latents; see Figures 26–29. The implementations manage to reduce the number of TB-associated deaths over time, where the best results are reported by the 1HP. Table 4 shows TB mortality using Eq (3.22), and we see the impact it has when we do not implement new therapies and how new therapies reduce mortality.

Both incidence and mortality are reduced with the implementation of new therapies. Table 5 shows the difference in incidence and mortality between the scenario with only the classical therapy and with the 3HP and 1HP. This information is the number of cases avoided and deaths avoided with these therapies, and we can see that both implementations reduce the impact of TB concerning the number of cases avoided, and deaths avoided and the 1HP showed better results.

It is important to mention that by applying the 3HP and 1HP implementations, the impact on the population is significantly reduced, but with them only in the scenario under study, we did not manage to eliminate TB in 15 years, so it is necessary to make use of other implementations and control strategies.

5.   Conclusions

In this paper, we present a mathematical model for the study of TB treatment, taking into account the treatment of latent TB and the new implementations of the 3HP and 1HP. The objective is to test the impact of 3HP and 1HP on the dynamics and to compare these implementations for the reduction of active TB cases. We use fractional order derivatives in the Caputo sense to take advantage of the benefits of this technique, mainly in the memory effect.

We performed computational simulations to validate our model for different fractional orders. We can conclude that in the studied scenario, the 3HP and 1HP implementations reduce the number of cases of active TB (drug-sensitive and drug-resistant) and TB mortality significantly. However, the 1HP showed better results. In the study of the basic reproduction number concerning the selection rate for the new implementations, 3HP and 1HP, and their success, we found that all values of the fractional derivative are less than unity, so these implementations reduce the impact of TB in the population.

We conclude that although these implementations reduce the impact of TB to eradicate TB in society, other strategies are needed.

This work can contribute to decision-making regarding the implementation of these therapies in a population by the validation of the model in the conditions of the studied setting.

Conflict of interest

The authors declare that they have no conflicts of interest in this paper.