Ergodic stationary distribution and extinction of stochastic pertussis model with immune and Markov switching

1.   Introduction

Pertussis, commonly recognized as whooping cough, is an extremely contagious respiratory infection caused by Bordetella pertussis bacteria ^[1]. The infectious disease is primarily spread through direct contact, such as respiratory droplets. It primarily affects the respiratory system and is characterized by severe coughing. Pertussis predominantly affects infants and young children but can also impact adolescents and adults, often leading to severe complications, particularly in vulnerable populations. The early signs of pertussis resemble common cold symptoms and are easily misdiagnosed. Therefore, it is difficult to quickly identify pertussis in the early stages. The treatment of pertussis mainly relies on antibiotics and symptomatic treatment. Currently, there is no specific drug for pertussis, resulting in no effective elimination of the pathogen. According to data from the Chinese Center for Disease Control and Prevention, there were nearly 500,000 cases of pertussis infections in China in 2024.

It is worth noting that over the past six decades, injectable vaccines have become an effective strategy to prevent diseases. Before the pertussis vaccine was widely available, pertussis caused a significant number of childhood deaths around the world every year. As the pertussis vaccine has been promoted globally, the mortality rate has decreased dramatically. However, in recent decades, there has been a resurgence in pertussis cases ^{[2,3,4,5,6,7,8]}. In 2014, $24.1$ million cases of pertussis were reported worldwide, resulting in 160,700 deaths among children under $5$ years old, particularly in developing countries ^[9]. This trend can be attributed to two main factors. On the one hand, vaccination rates in developing countries are generally low, resulting in insufficient vaccine coverage to form effective herd immunity. Herd immunity is the idea that when enough of the population is vaccinated, the entire community becomes more resistant to the disease, thereby protecting individuals who cannot be vaccinated, such as newborns or those with immune deficiencies. However, if vaccine coverage is insufficient, the virus can still spread among people, leading to outbreaks ^[10]. On the other hand, the immunity provided by the vaccine is not permanent, and the effectiveness of the vaccine diminishes over time, especially for adult and adolescent populations, where the protective effect is relatively weak. This means that even if children are fully vaccinated, they may lose their immunity to pertussis over time, increasing the risk of reinfection. It is also important to note that individuals who recover from an infection may gain immunity, but that does not offer complete protection ^[11]. The impact of immunity on disease transmission has gained great attention, while there are relatively few works related to immune waning. Based on the above background, this paper focuses on a pertussis model with long-term immune waning and natural immune boosting.

Mathematical models are essential for characterizing transmission dynamics of infectious diseases. Over the past few decades, mathematical models of pertussis have provided important insights into the occurrence of new pertussis cases. Tian and Wang ^[12] proposed a pertussis model incorporating recessive infection, thereby laying a solid foundation for understanding the impact of latent infection on pertussis recurrence. Given the multifaceted nature of pertussis infection, a model accounting for multiple infection pathways was given by ^[13]. Research on the effects of vaccination on pertussis dynamics has also garnered significant attention. Rozhnova and Nunes ^[2] focused on the long-term trends of pertussis prior to vaccination, aiming to minimize the number of free parameters. Safan et al.^[14] established a pertussis model based on data from newborns, which concluded that symptomatic infections can be completely eradicated through vaccination. Conversely, Águas proposed that the resurgence of pertussis may result from reduced transmission, independent of vaccination ^[15]. Most existing studies ^{[1,12,13,14,15]} have concentrated on deterministic models, often neglecting the influence of environmental noise and its role in pertussis dynamics. Meanwhile, the dynamic behavior of most epidemics (such as pertussis) is greatly influenced by random factors ^[16,17]. Due to stochastic models having advantages in adapting to complex future conditions ^{[18,19,20,21,22,23,24,25]}, we apply nonlinear perturbations ^[21] to a deterministic pertussis model. During the modeling process, we mainly considered the influence of environmental factors such as temperature and humidity changes on the transmission of pertussis. The transmission of pertussis has seasonal variations ^[26,27], mainly reflected in terms of temperature and humidity. In order to more accurately capture the seasonal transmission characteristics of pertussis, i.e., peak of infection in summer and winter, while infections in spring and autumn are relatively rare, we introduced a Markov switching process ^[25] in the model, aiming to simulate the influence of seasonal fluctuations of environmental factors (such as temperature and humidity) on the transmission of pertussis.

