Conditions for extinction and ergodicity of a stochastic Mycobacterium tuberculosis model with Markov switching

1.   Introduction

Tuberculosis (TB), caused by mycobacterium tuberculosis (Mtb), remains one of the leading causes of death worldwide, surpassing even acquired immune deficiency syndrome (AIDS) ^[1,2]. Approximately $25\%$ of the global population carries Mtb, with most progressing to latent infection. This latent state can persist for life or re-emerge as active disease, underscoring the need to understand Mtb-host dynamics. As a result, many studies were dedicated to exploring these interactions ^[3,4]. For instance, Ibarguen-Mondragon et al. ^[5] proposed a mathematical model describing the growth of Mtb populations:

Here, $\bar{M}_U(t), \; \bar{M}_I(t), \; \bar{B}(t)$ and $\; \bar{T}(t)$ represent the population densities of normal macrophages, infected macrophages, extracellular Mtb and T cells, respectively. Some main parameters of system (1.1) are summarized in Table 1.

To simplify the model, they introduce the following change of variables:

Replacing the new variables, the system (1.1) becomes

where

In ecosystems, many of the main parameters are affected by environmental white noises such as drought-fire interactions and species invasions, and therefore should generally display stochastic disturbances ^[6,7,8,9,10]. Stochastic models have been widely employed to capture the dynamics of various infectious diseases, including measles, malaria, tuberculosis, smallpox, and AIDS. However, few stochastic models have explored the impact of Mtb growth on infection outcomes.

However, the majority of ecosystems will eventually change due to many natural elements like temperature, precipitation, and PH. Furthermore, we see that during the warm season, the recruitment and death rates of both healthy and infected macrophages will change significantly from those during the cold season. Similarly, changes in nutrition or food resources commonly impact the intrinsic reproduction rate. Colored noise (or telegraph noise) is often used to describe the transition between different environmental states, such as from the rainy season to the dry season. The switching is memoryless and the waiting time for the next switching is exponentially distributed. Therefore, a continuous-time Markov chain $\varpi(t), \; t\geq 0$ with finite-state space $\mathbb {S} = \{1, 2, \cdots, N\}$ is used to represent random switches between environmental states ^{[11,12,13,14,15]}. Taking into account the sensitivity to environmental states, let us investigate time-varying parameters with various discrete values affected by colored noises. We will consider time-varying parameters influenced by both white and colored noise and introduce this noise into system (1.2) as follows:

where $B_1(t), \, B_2(t), \, B_3(t)\, \text{and}\, B_4(t)$ are mutually independent standard Brownian motions and the Markov chain $\varpi(t), \; t\geq 0$ with values in a finite state space $\mathbb {S} = \{1, 2, \cdots, N\}$. We assume that Brownian motion and Markov chain are independent. The generator matrix $\varGamma = (\gamma_{ij})_{N\times N}$ of the Markov chain is given by

where $\triangle t > 0, \gamma_{ij} > 0$ denotes the transition rate from $i$ to $j$ if $i\neq j$ while $\sum\limits_{j = 1}^{N}\gamma_{ij} = 0$. In addition $(\varpi(t))_{t\geq0}$ is irreducible and has a unique stationary distribution $\pi = (\pi_1, \pi_2, \cdots, \pi_N)$ satisfying the conditions $\pi\varGamma = 0, \; \sum\limits_{k = 1}^{N}\pi_{k} = 0.$

This paper aims at establishing some criteria for the existence of ergodic stationary distribution and extinction of mycobacterium tuberculosis model, which is almost a void in this area. As far as we know, this type of model has received little attention. There is not much research on the stochastic epidemic model in the literature because of how difficult it is to handle discrete Markov switching and remove linear perturbation terms. Unlike deterministic models, it is difficult to analyze the disease persistence and extinction of system (1.3) because of the stochastic fluctuations of each compartment in disease transmission; the stable equilibrium of system (1.3) will no longer exist. In this way, analyzing the persistence and extinction of tuberculosis disease is a challenging task. We will provide the relevant threshold dynamics and ergodic properties of the system (1.3) to the best of our ability.

