Asymptotic behavior of a stochastic hybrid SIQRS model with vertical transmission and nonlinear incidence

1.   Introduction

Every outbreak of infectious diseases will endanger people's lives and have extremely negative impacts on social economy. For example, the COVID-19 broke out in 2019 and swept the world. By May 3, 2023, the number of confirmed cases in the world was close to 765 million, of which more than 6.92 million people died due to infection ^[1]. Because of the fast mutation of the virus and the inability to develop effective drugs in a timely manner, isolation is considered a valid approach to reduce the spread of contagious diseases ^[2]. People can take practical and feasible measures such as self isolation to prevent the rapid spread of the disease, so as to reduce the pressure of the medical department and cause the infected people obtain effective treatment. However, the policy of isolation has also had some negative effects, which will lead to recession in the economy to varying degrees and unemployment due to reduced demand for labor. Hence, it is an interesting topic to find isolation strategies to prevent the further development of infectious diseases and minimize the negative impact of the epidemics.

The epidemic model plays a crucial role in countering infectious diseases. It incorporates various factors that affect the spread of diseases into the dynamic system to deeply understand the transmission mechanism of diseases, so that the potential impact of different factors can be better evaluated. Isolation is a very effective measure of controlling disease ^[3,4,5,6,7]. In ^[7], the authors have discussed a SIQRS model with isolation as the following form,

where $S_t$, $I_t$, $Q_t$ and $R_t$ represent the numbers of the susceptible, the infected, isolated, and the removed people, respectively. $\Lambda$ denotes the recruitment rate due to immigration, $\mu$ is the natural mortality, $\beta$ stands for the disease transmission rate, $r_1$ shows the rate of the recovered who lost immunity and returns to the susceptible, $\alpha_1$ and $\alpha_2$ express the cure rates of the infected $I_t$, and the isolated $Q_t$ to $R_t$. $\delta$ indicates the isolation rate of the infected. $\gamma_1$ and $\gamma_2$ mean the disease-caused mortality in classes $I_t$ and $Q_t$. These parameters are assumed to be positive. To facilitate writing, let $a_1 = \alpha_1+\delta+\mu+\gamma_1$ and $a_2 = \mu+\alpha_2+\gamma_2$.

Various systems in life, including infectious diseases, will inevitably be disturbed by stochastic factors, which will alter the trajectory of the system more or less. Thus, the epidemic models with stochastic factors have been discussed widely due to their more applicability and richer research contents ^{[4,8,9,10,11,12]}. White noise characterized by Brownian motion is a common stochastic disturbance, which is often introduced into infectious disease models. In addition, color noise represented by Markov chain is another important stochastic factor, which can portray the switching between different environments, states, or temperatures ^{[12,13,14,15]}. In this paper, we will discuss the above epidemic model containing these two types of stochastic noise.

In epidemiology, the incidence rate shows the cases of second-generation infected persons per unit time. In many literatures, the incidence function used in epidemic models is the bilinear function denoted by $\beta SI$ ^[11,16,17]. This function is based on the fact that the population is evenly mixed and everyone is equally likely to be infected. Owing to this assumption, the nonlinear incidence rates have a wider application and have attracted a large number of scholars to study ^{[10,12,14,18,19,20]}. The authors have discussed a SIQS model with the incidence rate $\frac{\beta SI}{1+rI}$ and obtained a value $R_S$ to determine the extinction and persistence of the model ^[19]. Guo-Luo have investigated a hybrid SIR model with Beddington-DeAngelis function ^[14] and the authors have studied the epidemic model with the incidence rate $\beta f(S)g(I)$ ^[10]. When incorporating random factors above and the nonlinear incidence rate, (1.1) becomes

where $W_t^{(j)}, \; j = 1, 2, 3, 4$ are mutually independent Brownian motion defined on the complete space $(\Omega, \mathfrak{F}, \{\mathfrak{F}\}_t, \mathbb{P})$. $\sigma_j(l)$, $j = 1, \cdots, 4$, $l = 1, \cdots, M$ express the intensities of stochastic disturbances and $\{\theta_t\}_{t\ge0}$ denotes the continuous time Markov chain, which is independent of $W_t^{(l)}$, taking values in the state space $\mathcal{S} = \{1, 2\cdots, M\}$ and the generator $Q = (\gamma_{ij})_{M\times M}$ satisfies

for $\varepsilon \downarrow 0$. $\gamma_{jl} > 0$ for $j\neq l$ and $\sum_{l = 1}^M\gamma_{jl} = 0$ for any $j\in \mathcal{S}$. The general incidence function $G(S, I, \theta)$ has the following assumption:

Assumption 1. For the variables $S$ and $I$, the function $G(S, I, \theta)$ is locally Lipschitz continuous. For each $l\in \mathcal{S}$, $G(S, I, l)$ is non-increasing in $I$ and non-decreasing in $S$ with $G(0, I, l) = 0$. Moreover, the function $G$ is continuous uniformly at $I = 0$, that is

Assume further that there exist positive constants $c(l)$ and $c_1(l)$ such that $\frac{\partial G(S, I, l)}{\partial S}\leq c(l)$ and $\frac{\partial G(S, I, l)}{\partial I}\ge -c_1(l)$ for any $I$ and $l\in \mathcal{N}$. Therefore, $G(S, I, l)\leq c(l)S$ holds due to $G(0, I, l) = 0$.

For the incidence function $G(S, I, l)I$ above, it contains many types that appear in other literature, such as the bilinear form $\beta(l)S{I}$, saturated rate $\frac{\beta(l)S{I}}{1+aI}$, the rate $\frac{\beta(l)S{I}}{m+S}$, Beddington-DeAngelis rate $\frac{\beta(l)S{I}}{1+m_1(l)S+m_2(l)I}$ and other forms.

As we know, besides the contact spread of disease, there is also a vertical transmission, in which the disease is transmitted from the infected mother to the newborn. Vertical transmission is considered as an important mode of AIDS transmission. Therefore, many scholars have discussed the epidemic models introducing the vertical transmission ^{[21,22,23,24]}. The authors in ^[24] have concerned a SIR model with the birth rate $b$ and vertical transmission rate $p$ from the infected mother. We utilize these symbols to express the same meanings. Assume that the newborns of the classes $S$, $Q$, $R$ all become susceptible and $\mu > b$ in this paper.

In view of the above discussion, we study the following stochastic hybrid SIQRS model with nonlinear incidence rate and vertical transmission

These factors can reflect different aspects of actual problems and increase the difficulty of study. Since the term $R_t$ appears in the first equation in (1.4), we cannot omit $R_t$ to reduce the dimension of the model as in ^[9] and need to study the system with four components. This paper is constructed as follows: Section 2 gives the existence and uniqueness of positive solutions of model (1.4) and some properties, which are used later. Section 3 presents the major results of the paper, that is, we obtain a threshold which can be used to decide the extinction of model (1.4) and the existence of invariant measure. Section 4 aims to prove the major results of Theorem 3.1, and Section 5 provides some remarks and compares our results with those of other studies. Section 6 constructs some examples and presents numerical simulations to test the results. Section 7 summarizes this article.

2.   Preliminaries

In this paper, $\mathbb{R}_+^4: = \{(a_1, a_2, a_3, a_4)|a_k\ge0, k = 1, 2, 3, 4\}$ and $\mathbb{R}_+^{4, o}: = \{(a_1, a_2, a_3, a_4)|a_k > 0$, $k = 1, 2, 3, 4\}$. $\mathbb{E}_{s, i, q, r, l}$ represents the expectation and $\mathbb{P}_{s, i, q, r, l}$ indicates the probability with initial value $(s, i, q, r, l)\in \mathbb{R}_+^4\times \mathcal{S}$. Assume that $\check{\alpha}: = \max_{l\in \mathcal{S}}\{\alpha(l)\}$, $\hat{\alpha}: = \min_{l\in \mathcal{S}}\{\alpha(l)\}$ and $\alpha_1\vee \alpha_2 = \max\{\alpha_1, \alpha_2\}$. Similar symbols for other variables are defined identically. Take into account the general hybrid stochastic differential equations(short for SDEs),

For the function $V(X_t, \theta_t)$, the operator $\mathcal{L}V(X, l)$ is defined by

Then the generalized It$\rm\hat{o}$'s formula is presented as

We recommend the Theorem 1.45 in ^[25] to grasp the details on the measure $\mu(ds, dl)$ and the function $\nu$.

We are going to lay out the following theorem to get the properties of the solution to model (1.4).

Theorem 2.1. For any initial condition $(S_0, I_0, Q_0, R_0, \theta_0)\in \mathbb{R}_+^4\times \mathcal{S}$ in (1.4), the following statements hold true: (1) model (1.4) has the unique solution $(S_t, I_t, Q_t, R_t, \theta_t)$, which stay in $\mathbb{R}^4_+\times \mathcal{S}$ with probability 1. In addition, the five-component solution $(S_t, I_t, Q_t, R_t, \theta_t)$ is the Markov-Feller process. (2) For any $0 < \alpha < \vartheta < 1$, there exist constants $A_1 > 0$ and $A_2 > 0$ satisfying