As far as we are aware, up to now, a pertussis model taking both immune and random factors into consideration has not received any attention. In order to fill this gap, this paper proposes a stochastic pertussis model with Markov switching. We apply Lyapunov functions and the Has'minskii theorem to analyze the pertussis model. The primary contributions and innovations of this paper are described below:

● The transmission rate is more sensitive to the external environment. A stochastic pertussis model with Markov switching of transmission rate is discussed in this paper. This model is an extension of Lavine et al. ^[10].

● In response to the recurrence and state-switching phenomena of pertussis, it is essential to implement appropriate measures. To this end, sufficient conditions are presented to ensure the extinction of a stochastic pertussis model with Markov switching. Additionally, considering persistence of the pertussis is equally important. We derive persistence conditions for pertussis by proving ergodicity.

The remainder of this paper is structured as follows: Section 2 introduces the mathematical model and necessary lemmas. In Section 3, we establish sufficient conditions for disease extinction. The existence of a unique ergodic stationary distribution (ESD) is demonstrates in Section 4. Section 5 shows numerical simulations to validate our theoretical results and explore the effects of random perturbations and natural immune boosting on disease transmission. Section 6 offers a discussion of our conclusions.

2.   Model formulation

Lavine et al.^[10] proposed a deterministic model that is outlined as follows:

where $S(t)$, $I(t)$, $W(t)$, and $R(t)$ are the susceptible compartment, the infectious compartment, the waning compartment, and the recovered compartment at time $t$, respectively. It is worth noting that $R(t)$ consists of individuals who have recently recovered or been boosted, possessing a strong immune response. $W(t)$ includes individuals who maintain immunity against infection but can further strengthen their immune response upon re-exposure $\kappa \beta I(t)W(t)$. If individuals in $W(t)$ are not re-exposed, they will eventually lose their immunity entirely and revert to the susceptible group $\eta W(t)$. Besides that, $\mu$ represents the death rate, $\nu$ denotes the vaccination probabilities specific to different age groups, $\beta$ characterizes the transmission rate, $\eta$ represents the rate at which immunity wanes when there is no immune boosting, $\gamma$ represents the rate of recovery, and $\kappa$ represents the coefficient for boosting. Denoting $N(t) = S(t)+I(t)+W(t)+R(t)$ as the total polulation, and on the basis of (2.1), we can obtain $dN(t) = (\mu-\mu N(t))dt$. By integrating, we can get

it will come out $N(t) \leq 1$, while $N(0)\leq 1$.

Denote the possible region of the determinist system (2.1)

The system (2.1) possesses a disease-free equilibrium point $E_0 = (S_0, I_0, W_0, R_0) = (\frac{\mu (1-\nu)+\eta W^0}{\mu}, 0, \frac{\eta\mu \nu }{(\mu+\eta)^2}, \frac{\mu \nu}{\mu+\eta})$. Combining the theories of Ma et al.^[18] and Sahu et al.^[19], the basic reproduction number $\mathcal R_0$ of the system (2.1) is obtained as follows:

By applying the global stability theory of equilibrium, system (2.1) has the following threshold dynamics:

● $E_{0}$ is globally asymptotically stable ^[10], which implies the number of infectious will decrease to zero.

● If $\mathcal R_0 > 1$, the endemic equilibrium $E^* = (S^*, I^*, W^*, R^*)$ is globally asymptotically stable and satisfies the following equation:

The proof of asymptotic stablility of $E^*$ is similar to that in Theorem 3.2 of ^[20], so we omit it here.