The structure of the paper is organized as follows: Section 2 introduces necessary notations and auxiliary lemmas. Section 3 investigates the conditions for the existence and uniqueness of a global positive solution to system (1.3). Section 4 applies stochastic Lyapunov methods to establish the ergodicity and positive recurrence of the stochastic Mtb model under regime switching. Finally, we derive the sufficient condition for extinction.

2.   Preliminaries

In this section, we will introduce three important lemmas for the subsequent dynamical analyses.

Lemma 2.1. (Has'minskii ^[16]) Assume that for any $i\neq j \in \mathbb {S},$ such that $\gamma_{ij} > 0.$ If the following conditions are satisfied:

(I) For any $k\in \mathbb {S}$ and for all $Y\in \mathbb {R}^n,$ $C(Y, k)$ is symmetric and

with some constant $\rho\in (0, 1].$

(II) There exists a nonempty bounded open set $U\in \mathbb {R}^n$ with compact closure, satisfying that, for each $k\in \mathbb {S},$ there exists a nonnegative function $V(\cdotp, k): U^c\times \mathbb {S}\rightarrow \mathbb {R}$ such that $V(\cdotp, k)$ is twice continuously differentiable and for some $\varrho > 0,$

then the solution $(Y(t), \varpi(t))$ of system (2.1) is positive recurrent and ergodic. It shows that $(Y(t), \varpi(t))$ has a unique stationary density $\mu (\cdotp, \cdotp),$ and for any Borel measurable function $\varphi (\cdotp, \cdotp): \mathbb {R}^n\times\mathbb {S}\rightarrow \mathbb {R}^n$ such that $\sum\limits_{k \in \mathbb {S}}\int_{\mathbb {R}^n}|\varphi(y, k)|\mu(\mathrm{d}y, k) < \infty,$ we have

Then, the ergodicity of Markov chain $\varpi (\cdotp)$ implies that $\lim\limits_{t \rightarrow \infty}\frac{1}{t}\int_0^t\varphi(\varpi(s))\mathrm{d}s = \sum\limits_{k \in \mathbb {S}}\pi_k\varphi(k)\; \; \; a.s.$

Lemma 2.2. The following linear system

where

then (2.2) has a unique solution $(g_1(1), \cdots, g_1(N), g_2(1), \cdots, g_2(N), g_3(1), \cdots, g_3(N))^T\gg 0.$

Proof. System (2.2) can be rewritten in the following form,

where $G\in \mathbb {R}^{3N},$ $H = (\beta(1)+c_2\upsilon(1), \cdots, \beta(N)+c_2\upsilon(N), \gamma_U(1), \cdots, \gamma_U(N), \alpha_T(1), \cdots, \alpha_T(N))^T$ and

Clearly, $A\in \mathbb {Z}^{3N\times 3N}.$ For each $k = 1, \ldots, N,$ consider the leading principal submatrix

and

Obviously, $A_k, \, A_{N+k, }\, A_{2N+k} \in \mathbb {Z}^{k\times k}.$ Additionally, each row of submatrix $A_k$ has the sum

For submatrix $A_{N+k},$

And for submatrix $A_{2N+k},$

By applying Lemma 5.3 of ^[17], we get $\text{det} A_k > 0, k = 1\ldots, 3N.$ In other words, we've shown that all the leading principal minors of $A$ are positive. Using Theorem 2.10 in ^[16] indicates that $A$ is a nonsingular M-matrix and for the vector $H \geq 0 \in \mathbb {R}^{3N},$ the linear Eq (2.2) has a unique solution $G = (g_1(1), \cdots, g_1(N), g_2(1), \cdots, g_2(N), g_3(1), \cdots, g_3(N))^T\gg 0.$ On the other hand, by system (2.2), we can easily observe that $g_1(k), g_2(k)$ and $g_3(k)$ should be positive, $k = 1, \ldots, N.$

Lemma 2.3. (^[18]) Let $Z(t)$ be the solution of the auxiliary stochastic differential equation

with the initial value $Z(0) = M_U(0) > 0,$ Then $M_U(t)\leq Z(t)$ for any $t\geq 0,$ a.s. Moreover, $(Z(t), \varpi(t))$ is positive recurrent and has the following property:

where

is determined by the following linear equation

3.   Existence and uniqueness of the global positive solution

When studying the dynamical behavior of an epidemic model, it is important to consider whether the solution is global and positive.

Theorem 3.1. For any initial value $(M_U(0), M_I(0), B(0), T(0), \zeta(0))\in \mathbb{R }^4_+\times \mathbb{S },$ there exists a unique solution $(M_U(t), M_I(t), B(t), T(t), \varpi(t))$ to system (1.3) on $t\geq 0$ and the solution will remain in $\mathbb{R }^4_+\times \mathbb{S }$ with probability one (a.s.).

Proof. Since the coefficients of $(1.3)$ satisfy the locally Lipschitz continuous condition, thus the system (1.3) has a unique local solution $(M_U(t), M_I(t), B(t), T(t), \varpi(t))\in \mathbb{R }^4_+\times \mathbb{S }$ on $t \in [0, \tau_e],$ where $\tau_e$ is an exposure time. Next, we claim that the solution is global, i.e $\tau_e = +\infty.$ Similar to the proof of Zu et al. ^[19] and Liu et al. ^[20], we will only show the key step is to construct a nonnegative Lyapunov function $Q_0: \mathbb{R }^4_+\rightarrow \mathbb{R }_+$ satisfying

where $\Theta$ is a positive constant. Define

where $a, b \, \text{and } \, d$ are positive constants to be defined later. Based on the basic inequality $u-1-\ln u \geq 0,$ for any $u > 0$, we have

Making use of Itô's formula to $Q_0$, we obtain

where $LQ_0: \mathbb{R }^4_+\rightarrow \mathbb{R }$ is defined by

Choose $a = \max\{\frac{\check{r}+2\check{k}_I}{\hat{\mu_I}}, \; \frac{\check{\gamma}_U}{\hat{\mu}_U}\}, b = \frac{\hat{\mu}_B}{\check{\beta}}, \; d = \frac{\hat{\mu}_T}{\check{\alpha}_T},$ such that $-a\hat{\mu}_I+\check{r}+2\check{k}_I\leq0, \; \; \; -a\hat{\mu}_U+\check{\gamma}_U\leq0,$ then,

where $\Theta$ is a positive constant. The following proof is almost the same as those in Theorem 2.1 of Li et al. ^[21]. Hence, we omit it here.

4.   Ergodic stationary distribution

In this part, we demonstrate the existence of a unique ergodic stationary distribution, which suggests that the virus is widespread, based on the theory presented in Lemma 2.1.

Define the critical condition

where $g(k) = (g_1(k), g_2(k), g_3(k)) ^T$ is the solution of the linear system (2.2) and $c_2$ is defined in Lemma 2.2.

Theorem 4.1. If $R_0^* < 0$ and $\check{\sigma}_1^2 < 2\hat{\mu}_U$, $\check{\sigma}_2^2 < \hat{\mu}_I$, $\check{\sigma}_4^2 < 2\hat{\mu}_T$ are satisfied, then for any initial value $(M_U(0), M_I(0), B(0), T(0), \zeta(0))\in \mathbb{R }^4_+\times \mathbb{S },$ the solution $(M_U(t), M_I(t), B(t), T(t))$ of system $(1.3)$ is positive recurrent and has a unique ergodic stationary distribution $\phi (\cdotp, \cdotp).$

Proof. Since the diffusion matrix

is positive definite, which implies that condition (I) in Lemma 2.1 is satisfied. Next we will prove the condition (II) holds. Define a $C^2-$function $\widetilde{Q}: \mathbb{R }^4_+\times \mathbb{S }\longrightarrow \mathbb{R }$ as follows:

It is clear that there is a unique point $(\bar{M}_U(k), \; \bar{M}_I(k), \; \bar{B}(k), \; \bar{T}(k), k),$ which is the minimum value of $\widetilde{Q}(M_U, M_I, B, T, k).$ Define a nonnegative $C^2-$ function $Q: \mathbb{R }^4_+\times \mathbb{S }\longrightarrow \mathbb{R }$ in the following from

where $(M_U, M_I, B, T, k)\in (\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times \mathbb{S },$ and $n > 1$ is a sufficiently large integer, $Q_1 = -\ln M_U-c_1\ln M_I-c_2\ln B,$ $Q_2 = g_1(k)B+c_2g_2(k)M_U+c_1g_3(k)T,$ $Q_3 = -\ln M_U-\ln B-\ln T,$ $Q_4 = \frac{1}{2}\bigg(M_U+M_I+\frac{\hat{\mu}_I}{4\check{r}}B+\frac{\hat{\mu}_I}{4\check{k}_I}T\bigg)^2$ and

$g(k): = (g_1(1), \cdots, g_1(N), g_2(1), \cdots, g_2(N), g_3(1), \cdots, g_3(N))^T,$ is the unique solution of system (2.2), $\omega (k): = (\omega(1), \cdots, \omega(N))^T$ will be determined later and $M_0 > 0$ is a sufficiently large number satisfying the following condition,

where

Applying the Itô's formula to $Q_1$ and $Q_2$, we have

and

where $g_1(k), g_2(k), g_3(k)$ are defined in Lemma 2.2. In view of (4.5), (4.6) and $dx-ex^2\leq\frac{d^2}{4e}\; (e > 0), \; \; \; \forall x\in \mathbb{R },$ we obtain

where $R_0(k) = -3\sqrt [3]{c_1c_2\mu_U(k)\beta(k)r(k)}+c_1\bigg(\mu_I(k)+\frac{1}{2}\sigma_2^2(k)\bigg)+c_2\bigg(\mu_B(k)+\frac{1}{2}\sigma_3^2(k)\bigg)+\mu_U(k)+\frac{1}{2}\sigma_1^2(k)+\frac{1}{4}\upsilon(k) g_1(k)+c_2g_2(k)\mu_U(k).$

Since the generator matrix $\Gamma$ is irreducible, for $R_0 = (R_0(1), \ldots, R_0(N))^T,$ there exists $\omega = (\omega(1), \ldots, \omega(N))^T$ satisfying the following Poisson system

which implies that

Substituting this equality into $(4.7)$

Employing the Itô's formula to $Q_3$ and $Q_4$, one has

and

where

It follows from (4.8)–(4.10) that

Next, we construct a compact subset $U$ such that the condition (II) in Lemma 2.1 holds. Define the following bounded set

where $0 < \varepsilon < 1$ is a sufficiently small constant. In the set $\mathbb{R }^4_+\setminus U = U^C$ we choose $\varepsilon$ satisfying the following condition

where

For convenience, we can divide $\mathbb{R }^4_+ \setminus U = U^C$ into eight domains,

Obviously, $U^C = U_1\cup U_2 \cup U_3\cup U_4\cup U_5\cup U_6\cup U_7\cup U_8.$ Next, we will prove that

Case 1. For any $(M_U, M_I, B, T, k) \in U_{1}\times \mathbb{S },$ according to (4.4) and (4.11), we obtain that

Case 2. For any $(M_U, M_I, B, T, k) \in U_{2}\times \mathbb{S },$ in view of (4.12), we can obtain

Case 3. For any $(M_U, M_I, B, T, k) \in U_{3}\times \mathbb{S },$ according to (4.12), we can deduce that

Case 4. For any $(M_U, M_I, B, T, k) \in U_{4}\times \mathbb{S },$ from condition (4.12), we can deduce that

Case 5. For any $(M_U, M_I, B, T, k) \in U_{5}\times \mathbb{S },$ in view of (4.13), we can deduce that

Case 6. For any $(M_U, M_I, B, T, k) \in U_{6}\times \mathbb{S },$ by condition (4.14), we can conclude that

Case 7. For any $(M_U, M_I, B, T, k) \in U_{7}\times \mathbb{S },$ in view of (4.15), we obtain that

Case 8. For any $(M_U, M_I, B, T, k) \in U_{8}\times \mathbb{S },$ from (4.16), it is deduced that

Then the assertion (4.17) is verified, i.e., condition (II) of Lemma 2.1 holds.

5.   Exponential extinction

Define

where $l_1(k) = \frac{\beta(k)+c_2\upsilon(k)}{\beta(k)}, \, \, \, l_2(k) = \frac{\beta(k)+c_2\upsilon(k)}{2\gamma_U(k)},$ $c_2 = \frac{\bigg(\sum\limits_{k = 1}^N\pi_k(\mu_U(k)\beta(k)r(k))^{\frac{1}{3}}\bigg)^3} {\bigg(\sum\limits_{k = 1}^N\pi_k(\mu_I(k)+\frac{1}{2}\sigma_2^2(k))\bigg)\bigg(\sum\limits_{k = 1}^N\pi_k(\mu_B(k)+\frac{1}{2}\sigma_3^2(k))\bigg)}$ and $g_1(k)$ is the solution of the linear system (2.2).

Theorem 5.1. Assume that $R_0^E < 0$ for any initial value $(M_U(0), M_I(0), B(0), T(0), \zeta(0)) \in \mathbb{R }^4_+\times \mathbb{S },$ the solution $(M_U(t), M_I(t), B(t), T(t), \varpi(t))$ of system (1.3) will follow

which means that the disease of system (1.3) will exponentially go to extinction with probability one, where

Proof. Define a $C^2$-function $H: \mathbb{R }^4_+\times \mathbb{S }\rightarrow \mathbb{R }$ as follow,

Employing Itô's formula to $H$ and applying (5.1), one has

We use the following relationship,

then,

Integrating from $0$ to $t$ and dividing t on both sides of (5.2), we have

By Lemma 2.3, we have

Using the strong law of large number for local martingale, one has

Taking the superior limit on both sides of (5.3) and applying the ergodicity of Markov chain $\varpi(t)$, we obtain

Hence, we can equivalently obtain

This completes the proof.

6.   Numerical simulations

In this section, some numerical examples are provided to support our theoretical findings. Using Milstein's high-order method, the corresponding discretization equation of system (1.3) is

here $\eta_{1, j}, \, \eta_{2, j}, \, \eta_{3, j}, \, \eta_{4, j}$ are $N(0, 1)$ distributed independent Gaussian random variables.

Let $N = 2$ and the generator $\Gamma = (\gamma_{ij})_{2\times 2}$ of the Markov chain be

By solving $\pi \Gamma = 0,$ the stationary distribution $\Gamma$ follows $\pi = (\pi_1, \pi_2) = (\frac{8}{13}, \frac{7}{9}).$

Example 6.1. Take initial value $(M_U(0), M_I(0), B(0), T(0)) = (1, 5, 3.5, 0.65)$ and

Case 1. Choose $\big(\sigma_1(1), \sigma_1(2)\big) = (0.03, 0.01), \; \big(\sigma_2(1), \sigma_2(2)\big) = (0.002, 0.001),$ $\big(\sigma_3(1), \sigma_3(2)\big) = \big(\sigma_4(1), \sigma_4(2)\big) = (0.006, 0.005),$ then $R_0^* = -1.873 < 0.$ By Theorem 4.1, we obtain that there exists a unique ergodic stationary distribution of system (1.3). Our simulations confirm these results: The sample paths of $M_U(t), M_I(t), B(t), T(t),$ and their corresponding probability density function (PDF) are shown in Figure 1. Figure 2 shows the corresponding movement of Markov chain $\big(\varpi(t)\big)_{t\geq 0}$ in the state space $\mathbb {S} = \{1, 2\}.$

Case 2. Choose $\big(\sigma_1(1), \sigma_1(2)\big) = (0.3, 0.1), \; \big(\sigma_2(1), \sigma_2(2)\big) = (0.02, 0.01),$ $\big(\sigma_3(1), \sigma_3(2)\big) = \big(\sigma_4(1), \sigma_4(2)\big) = (0.006, 0.005).$ Simple computation $R_0^* = -1.7581 < 0.$ Then from Theorem 4.1 it follows that system (1.3) has a unique stationary distribution. Simulations are presented in Figure 3. By comparing with Figure 1, the numbers of $M_U(t), M_I(t), B(t), \text{and}\; T(t)$ are largely fluctuated by the stochastic noises.

Case 3. Choose $\big(\sigma_1(1), \sigma_1(2)\big) = (0.6, 0.2), \; \big(\sigma_2(1), \sigma_2(2)\big) = (0.5, 0.4),$ $\big(\sigma_3(1), \sigma_3(2)\big) = \big(\sigma_4(1), \sigma_4(2)\big) = (0.9, 0.8).$ We can easily obtain $R_0^* = 0.3042 > 0,$ we can not determine whether there exists an ergodic stationary distribution. From Figure 4, we can see that the disease of system (1.3) will be extinct in a long time. From Figures 1, 3 and 4, we can find that when white noise intensity $\sigma^2(k)$ increases, infected populations tend to go extinct faster.

On the left column of Figure 5, the red, blue and green lines represent the sample paths of $M_U(t), M_I(t),$ $B(t), \text{and }\; T(t)$, when there is only one state $k = 1, k = 2$ and switching between states $k = 1, 2$. Similarity, On the right column of Figure 5, the red, blue and green lines represent the PDF of $M_U(t), M_I(t), B(t), \text{and }\; T(t)$. It is displayed directly that the green line is located between the red and the blue lines. That is to say the switching state is located between states $k = 1$ and $k = 2$.

Example 6.2. We choose $\big(\alpha_T(1), \alpha_T(2)\big) = (0.015, 0.0097), \; \big(\mu_I(1), \mu_I(2)\big) = (0.7, 0.8), \; \big(r(1), r(2)\big) = (0.05334, 0.04), \; \big(k_I(1), k_I(2)\big) = (0.3636, 0.32),$ and $\big(\sigma_i(1), \sigma_i(2)\big) = (0.01, 0.02), \, i = 1, 2, 3, 4.$ Other coefficients are the same as in Example 6.1. By direct calculation, we derive $R_0^E = -0.1032 < 0.$ Then the disease of system (1.2) will be extinct in a long time, which can be verified in Figure 6.

7.   Conclusions

This paper is devoted to studying a stochastic mycobacterium tuberculosis model, that is perturbed by white and colored noises. First, we show that the unique solution of system (1.3) is global and positive with probability one. In order to establish the existence of an ergodic stationary distribution, we construct a stochastic Lyapunov function with regime switching. Different switching parameters correspond to different peaks in the distribution function, and each peak represents the equilibrium value. Further, we can infer from Example 6.1 that large perturbations can change population dynamics, whereas smaller perturbations can lead to disease persistence.

Some interesting topics deserve consideration. Such as considering mean-reverting Ornstein-Uhlenbeck processes, non-Gaussian Levy noise, and impulsive perturbations on system (1.2). We can also use the method of this paper to study other epidemic models. We leave these cases for our work.

Author contributions

Ying He: Conceptualization, Investigation, Formal analysis, Writing – review and editing. Bo Bi: Formal analysis, Writing – review and editing, Numerical simulation. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work is supported by Hainan Provincial Natural Science Foundation of China (No. 122RC679, 121RC554), Talent Program of Hainan Medical University (No. XRC2020030) and Northeast Petroleum University Special Research Team Project (No. 2022TSTD-05).

Conflict of interest

The authors declare there is no conflict of interest.