Proof. The solution must satisfy the changing characteristics of the model (1.4). We primarily pay attention to (2) due to (1) is analogous to the proof of Theorem 2.2 in ^[26]. Construct the function $V_1(S, I, Q, R): = (S+I+Q+R)^{1+\vartheta}+S^{-\alpha}$ and $\sigma(l): = \max_{i = 1, 2, 3, 4}\{\sigma_i(l)\}$, then it has

Choose sufficiently small $\alpha > 0$ such that $A_1: = \min_{l\in \mathcal{S}}\{\mu(l)-b(l)-\frac{\alpha\sigma^2(l)}{2}\} > 0$, $\forall l\in \mathcal{S}$. Because

for $0 < \alpha < \alpha_3 < \vartheta < 1$, we have

Hence, owing to $\frac{\alpha}{\alpha_3} < 1$ and $\alpha_3 < \vartheta$, it yields $\mathcal{L}V_1(S, I, Q, R)+A_1V_1(S, I, Q, R)\leq A_2$, where

Calculating $e^{A_1t}V_1(S+I+R+Q)$ by the It$\rm\hat{o}$'s formula leads to

Integrating from $0$ to $t$ and taking expectation, it has

This proves the assertion. □

3.   Main results

In this section, we will present the main conclusions of this paper. Before this, we briefly discuss the generation of the threshold of disease extinction in model (1.4).

Consider the first equation in model (1.4) on the boundary $I_t = 0$, $Q_t = 0$ and $R_t = 0$, we have

For the initial value $s$ of Eq (3.1), let $\widetilde{S}_t^s$ be its solution. Direct calculation to the non-negative function $\widetilde{S}-\ln\widetilde{S}-1$ and exploiting the results in ^[27] say that non-degenerate system (3.1) is positive recurrent, thus, the unique invariant measure $\chi_0(\cdot, \cdot)$ for (3.1) on $[0, \infty)\times \mathcal{S}$ satisfying $\chi_0([0, \infty), \mathcal{S}) = 1$ exists. Moreover, the stationary distribution $\pi$ of $\{\theta_t\}_{t\ge0}$ is the marginal distribution of $\chi_0(\cdot, \cdot)$. Due to Theorem 2.1, it has

Hence, the value

is well-defined.

Using the It$\rm\hat{o}$'s formula to $\ln I_t$ and dividing by $t$, one has

If $\limsup_{t\to\infty}\frac{\ln I_t}{t} < 0$, then $\lim_{t\to\infty}I_t = 0$. Using the Fatou lemma implies $\lim_{t\to\infty}Q_t = 0$ and $\lim_{t\to\infty}R_t = 0$. Thus, for $t$ sufficiently large, $I_t\approx 0$ and $S_t$ will approach $\widetilde{S}_t$ on the boundary, then

and $\limsup_{t\to\infty}\frac{\ln I_t}{t}$ will be near to the threshold $\Delta$.

Sketchily, when $\Delta < 0$, for the initial condition $(s, i, q, r, l)$ with small enough $i$, it yields $\limsup_{t\to\infty}\frac{\ln I_t}{t}\approx \Delta < 0$, that is, the disease will die out. Conversely, when $\Delta > 0$, $\limsup_{t\to\infty}\frac{\ln I_t}{t}\approx \Delta > 0$ will let $I_t$ be not small in the long term. This procedure seems simple, but the strict proof is not simple and needs scrupulous treatment.

Now, we present our main conclusions, in which $\Delta$ will be proved to distinguish different behaviors of disease in model (1.4).

Theorem 3.1. For $\Delta$ in (3.2), we have

(1). When $\Delta < 0$, the solution $(S_t, I_t, Q_t, R_t, \theta_t)$ of model (1.4) with initial condition $(s, i, q, r, l)\in \mathbb{R}^{4, o}_+\times \mathcal{S}$ has that

which means the disease will become extinct in exponential form with rate $\Delta$.

(2). For model (1.4), when $\Delta > 0$, there exists a constant $\nabla > 0$ such that

which signifies the disease $I_t$ is persistent in the mean and the model has the unique invariant measure $\chi^*$.

4.   The proof of Theorem 3.1

In this section, we will prove the two conclusions of Theorem 3.1, and implement them separately in two subsections.

4.1.   The proof of Part 1 of Theorem 3.1