Environmental noise is ubiquitous and can significantly impact the transmission of pertussis. The deterministic model (2.1) has inherent limitations in accurately predicting future dynamics. Linear perturbation, a simplifying assumption commonly employed, is commonly used to describe disease spread ^[29,30,31]. Based on ^{[21,23,25,32]}, we use nonlinear perturbation for the system (2.1). This study investigates the natural mortality rate of the system (2.1) under nonlinear perturbations, such as $\mu S(t) \rightarrow \mu +(\sigma_{ii}+\sigma_{ij}S(t))S(t)dB_{i}(t)$, which is given by

where $B_{i}(t)(i = 1, 2, \ldots, 4)$ refers to four standard Brownian motions, each independent. $\sigma_{i1}^{2} > 0$, $\sigma_{i2}^{2} > 0$ describe the intensities of nonlinear and linear perturbations on $B_{i}(t)$, respectively, $i = 1, 2, \ldots, 4.$ According to Mao ^[18], the processes $B_{i}(t)(i = 1, 2, \ldots, 4)$ are defined on a complete probability space $\{\Omega, \mathcal{F}$, $\{\mathcal{F}_t\}_{t > 0}$, $\mathbb{P}\}$, where the filtration $\{\mathcal{F}_{t}\}_{t\geq0}$ is right continuous and increasing, provided that $\mathcal{F}_{0}$ contains all $\mathbb{P}\mathrm{-null}$ sets.

In the actual pertussis transmission model, changes in the external environment, especially alterations in temperature and humidity, can lead to variations in parameters. In the modeling process, fixing parameters cannot represent such state changes. Similarly, introducing random noise $\sigma$ cannot accurately depict such changes. The mortality rate $\mu$, the immunity wanes $\eta$, and vaccination probabilities $\nu$ are less affected by environmental changes, while the contact rate $\beta$ is more sensitive to environmental changes ^[33]. In this paper, we introduce Markov switching to describe the state transitions caused by environmental changes and simplify the model by only considering the influence of the state change on $\beta$, described in detail as follows:

where $(r(t))_{t\geq0}$ is a right-continuous Markov chain taking values in a finite space $\mathcal N = \{1, 2, \ldots, N\}$. The following is transition probability:

where the time increment $\Delta > 0$, and $q_{ij} > 0$ is the transition coefficient from $i$ to $j$ if $i\neq j$. For any $i\in\mathcal{N}$, if $j = i$, $q_{ii} = -\sum_{i \neq j}q_{ij}$. $o(\Delta)$ signifies $\lim_{\Delta \rightarrow 0}o(\Delta)/ \Delta = 0$. In addition, $(r(t))_{t\geq0}$ is irreducible and has a unique stationary distribution $\pi = \{\pi_{1}, \pi_{2}, \ldots, \pi_{N}\}\in \mathbb R^{1\times n}$, which is determined by $\pi\Gamma = 0$ subject to $\sum_{k = 1}^N\pi_k = 1$, $\pi_k > 0$.

We need to give several mathematical notations that will be used in this paper. Let an N-dimensional vector $\overrightarrow{A} = (A(1), A(2), \ldots, A(N))^T$, assuming that $\check{A} = \max_{k\in\mathbb{S}}\{A(k)\}$, $\hat{A} = \operatorname*{min}_{k\in\mathbb{S}}\{A(k)\}$.

To investigate the dynamical behavior of the stochastic epidemic model (2.5), it is essential to establish the existence of a global positive solution. Since $S(t)$, $I(t)$, $W(t)$, and $R(t)$ represent numbers of different kinds of individuals, these values should remain non-negative. And based on the background of pertussis transmission, we assume $\mu\geq 0, \nu\geq 0, \beta \geq 0, \gamma \geq 0, \kappa\geq 0, \eta\geq 0$. In this section, we will show that system (2.5) possesses a unique global positive solution.

Theorem 2.1. For any initial condition $(S(0), I(0), W(0), R(0), r(0))\in \mathbb R_+^4\times \mathcal{N}$, system (2.5) admits a unique solution $(S(t), I(t), W(t), R(t), r(t))$ on $t\geq0$, and the solution will remain in $\mathbb R_+^4\times \mathcal{N}$ with probability one, a.s..

Proof. We omit the proof, which is similar to Zu et al. ^[24] and Zhou et al. ^[25]. □

Remark 2.1. The existence and uniqueness of the global positive solution are essential for ensuring the reliability of predictions in the stochastic pertussis model. Without establishing these fundamental properties, the model's predictive power would be compromised, as multiple or non-positive solutions could arise, leading to ambiguous or incorrect interpretations of the disease dynamics.

3.   Extinction of the system (2.5)

The extinction of diseases is an important topic discussed in the field of biomathematics. In this section, we seek to establish the sufficient conditions for the extinction of random systems (2.5). To this end, a value $\tilde{\mathcal {R}}_0$ associated with $\mathcal R_{0}$ is defined below

Theorem 3.1. For any initial value $(S(0), I(0), W(0), R(0))\in \mathbb R_{+}^{4}\times\mathcal N$, let $(S(t), I(t), W(t), R(t))$ be the solution of stochastic system (2.5). If $\tilde{\mathcal {R}}_0 < 1$ holds, then the disease of the stochastic system (2.5) will go extinct.

Proof. Utilizing the It$\hat o$'s formula yields

Integrate both sides of Eq (3.1) from $0$ to $t$, we can obtain

where $\Phi(t) = \int_{0}^{t}(\sigma_{21}(r(u))I+\sigma_{22}(r(u)))dB_{2}(u)$ is a continuous real-valued local martingale and its quadratic variation is

Applying the exponential martingale inequality of [18, Theorem 7.4], one gets

wherein $0 < \epsilon < 1$, $k$ is a random positive number. Utilizing the Borel-Cantelli's lemma [18, Lemma 2.4] for every $\omega \in \Omega$, there exists a integer $k_{0}$ such that all $k\geq k_{0}$,

Combining (3.2) and (3.3) for all $t \in (k-1, k]$, and letting $\epsilon \rightarrow 0$, one has

If $k\rightarrow +\infty$, then $t\rightarrow +\infty$, one obtains that $\frac{\ln I(0)}{t}\to 0$, $\frac{\ln k}{k-1}\to 0$. By combining the dominated convergence theorem and the ergodic properties of Markov chains, one can gain

Combining (3.4) and $\tilde{\mathcal {R}}_0 < 1$, we have

which shows $\lim_{t\to \infty}I(t) = 0$. □

Remark 3.1. In Theorem 3.1, we utilize the property of $S(t)\leq N(t)$ almost surely during stochastic perturbations. When the external environment remains stable and the transmission coefficient does not undergo sudden changes $N = 1$, $\mathcal R_0$ will be a threshold result common to deterministic system (2.1) and stochastic system (2.5). Based on $\tilde{\mathcal {R}}_0$, the measures of isolating close contacts and wearing masks can effectively reduce transmission rate and curb the spread of pertussis.

4.   Ergodic stationary distribution

Whether a disease will persist is a hot topic of current research ^{[23,24,25,36,37,38]}. Based on the Has'minskii theorem ^[22], an endemic equilibrium does not exist in the stochastic epidemic model, but an ESD can still exist. The ESD reflects the persistence of disease transmission within the community ^[16,39].

Denote $\overline\mu_{1} = \mu+\sum_{k = 1}^{N}\frac{\pi_{k}\sigma_{12}^2(k)}{2}+2\sqrt{\mu\sum_{k = 1}^{N}\pi_{k}\sigma_{11}(k)\sigma_{12}(k)} +2\sqrt[3]{\mu^2\sum_{k = 1}^N\pi_{k}\sigma_{11}^2(k)}$, $\mu_{2} = \mu+\gamma+\sum_{k = 1}^{N}\frac{\pi_{k}\sigma_{22}^2(k)}{2}$, $\overline\mu_{3} = \mu+\eta+\sum_{k = 1}^{N}\frac{\pi_{k}\sigma_{32}^2(k)}{2}$, $\mu_{4} = \mu+\eta+\sum_{k = 1}^{N}\frac{\pi_{k}\sigma_{42}^2(k)}{2}$, $\mathcal{R}_0^e = \frac{\hat\beta {S^e}}{\mu_{2}}$.