We shall first prove the Part 1 in Theorem 3.1 in this subsection. Let's start with the following lemma.

Lemma 4.1. If $\Delta < 0$, for any $K > 0$ and $\epsilon > 0$, the constant $k_1 > 0$ can be found so that for any initial value $(s, i, q, r, l)\in[0, K]\times [0, k_1]^3\times \mathcal{S}$(where $[0, k_1]^3$ denotes $[0, k_1]\times [0, k_1]\times [0, k_1]$), it yields

Proof. The idea of this proof is that when the initial values of $I_t$, $Q_t$, and $R_t$ are all very small and under the condition of $\Delta < 0$, $I_t$, $Q_t$, and $R_t$ shall always be small enough. Define a constant $\nu$ by $\nu = \min\{-\Delta, \hat{a}_2+\frac{\hat{\sigma}_3^2}{2}, \hat{\mu}+\hat{r}_1+\frac{\hat{\sigma}_4^2}{2}\}$. Thus, $\nu > 0$. For $\widetilde{S}_t$ in (3.1), due to the existence of terms $qb(\theta_t)I$, $b(\theta_t)Q$ and $(r_1(\theta_t)+b(\theta_t))R$ in the first equation of (1.4), we can't use the comparison theorem to get $S_t\leq \widetilde{S}_t$ with the same initial condition. Take into account the equation

Let $\widetilde{S}^{(k)}_t$ be the solution of (4.2) with initial condition $s\in[0, K]$. Similar to (3.1), (4.2) admits the unique invariant measure denoted by $\chi_k$. Lemma 3.1 in ^[28] says that there is $k_0$ satisfying $\widetilde{\Delta}\leq \Delta+\frac{\nu}{9}$ with

Due to (1.3), for any $s > 0$, $l\in \mathcal{S}$ and $0\leq i\leq k_0$, it has

Consider (4.2) with $k$ replaced by $k_0$ above, by virtue of the ergodicity of $\widetilde{S}^{(k_0)}_t$, one has

Thus, there exist constants sufficiently small $\epsilon > 0$ and $T_1 > 0$ so that $\mathbb{P}_{K, l}(\Omega_1)\ge 1-\frac{\epsilon}{4}$ for $\forall t\ge T_1$, where

Here $(K, l)$ in symbol $\mathbb{P}_{K, l}$ denotes the initial condition of (4.2) with $k_0$. For $s\leq K$, $t\ge0$, it has $\widetilde{S}^{(k_0)}_{s, l}(t)\leq \widetilde{S}^{(k_0)}_{K, l}(t)$ by the uniqueness of solution. This makes $\mathbb{P}_{s, l}(\Omega_1)\ge 1-\frac{\epsilon}{4}$.

Assume $M^{(k)}_t: = \int_0^t\sigma_k(\theta_u) dW_u^{(k)}$, $k = 1, 2, 3, 4$. According to $\lim_{t\to\infty}\frac{M^{(k)}_t}{t} = 0, \; a.s.$, then there is $T_2 > 0$ so that $\forall t\ge T_2$, $\mathbb{P}(\Omega_2)\ge 1-\frac{\epsilon}{4}$ with

Assume $T = \max\{T_1, T_2\}$, let

and

then Theorem 2.1 and the fact that $G(S, I, l)\leq \check{c}S$ lead to $\mathbb{P}(\Omega_3)\ge 1-\frac{\epsilon}{4}$ and $\mathbb{P}(\Omega_4)\ge 1-\frac{\epsilon}{4}$ for a sufficiently large $M_1$.

Let $C_1: = e^{M_1}+\check{\delta}e^{4M_1}T$, choose the constant $k_1 > 0$ to be small enough so that

where the constants $C_2 > 0$, $C_3 > 0$ will be found in (4.17) and (4.20). Define a stopping time $\tau_1$ as

By the expressions of $I_t$, $Q_t$, $R_t$ in (1.4) with initial data $(I_0, Q_0, R_0) = (i, q, r)\in[0, k_1]^3$, using the method of constant variation yields

and

where $\Upsilon_1(t) = e^{-\int_0^t(a_2(\theta_v)+\frac{\sigma_3^2(\theta_v)}{2})dv+\int_0^t\sigma_3(\theta_v)dW_v^{(3)}}$, $\Upsilon_2(t) = e^{-\int_0^t(\mu(\theta_v)+r_1(\theta_v)+\frac{\sigma_4^2(\theta_v)}{2})dv+\int_0^t\sigma_4(\theta_v)dW_v^{(4)}}$.