Theorem 4.1. Suppose that $\mathcal {R}_0^e > 1$. Then the solution $(S(t), I(t), W(t), R(t), r(t))$ of system (2.5) has a distinct invariant distribution regardless of the initial state $(S(0), I(0), W(0), R(0), r(0))\in \mathbb R_{+}^4\times \mathcal N$.

Proof. The diffusion matrix related to system (2.5) is illustrated by

Opt for $\tilde{M} = \min_{(S, I, W, R)\in \tilde D_{\epsilon}\subset R_{+}^4}\{(\sigma_{11}(k)S^2+\sigma_{12}(k)S)^2, \; (\sigma_{21}(k)E^2+\sigma_{22}(k)E)^2, \; (\sigma_{31}(k)I^2+\sigma_{32}(k)I)^2, \; (\sigma_{41}(k)A^2+\sigma_{42}(k)A)^2\} > 0$, ensuring

wherein $\zeta: = (\zeta_{1}, \zeta_{2}, \zeta_{3}, \zeta_{4})$. Inspired by ^[25], we will employ the generalized $\theta$-stochastic criterion method. To verify condition (ⅱ), we have the following four steps.

Step 1. For all $\tau\in(0, 1)$, here are some continuous $\tau$-stochastic Lyapunov functions to consider:

We observe that $\mu_{i}$ (i = 1, $\ldots$, 4) are all functions of the given variable $\tau\in(0, 1)$ that increase monotonically. Furthermore, it is found that

It should be pointed out that $\lim_{\tau\rightarrow {0^+}}S^e = \overline {S^e}$.

Step 2. To counteract the effects of nonlinear perturbations and Markov switching, we propose the following Lyapunov functions:

where $\zeta_{1}, \zeta_{2}, \zeta_{3}$, and $\xi_{1}, \xi_{2}$ are determined by (4.2), (4.10), and (4.5). $\omega_{i}(k), i = 1, 2, 3, 4$ are determined subsequently.

Utilizing the It$\hat o's$ formula to $V_{1}$, combining

and $(a+b)^2\leq 2(a^2+b^2)$, we get

Choosing

we gain

Let $Z_{1}(k): = \frac{\sigma_{12}^2(k)}{2}+\frac{2\mu}{(1-\tau)\xi_{1}}+\frac{\xi_{1}\sigma_{11}(k)\sigma_{12}(k)}{2}+\frac{8\mu}{3(1-\tau)\xi_{2}}+\frac{\xi_{2}^2\sigma_{11}^2(k)}{6}$. Motivated by ^[38], due to $\Gamma$ being irreducible, for $\overrightarrow Z_{1} = (Z_{1}(1), Z_{1}(2), \ldots, Z_{1}(N))^{\tau}$, we can get a vector $\overrightarrow \omega_{1} = (\omega_{1}(1), \omega_{1}(2), \ldots, \omega_{1}(N))^{\tau}$ so that the following Poisson system $\Gamma\overrightarrow \omega_{1} = \sum_{l = 1}^N\pi_{l}Z_{1}(l)-\overrightarrow Z_{1}$ is true, which indicates

Combining (4.3) and (4.4), it is clear that

In order to ensure $\frac{2\mu}{(1-\tau)\xi_{1}}+\frac{\xi_{1}\sum_{k = 1}^{N}\pi_{k}\sigma_{11}(k)\sigma_{12}(k)}{2}+\frac{8\mu}{3(1-\tau)\xi_{2}}+\frac{\xi_{2}^2\sum_{k = 1}^{N}\pi_{k}\sigma_{11}^2(k)}{6}$ reaches a minimum value, we opt for

We obtain

Employing the It$\hat o$'s formula to $V_{2}$ and referring to the similar methods described earlier, we get

where $Z_{2}(k) = \mu+\gamma+\frac{\sigma_{22}^2(k)}{2}$, a vector $\overrightarrow\omega_{2} = (\omega_{2}(1), \omega_{2}(2), \ldots, \omega_{2}(N))^{\tau}$ is chosen to meet the Poisson system $\Gamma\overrightarrow\omega_{2} = \sum_{l = 1}^{N}\pi_{l}Z_{2}(l)-\overrightarrow Z_{2}$, where$\overrightarrow Z_{2} = (\overrightarrow Z_{2}(1), \overrightarrow Z_{2}(2), \ldots, \overrightarrow Z_{2}(N))^{\tau}$. Therefore

In light of the similar steps mentioned in (4.7), we can derive

We select

and then derive

Let $Z_{3}: = \mu+\eta+\frac{\sigma_{32}^2(k)}{2}+\frac{\tau^2\sigma_{31}^2(k)}{6}$, a vector $\overrightarrow\omega_{3} = (\omega_{3}(1), \omega_{3}(2), \cdots, \omega_{3}(N))^\tau$ can be found that satisfies the Poisson system $\Gamma\overrightarrow\omega_{3} = \sum_{l = 1}^{N}\pi_{l}Z_{3}(l)-\overrightarrow Z_{3}$, wherein $\overrightarrow Z_{3} = (Z_{3}(1), Z_{3}(2), \cdots, Z_{3}(N))^{\tau}$. Hence

Applying the It$\hat o$'s formula to $V_{4}$, we can get

Define $Z_{4} = \mu+\eta+\frac{\sigma_{42}^2(k)}{2}$. Using the same method mentioned in (4.4), one can obtain

Step 3. Define two Lyapunov functions $Y_{1}, Y_{2}:\mathcal R_{+}^5\times \mathcal N \rightarrow \mathbb R$ that play important roles

wherein the constants $(\alpha_{1}, \alpha_{2}, \alpha_{3})$ and $\phi$, which are positive, are derived from (4.25) and (4.28), respectively. Suppose $S^e, I^e, W^e, R^e$ satisfy the following equations:

Let $I^e = 1$, meanwhile, we can get

Take a standard transformation of variables $(S, I, W, R)$

Considering that

From the inequality $\ln a\leq a-1\; \; (\forall a > 0)$, using the formulas (4.6), (4.16), and (4.17), we get

Using the same method as above, together with formula (4.8) and (4.18), we can obtain

Based on (4.12) and (4.16), we derive

According to (4.14) and (4.16), we have

Substituting (4.19)–(4.22) into $Y_{1}$, we obtain

Let

We derive through calculation

For the sake of convenience, we define $\chi_{1} = \frac{\check\sigma_{21}\check\sigma_{22}}{S^e}+\alpha_{1}\check\beta S^e+\alpha_{2}\kappa\check\beta W^e, \; \chi_{2} = \alpha_{2}\check\sigma_{31}\check\sigma_{32}W^e, \; \chi_{3} = \alpha_{3}\check\sigma_{41}\check\sigma_{42}R^e+\alpha_{2}\frac{8\eta W^e}{3\tau(1-\tau)}$. Combining (4.23) and (4.24), we get

By (4.15) and (4.26), we have

wherein $\phi$ fulfills the following two equations:

Substituting (4.27) into (4.28), we obtain

Step 4. Define the following two functions $Y_{3}$, $Y_{4}$:

where

Constructing the following form

wherein $M_{0} > 0$ is large enough to fulfill the following inequality:

We define the following compact set:

where $\epsilon < 1$ is a small enough positive constant that satisfies the following inequalities:

where $Q_{2} = \sup_{(I, R)\in\mathbb R_{+}^2}\{M_{0}[(\chi_{1}+\phi\gamma)I+\frac{\check\sigma^2_{41}}{2}R^eR^2+\frac{\check\sigma^2_{21}}{2S^e}I^2]\}$. Next, we partition $\mathbb R _{+}^4\backslash\mathcal D_{\epsilon}$ into eight domains.

It is clear that $\mathcal D_{\epsilon}^c = \cup_{i = 1}^8\mathcal D_{i}^c$. After this, to demonstrate $\mathcal LV\leq-1$.

Case 1. Suppose $(S, I, W, R)\in \mathcal D_{1}^c$, in consideration of (4.32) and (4.33), we gain