Therefore, by virtue of (4.8), we get with $\omega\in \Omega_3\cap\Omega_4$ and $t\in [0, T]$ that

For $t\leq T$, the expressions of $\Upsilon_1(t)$ and $\Upsilon_2(t)$ with $\omega\in \Omega_4$ result in

Using the results above and (4.9), we get

In addition,

Hence, for almost every $\omega\in \cap_{i = 3}^4\Omega_i$ and $t\leq T$, (4.7) and (4.11)–(4.13) can deduce that $\max\{I_t, Q_t, R_t\} < k_0$, which implies $T < \tau_1$.

$\tau_1 = \infty$ will be proved next for almost every $\omega\in \cap_{l = 1}^4\Omega_l$.

Observe that for $S_t = \widetilde{S}_0^{(k_0)} = s$ in (1.4) and the Eq (4.2) with $k$ replaced by $k_0$, $S_t\leq \widetilde{S}_t^{(k_0)}$, $\forall t < \tau_1$ is established due to $\max\{I_t, Q_t, R_t\} < k_0$ and the comparison theorem. Thus, when $t\in[T, \tau_1)$ and almost every $\omega\in \cap_{l = 1}^4\Omega_l$, we obtain from (4.4) and (4.8) that

For $Q_t$ on $t\ge T$, (4.9) can be reorganized as

For the second term in (4.15), we get that

For the first term of (4.15), one has

This as well as (4.16) results in

for some constant $C_2 > 0$.

Now, consider $R_t$ in (4.10), one has

For the second expression in (4.18), we have

Similar to the first term in (4.15), we have from the first term in (4.18) that

This result, combined with (4.19), leads to that

for some constant $C_3 > 0$.

Let a positive integer $n_0 > T$. By virtue of (4.7), (4.14), (4.17) and (4.20), it easy to get that for $t\in [0, \tau_1\wedge n_0)$ and almost every $\omega\in \cap_{l = 1}^4\Omega_l$, $I_t\leq k_1(e^{2M_1}+1)\leq k_0$, $Q_t\leq k_1(C_1+C_2)\leq k_0$ and

Hence, $\max\{I_t, Q_t, R_t\}\leq k_0$ implies $\tau_1 > n_0$. Due to $n_0$ is arbitrary, we have $\tau_1 = \infty$, which means that $\lim_{t\to\infty}\frac{\ln I_t}{t}\leq -\frac{5\nu}{9} < 0$, $\lim_{t\to\infty}\frac{\ln Q_t}{t}\leq -\frac{\nu}{3} < 0$ and $\lim_{t\to\infty}\frac{\ln R_t}{t}\leq -\frac{\nu}{9} < 0$. It's easy to figure out that $\mathbb{P}(\cap_{l = 1}^4\Omega_l)\ge 1-\epsilon$. Therefore, (4.1) is proved. □

With Lemma 4.1, the following proof when $\Delta < 0$ is analogous to Section 2 of Theorem 2.2 in ^[28]. In this way, we have proved Part 1 of Theorem 3.1.

4.2.   The proof of Part 2 in Theorem 3.1

Next, we will prove Part 2 in Theorem 3.1. We first prove the persistence of the disease in (1.4) by taking advantage of a new way when $\Delta > 1$.