Case 2. If $(S, I, W, R)\in \mathcal D_{2}^c$, combining (4.32) and (4.33), we have

Case 3. If $(S, I, W, R)\in \mathcal D_{3}^c$, in consideration of (4.32) and (4.34), we obtain

Case 4. If $(S, I, W, R)\in \mathcal D_{4}^c$, given (4.32) and (4.35), we get

Case 5. If $(S, I, W, R)\in \mathcal D_{5}^c$, by (4.32) and (4.36), we gain

Case 6. If $(S, I, W, R)\in \mathcal D_{6}^c$, combining (4.32) and (4.36), one gets

Case 7. If $(S, I, W, R)\in \mathcal D_{7}^c$, in view of (4.32) and (4.36), we acquire

Case 8. If $(S, I, W, R)\in \mathcal D_{8}^c$, based on (4.32) and (4.36), we have

To summarize, we establish that $\mathcal LV\leq-1$ for all $(S, I, W, R)\in\mathcal D_{\epsilon}^c$. Thus, the second condition of Has'minskii theorem is met. Given that $\mathcal R_{0}^e > 1$, the solution of system (2.5) possesses a unique ESD $\pi(\cdot)$. □

Remark 4.1. We can calculate $\mathcal R_{0}^e \leq \mathcal R_{0}$. Theorem 4.1 implies that as the noise intensity approaches zero, there is not infection and $N = 1$, system (2.5) transforms into system (2.1), and the basic production number $\mathcal R_{0}$ serves as a common threshold. We obtain that when the noise intensity is relatively low, pertussis still prevails. Integrating the aforementioned Theorems 4.1 and 3.1, it is sensible to strengthen external intervention, wear masks, and appropriately control the number of students in the classroom to mitigate the spread of pertussis.

5.   Numerical simulations

In this section, numerical simulations are utilized to validate our conclusions. Using Milstein's higher-order method ^[40], the discretization equation associated with system (2.5) is given by:

of which $k\in\mathcal N$, $\Delta t > 0$ represents the size of a single iteration step and random variables $\xi_{j}$, $\eta_{j}$, $\zeta_{j}$ and $\vartheta_{j}$ are four independent variables that adhere to a Gaussian distribution $\mathbb N(0, 1)$ for $j = 1, 2, \ldots, n$, respectively. We consider $r(t)$ as a right-continuous Markov chain defined on the state space $\mathcal N = \{1, 2\}$, with the corresponding generator

By solving $\pi \Gamma = 0$ and $\pi_1+\pi_2 = 1$, the stationary distribution of $\Gamma$ is obtained as $\pi = (\pi_1, \pi_2) = (\frac{5}{9}, \frac{4}{9})$. The initial values are chosen as $(S(0), I(0), W(0), R(0)) = (0.45, 0.35, 0.15, 0.05)$ and the primary parameters in system (2.5) are set as follows:

It is worth noting that, we estimate the parameters of model (2.1) by using the least squares method, and take the year as the unit of measurement. Based on ^[2] and ^[10], $\mu = 0.02$, $\nu = 0.8$, $\eta = 0.16$ can be chosen, respectively. Then, by analyzing the data of pertussis infections from $2006$ to $2019$ (selected from the National Health Commission of the People's Republic of China, see Figure 1), $\kappa = 0.03$ and $\gamma = 0.5$ have been selected. In Figure 1, the susceptible initial value is set to $1030$, the infectious initial value is set to $2100$, the waning initial value is set to $3000$, and the recovered initial value is set to $3000$.

By direct calculation, $\mathcal R_{0} = 3.2004 > 1$. We focus on three aspects, as follows:

● The system (2.5) shows both a stationary distribution and ergodic behavior when $\mathcal R_0^e > 1$.

● The dynamical characteristics of system (2.5) for $\mathcal R_0^e < 1$.

● The effect of stochastic perturbations on the disease dynamics within system (2.5).

Example 5.1. We opt for the transmission rates $(\beta_{1}, \beta_{2}) = (1, 2)$ and the stochastic noises $(\sigma_{ij}(1), \sigma_{ij}(2)) = (2 \times 10^{-4}, 6 \times 10^{-4})$ for any $i = 1, 2, 3, 4; j = 1, 2,$ respectively. The calculations show that $\mathcal R_0^e = 1.6654 > 1$. Moreover, since $\mathcal R_0 > 1$, we can obtain the existence of an endemic equilibrium for deterministic system (2.1). Figure 2 illustrates that the stochastic system (2.5) admits a unique global positive solution, which follows a unique ESD $\pi(\cdot)$, signifying the persistence of the disease in a community. On the other hand, we notice a slight difference between the stochastic solution and the solution without randomness. Meanwhile, Figures 3 and 4 depict the corresponding fitting density function and frequency histogram, as well as the movement of Markov chain, respectively.

Example 5.2. Let the noises $(\sigma_{ij}(1), \sigma_{ij}(2)) = (2 \times 10^{-1}, 6 \times 10^{-1})$ for each $i = 1, 2, 3, 4;$ $j = 1, 2$. By means of direct derivation, we get $\mathcal R_{0}^e = 0.4512 < 1$. Therefore, we cannot determine the existence of an ESD of system (2.5). As shown in Figure 5, we can obtain that disease of system (2.5) will become extinct in the long run. Figure 6 shows the corresponding fitted density functions.

Example 5.3. To demonstrate the influence of stochastic perturbations, we have selected three distinct noise levels $\sigma_{1} = (\sigma_{ij}(1), \sigma_{ij}(2)) = (0.002, 0.006)$, $\sigma_{2} = (\sigma_{ij}(1), \sigma_{ij}(2)) = (0.040, 0.080)$, $\sigma_{3} = (\sigma_{ij}(1), \sigma_{ij}(2)) = (0.300, 0.500)$, wherein $i = 1, 2, 3, 4, j = 1, 2$. Figure 7 shows the trajectories of $S(t), I(t), W(t), R(t)$, respectively. It is evident that the values of $S(t), I(t)$, $W(t), R(t)$ exhibit significant volatility as the noise intensity increases. Moreover, combining Figures 2, 5, and 7, we can observe that larger perturbations will drive the disease to extinction, whereas smaller perturbations may lead to its persistence within the community.

Example 5.4. (Sensitivity analysis) Under the conditions of satisfying 3.1 and 4.1, we can find that the basic reproduction number $\mathcal R_0$ is a crucial value for the eradication and persistence of pertussis disease. Therefore, Figure 8 shows sensitivity analysis on the reproductive number $\mathcal R_0$ of the disease. It is evident that transmission rates $\beta$ and loss of immunity rate $\eta$ are positively correlated with $\mathcal R_0$. Meanwhile, there is a negative correlation between recovery rate $\gamma$, vaccination probabilities $\nu$, death rate $\mu$, and $\mathcal R_0$. Physical exercise to enhance immunity and increase vaccination rates are beneficial for controlling the spread of pertussis.

6.   Conclusions

In this paper, we primarily focus on the dynamics of pertussis in nature. By integrating natural immune enhancement and long-term immune waning, a stochastic pertussis system with Markov switching (2.5) is constructed and analyzed. Following the ideas of papers ^[24,25], the existence and uniqueness of the global positive solution are proven. Next, the exponential extinction of pertussis is shown in Theorem 3.1. Moreover, inspired by ^[25], $\theta-$stochastic Lyapunov functions are used to establish the sufficient condition $\mathcal R_{0}^e > 1$, which guarantees the existence and uniqueness of an ESD for system (2.5) in Theorem 4.1. Finally, several examples are provided to confirm our theoretical findings (Figures 2–7). Based on these examples, we infer that stochastic noise is negatively correlated with the persistence of the disease. Furthermore, some questions deserve further study, for example, considering a stochastic pertussis model with semi-Markovian switching and impulses. The density function of a stochastic delay pertussis model with semi-Markov switching and Ornstein-Uhlenbeck process ^[41] may also be analyzed.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research was supported by the National Natural Science Foundation of China (Nos.12161068, 12261069, 12301630).

Conflict of interest

The authors declare there is no conflicts of interest.