Let $\bar{c} = (c(1), \cdots, c(M))^T$($c(l)$ appears in Assumption 1 and $K = diag(\mu(1)-b(1), \mu(2)-b(2), \cdots, \mu(M)-b(M))$, take into account the equation $(K-Q)\eta = \bar{c}$, then it has a unique positive solution (Theorem 2.10 in ^[25]). Assume that $\eta = (\eta(1), \eta(2), \cdots, \eta(M))^T$ is its solution, then $(\mu(l)-b(l))\eta(l)-\sum_{j\in \mathcal{S}}\gamma_{lj}\eta(j) = c(l)$.

Let $V_2: = -\ln I$ and $V_3: = \textbf{1}_{\{\widetilde{S}\ge S\}}(\widetilde{S}- S)$(where $\textbf{1}$ denotes the indicator function), then direct calculation by the It$\rm\hat{o}$'s formula to $V_2+\eta(l)V_3$ and using the monotonicity of $G(S, I, l)$ at $S$ result in

Integrate for (4.21) and divide by $t$ as well as take the limit, then the ergodicity of $\widetilde{S}$ means

In this way, we have proved the persistence of the disease. Next, we shall prove that model (1.4) has an invariant probability measure.

Let $V_4 = V_2+\eta(l)V_3-\check{\eta}I$, $V_5 = \frac{1}{1+\alpha_4}(S+I+Q+R)^{1+\alpha_4}$, $V_6 = -\ln S-\ln Q-\ln R$ and $\bar{V} = H_1V_4+V_5+V_6$, where $H_1 > 0$ and $\alpha_4\in (0, 1)$ will be detailed later. The continuity of $\bar{V}$ leads to that there is a minimum value $\bar{V}_*$ such that $\widetilde{V} = \bar{V}-\bar{V}_*$ is non-negative.

Let $\sigma(l)$ and $A_1$ be the same as in Theorem 2.1, using the It$\rm\hat{o}$'s formula to $V_5$ and $V_6$, one has

and

Hence,

where

From the assumption that $G(S, I, l)$ is continuous uniformly at $I = 0$, hence, when $I$ is sufficiently small and $H_1$ is sufficiently large, it has

and $V_7\leq -1$.

Next, when $S\to 0^+$ or $Q\to 0^+$ or $R\to 0^+$, $V_7\leq -1$ can be obtained due to the terms $-\frac{\hat{\Lambda}}{S}$, $-\frac{\hat{\delta}I}{Q}$, $-\frac{\hat{\alpha}_1I}{R}$. Moreover, the term $-\frac{A_1}{2}(S+I+Q+R)^{\alpha_4+1}$ leads to $V_7\leq -1$ when $S(or I, Q, R)\to \infty$. The detailed proof process is similar to Theorem 4.1 in ^[29].

For the sufficiently small constant $\varepsilon$, define $D_\varepsilon = \{(S, I, Q, R, l)\in \mathbb{R}_+^4\times \mathcal{S}:\; \varepsilon\leq S\leq\frac{1}{\varepsilon}, \varepsilon\leq I\leq\frac{1}{\varepsilon}, \varepsilon^2\leq Q\leq\frac{1}{\varepsilon^2}, \varepsilon^2\leq R\leq\frac{1}{\varepsilon^2}\}$, it can be concluded from the above that $V_7\leq -1$ in $\mathbb{R}_+^4\times \mathcal{S} \backslash D_\varepsilon$.

Due to the compactness of the set $D_\varepsilon$ and continuity of the function $V_7(S, I, Q, R, l)$, we get that there exists a constant $K_3 > 0$ such that $V_7(S, I, Q, R, l)\leq K_3$ for $(S, I, Q, R, l)\in D_\varepsilon\times \mathcal{S}$. Therefore,

The ergodicity of $\widetilde{S}$ and $\theta_t$ reaches to

From the non-negativity of $\widetilde{V}$ and taking the limit, we get

which means

For the Markov-Feller process $(S_t, I_t, Q_t, R_t, \theta_t)$, (4.22) and the compactness of the set $D_\varepsilon$ results in the invariant probability measure marked as $\chi^*$ by virtue of Theorem 2 in ^[30].

5.   Some remarks on the results

In the previous sections, we have presented and proved the threshold

to determine the different properties of the model we have established. However, the value cannot be calculated obviously, and in this section we examine another form of this value for some specific incidence functions. Let $G(S, I, \theta_t) = \frac{\beta(\theta_t)S}{f(I)}$, where $f(I)$ is increasing as $I$ with $f(0) > 0$. The functions satisfying these conditions have the forms with $\beta(l) S$, $\frac{\beta(l) S}{1+aI}$, $\frac{\beta(l)S}{1+aI^2}$, etc. Then, $\Delta$ can be expressed as

because $\pi$ is the marginal distribution of $\chi_0(\cdot, \cdot)$. Let us focus on another form of the first term and provide the following remark first.

Remark 5.1. $\sum_{l\in \mathcal{S}}\int_{(0, \infty)}\frac{\beta(l)s}{f(0)}\chi_0(ds, l) = \sum_{l\in \mathcal{S}}\pi_l\varrho(l)\Lambda(l)$, where $\varrho = (\varrho(1), \cdots, \varrho(M))^T$ satisfies the equation $(diag(\mu(1)-b(1), \cdots, \mu(M)-b(M))-Q)\varrho = (\frac{\beta(1)}{f(0)}, \cdots, \frac{\beta(M)}{f(0)})^T$. So,

Let $R_1^S = \frac{\sum_{l\in \mathcal{S}}\pi_l\varrho(l)\Lambda(l)}{\sum_{l\in \mathcal{S}}\pi_l(a_1(l)-pb(l)+\frac{\sigma_2^2(l)}{2})}$, then $\Delta < 0$ is equivalent to $R_1^S < 1$.

Proof. It is easy to see the equation $(diag(\mu(1)-b(1), \cdots, \mu(M)-b(M))-Q)\varrho = (\frac{\beta(1)}{f(0)}, \cdots, \frac{\beta(M)}{f(0)})^T$ has the nonnegative solution $\varrho = (\varrho(1), \cdots, \varrho(M))^T$, whose proof is similar to that in Subsection 4.2. This implies $(\mu(l)-b(l))\varrho(l)-\sum_{j = 1}^M\gamma_{lj}\varrho(j) = \frac{\beta(l)}{f(0)}$.

For $\widetilde{S}_t$ in (3.1), let $V_8(l): = \varrho(l)\widetilde{S}$, then we obtain that

Thus,

Dividing by $t$, taking the limit and combining with the ergodicity of the Markov chain bring about

By virtue of the ergodicity of $(\widetilde{S}_t, \theta_t)$, one has

Hence,

We arrive at the Remark 5.1. □

In what follows, we will compare the results with those of other papers.

Remark 5.2. The authors have studied the SIQR model in ^[4] with the incidence rate $G = \frac{\beta {S}I}{1+\alpha I}$ and no Markovian switching, obtained the value $R_0^S = \frac{\beta \Lambda}{\mu(\mu+\gamma+\delta+\theta+\frac{\sigma_2^2}{2})}$ to distinguish the disease extinction or persistence, and provided another different value $\widehat{R}_0^S$ to derive the stationary distribution of the discussed model. Notice that in order to obtain different properties of the model, another condition $\mu > \frac{\max_{i = 1}^4\sigma_i^2}{2}$ is necessary. While in this paper with $f(0) = 1$ and $b = 0$, we get the same value $\Delta$($\Delta > 0$ is equivalent to $R_0^S > 1$ in ^[4]) to distinguish different dynamics of the model without additional conditions.

Remark 5.3. In ^[21], the authors have discussed a SIS model with vertical transmission and provided two values $R_0^S$, $\tilde{R_0^S}$ to determine different dynamics, that is, when $R_0^S > 1$, the model admits a stationary distribution and the disease will continue, while $\tilde{R_0^S} < 1$, the disease will die out. Obviously, there is a certain interval between the values that determine two different behaviors, and the two values are not the same. However, the value $\Delta$ in this paper can be used to judge different dynamics of the SIQRS model with Markovian switching and vertical transmission.

Remark 5.4. Liu has investigated a hybrid SIS model with the bilinear rate $G(S, I, \theta_t) = \beta(\theta_t)SI$ and no vertical transmission, obtained the value $\bar{R}_0$ to determine the ergodic stationary distribution and extinction ^[16]. Through the discussion of Remark 5.1, it can be seen that the threshold $R_1^S$ in this paper is identical to $\bar{R}_0$ in the model with the bilinear rate and no vertical transmission. Hence, ours can be regarded as the generalization of ^[16].

Remark 5.5. From the above analysis, it can be inferred that when $\Delta < 0$, the disease will tend to be extinction, how to take measures to let $\Delta < 0$ hold true so as to achieve the goal of disease control is a practical problem. By the expressions of $\Delta$ (or $R_1^S$), some feasible measures in practice are as follows: (i) When the epidemic is severe, medical forces should be increased to improve the cure rate, and the isolation rate can be increased to separate different populations and reduce mutual infection. (ii) When epidemics spread vertically, the vertical transmission should be reduced to control the disease. Women who are willing to have children should undergo testing or treatment. They may prepare for pregnancy when they are not infected, and when infected, the birth rate of newborns should be reduced.

6.   Numerical examples

We will list some examples and show their simulations to check the theoretical results.

Example 1. We first check the persistence and extinction of model (1.4) under Markovian switching. Let $(\theta_t)_{t\ge0}$ be the Markov chain with space $\mathcal{M} = \{1, 2\}$, the $Q-$ matrix is

So, the stationary distribution $\pi = (\frac{b}{a+b}, \frac{a}{a+b})$. Let the function $G(S, I, \theta) = \frac{\beta(\theta_t)S}{1+2I}$ and the initial values are $S_0 = 2.1$, $I_0 = 1$, $Q_0 = 0.2$, $R_0 = 1.2$, assume that the values of each parameter are as follows: $a = 2$, $b = 1$, $q = 0.6$, $\Lambda = [0.04, 0.02]$ (the two numbers represent the values of $\Lambda$ in two environments, and the followings are similar), $\mu = [0.03, 0.05]$, $\beta = [0.4, 0.15]$, $r_1 = [0.25, 0.35]$, $\alpha_1 = [0.45, 0.5]$, $\alpha_2 = [0.5, 0.6]$, $\delta = [0.45, 0.5]$, $\gamma_1 = [0.03, 0.02]$, $\gamma_2 = [0.02, 0.01]$, $\sigma_1 = [0.05, 0.1]$, $\sigma_2 = [0.15, 0.1]$, $\sigma_3 = [0.2, 0.1]$, $\sigma_4 = [0.1, 0.2]$, then $\Delta$ of this paper in environment 1 denoted by $\Delta_1$ equals to $0.6368$$>$$0$, $\Delta$ in environment 2 denoted by $\Delta_2$ is $-0.713 < 0$ and $\Delta$ in the whole environment is $-0.522 < 0$. By virtue of Theorem $3.1$, it has that the disease will last in environment 1 (see Figure 1(a)), disappear in environment 2 (see Figure 1(b)), and will also go extinct in the whole environment (see Figure 1(c)).

Example 2. The example here discusses the impact of isolation on disease control through numerical simulation. For simplicity, we study only the situation in one environment, that is, there is no Markovian switching. Take the parameters in environment 1 in Example 1, except the isolation rate $\delta$. We take two different values to compare the size of $S_t$ and $I_t$ in the model. Let $\delta_1 = 0.3$ and $\delta_2 = 0.6$ respectively, then $\Delta = 0.7868 > 0$ and $\Delta = 0.4868 > 0$, the disease will go on. We see from Figure 2(a) that the the size of disease $I_t$ with $\delta = 0.6$ is less than size of disease with $\delta = 0.3$. We know that with the increase of the isolation rate, more and more infected people are isolated (depending on the severity, they can be isolated at home), which will reduce the transmission to varying degrees. In addition, the severity of symptoms of $I_t$ and $Q_t$ people may be different, then different treatment measures for $I_t$ and $Q_t$ will save a certain amount of medical resources, which can make people recover and increase the size of susceptible class $S_t$, see Figure 2(b).

Example 3. This example will verify the effect of vertical transmission rate $p$ on disease behavior. Similar to the situation in Example 2, we only discuss one environment. Assume that $b = 0.08$, $\Lambda = 0.07$, $\mu = 0.1$, $\beta = 0.4$, $r_1 = 0.25$, $\alpha_1 = 0.45$, $\alpha_2 = 0.5$, $\delta = 0.4$, $\gamma_1 = 0.03$, $\gamma_2 = 0.02$, $\sigma_1 = 0.05$, $\sigma_2 = 0.15$, $\sigma_3 = 0.2$, $\sigma_4 = 0.1$, the initial data $S_0 = 2.1$, $I_0 = 0.6$, $Q_0 = 0.2$, $R_0 = 1.2$. We compare the size of different classes of $p$ under two values, let $p = 0.8$ and $p = 0.1$, then $\Delta$ with $p = 0.8$ equals to $0.4727 > 0$ and $\Delta = 0.4167$ under $p = 0.1$, the disease will last. From Figure 3(a), we see that a higher vertical transmission rate will produce more infected people. Under the same isolation rate, $Q_t$ will also become larger, see Figure 3(b). The increase of vertical transmission rate $p$ makes the individuals of the susceptible in population smaller, see Figure 3(c).

7.   Conclusions and future research

In this article, we study a class of a stochastic hybrid SIQRS model with nonlinear incidence and vertical transmission and gives a threshold $\Delta$ to distinguish different behaviors of the model. The disease will die out when $\Delta < 0$. If $\Delta > 0$, the model we discuss admits an invariant measure. In proving the latter conclusion, we construct a new class of Lyapunov functions. The values obtained in this paper are the same, while many other studies differ in the values of different behaviors.

Some other issues are worthy of concern. Some diseases do not have symptoms at the initial stage of infection but in the latent period. Therefore, stochastic models with a latent period or time delay can be studied. Models with other types of noise such as L$\rm{\acute{e}}$vy noise can be discussed. In practice, measures such as media coverage and vaccination will be taken to control diseases, so introducing these measures into the model and analyzing their impacts can be further investigated in the future. Moreover, the optimal control problems of measures that appear in the model can also be discussed. We leave these issues for further discussion.

Use of AI tools declaration

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

Acknowledgments

This research was supported by the Science and Technology projects of Jiangxi Province Education Department (No. GJJ212716) and the National Key R & D Program of China (No. 2023YFC3008902).

Conflict of interest

The authors declare that they have no competing interests.