The Unique ergodic stationary distribution of two stochastic SEIVS epidemic models with higher order perturbation

1.   Introduction

Recent years have witnessed a rapid development in the understanding of the transmission mechanism and the prevalence of diseases in mathematical epidemiology by compartmental modeling. Common categories of compartmental epidemic models, such as SIS ^[1,2], SIR ^[3,4,5], SEIR^[6,7,8,9] and some other types of epidemic models ^{[10,11,12,13]}, have been extensively studied by many researchers. Considering that using a vaccine to inhibit the spread of disease is an effective and sustainable method, researchers try to explore the effect of the vaccination measure by introducing a vaccinated ($V$) class into a mathematical model. Some clinical outcomes ^[14,15,16] have demonstrated that vaccines only provide temporary immunity to the disease, i.e., a person who has received the vaccination will once again be susceptible to the disease after the vaccine wears off. This contributes significantly to improving our understanding of disease prevention and control. Thus, it is necessary to investigate the SEIVS epidemic model (see ^[17,18]). Recently, Wang et al.^[19] put forward the following SEIVS epidemic model with latency and temporary immunity described by a system of ordinary differential equations:

where $0 < \rho\leq1$, $0\leq\kappa\leq2$, $\alpha > 0$, $V$ stands for the vaccinated compartment and the recovered and all parameters are positive. Parameter $A$ is the recruitment rate, $p$ stands for the fraction of the newborns vaccinated, $\mu$ represents the natural death rate, $\beta$ measures the disease transmission coefficient, $e$ refers to the inhibitory effect, $\eta$ depicts the recovery rate of $E$ due to natural immunity, $1/\delta$ means the average time of immunity waning, $1/\varepsilon$ represents the latent period and $1/\tau$ represents the mean infectious period. In particular, if $\alpha = \kappa = \rho = 1$ for the nonlinear incidence rate $\beta I^\rho S^\alpha/(1+eI^\kappa)$ in model (1.1), then the regressive epidemic model with a saturated incidence rate becomes

Similar to ^[19], the matrices ${\bf F}$ and ${\bf V}$ of model (1.2) respectively take the form ${\bf F} = (\texttt{f}_{ij})_{2\times 2}$ and ${\bf V} = (\texttt{v}_{ij})_{2\times 2}$ where $\texttt{f}_{11} = 0$, $\texttt{f}_{12} = \beta S_0$, $\texttt{f}_{21} = 0$, $\texttt{f}_{22} = 0$, $\texttt{v}_{11} = \mu+\varepsilon+\eta$, $\texttt{v}_{12} = 0$, $\texttt{v}_{21} = -\varepsilon$, $\texttt{v}_{22} = \mu+\tau$ and $Q_0 = (S_0, 0, 0, V_0) = (\frac{A((1-p)\mu+\delta)}{\mu(\mu+\delta)}, 0, 0, \frac{pA}{\mu+\delta})$, so the basic reproduction number $R_c = \rho({\bf F}{\bf V}^{-1}) = \frac{\varepsilon\beta A((1-p)\mu+\delta)}{\mu(\mu+\tau)(\mu+\varepsilon+\eta)(\mu+\delta)}$ of model (1.2), which determines whether an epidemic will develop. Moreover, the corresponding dynamical results are summarized as follows:

$\bullet$ If $R_c\leq1$, the disease-free equilibrium $Q_0 = (S_0, 0, 0, V_0)$ of model (1.2) is globally asymptotically stable, which means that the disease cannot spread.

$\bullet$ If $R_c > 1$, the endemic equilibrium $Q_\ast = (S_\ast, E_\ast, I_\ast, V_\ast)$ of model (1.2) is globally asymptotically stable, which means that the disease always remains in a population.

In the real world, environmental noises are ubiquitous and may more and less affect the transmission of epidemics. To reflect this fact, it is more realistic to research the corresponding stochastic model by incorporating environmental noises into the deterministic model (1.2). Particularly, white noise and telephone noise are two common types of noises. On the one hand, many scholars have done numerous significant studies by considering white noise and obtained rich results^{[20,21,22,23]}. For example, an insightful work of Mao et al.^[20] clearly indicated that white noise could effectively inhibit the explosion of a potential population. Note that Lu et al.^[21] investigated an SEIQV epidemic model by considering higher-order perturbation and proved the ergodic stationary distribution (ESD) and the extinction, and they also discussed the equilibrium stability of the model. Rajasekar et al.^[22] derived sufficient conditions for the ESD and the extinction of a second-order perturbed SIRS epidemic model with relapse and media impact. Along with the interesting idea of higher-order white noise, we introduce it into model (1.2) and formulate a stochastic version:

where the initial values and all parameters are positive, $B_i(t)$ $(i = 1, \cdots, 4)$ are independent standard Brownians and $\sigma^2_{i1}$ and $\sigma^2_{i2}$ represent their intensities.

On the other hand, the system may switch from an environmental regime to another when it is affected by telegraph noise^[24,25,26]. The switching is generally memoryless and the waiting time for the next switch follows exponential distribution. Mathematically, this switching of environmental regimes can be characterized by a continuous-time Markov chain $(r(t))_{t\geq0}$ with a finite state space $\mathbb{S} = \{1, 2, ..., N\}$ and the generator matrix $\Gamma = (\gamma_{ij})_{N\times N}$, for $\Delta > 0, \gamma_{ij}\geq 0$:

where $\gamma_{ij} > 0$ is the transition rate from state $i$ to state $j$ satisfying $\sum_{j = 1}^N\gamma_{ij} = 0$ when $i\neq j$ (see ^[27,28]). It is worthy of attention that a mass of epidemic models with both white noise and telegraph noise have been studied^{[29,30,31,32]}. For instance, Han and Zhao investigated an SIRS epidemic model influenced by two types of noises and obtained the asymptotic stability of the disease-free equilibrium point of the corresponding deterministic model ^[29]. Zhou et al. constructed an SEQIHR epidemic model with media coverage and Markov switching, and obtained the extinction and the ESD ^[31]. Motivated by these arguments, we further formulate a new stochastic version of model (1.2) with white noise and Markov switching, as follows:

As a continuation of the work previously studied by Wang et al. ^[19], we introduce white noise and telegraph noise into the regression model (1.2) with saturation incidence and temporary immunity and further propose two stochastic SEIVS models (models (1.3) and (1.4)) with high-order perturbation. Furthermore, the primary contributions of this study aim to analyze the stationary distributions of models (1.3) and (1.4), which are concerned with the stochastic statistical characteristic of the long-term behaviors of the sample trajectories. To the best of our knowledge, the existence and uniqueness of the stationary distributions of models (1.3) and (1.4) have not been investigated in any of the published articles so far. We now provide a brief outline of the paper. First, we define some crucial mathematical concepts and lemmas in Section 2. Section 3 focuses on the theoretical analysis results of models (1.3) and (1.4) by using Has'minskii theory and Lyapunov functionals. Numerical examples are given in Section 4. In the end, this paper's conclusions and future directions are provided in Section 5.

2.   Preliminaries

This section provides several useful preliminaries and lemmas that will be applied to the proof of the dynamic behaviors of models (1.3) and (1.4). Define $\mathbb{R}_+^n = \{z\in\mathbb{R}^n|z_i > 0, 1\leq i\leq n\}$. For the Markov chain $r(t)$, assume that it is irreducible; then, there is a unique stationary distribution $\pi = \{\pi_1, \pi_2, ..., \pi_N\}$ which is determined by $\pi\Gamma = 0$ and $\sum_{k = 1}^N\pi_k = 1$ $(\pi_k > 0)$ for any $k\in\mathbb{S}$. For the bounded constant sequence $h(k)$, assign $\hat{h} = \min_{k\in\mathbb{S}}\{h(k)\}$ and $\check h = \max_{k\in\mathbb{S}}\{h(k)\}$. The Markov chain $r(t)$ is independent of the Brownian motion $B_i(t)(i = 1, 2, 3, 4)$. The initial values and coefficients $p(k)$, $A(k)$, $\mu(k)$, $e(k)$, $\beta(k)$, $\delta(k)$, $\varepsilon(k)$, $\eta(k)$, $\tau(k)$, $\sigma_{i1}(k)$ and $\sigma_{i2}(k)$ of model (1.4) are positive for any $k\in\mathbb{S}$.

Lemma 2.1. ^[33] For any $x\geq0$, one has

Let $X(t)$ be a regular time-homogeneous Markov process in $\mathbb{R}^n$ described by the following stochastic differential equation

The diffusion matrix of the Markov process $X(t)$ is given by

Lemma 2.2. ^[34] The Markov process $X(t)$ admits a unique ESD $\pi(\cdot)$ if there exists an open domain $\zeta\subset\mathbb{R}^n$ with the regular boundary $\hbar$ such that

$(\mathcal{H}_1)$: for all $x\in\zeta$, the diffusion matrix $\Lambda(x)$ is strictly positive definite;

$(\mathcal{H}_2)$: for any $\mathbb{R}^n\setminus\zeta$, $LU$ is negative, where $U$ is a nonnegative $C^2$-function.

Suppose that the diffusion process $(X(t), r(t))$ satisfies the following equation

Assign $f_2(\cdot, \cdot):\mathbb{R}^n\times\mathbb{S}\rightarrow \mathbb{R}^n$ and $\mathcal{G}(\cdot, \cdot):\mathbb{R}^n\times\mathbb{S}\rightarrow \mathbb{R}^{n\times n}$ such that $\mathcal{G}(x, k)\mathcal{G}^T(x, k) = (d_{ij}(x, k))$. Define the operator $L$ related to (2.2) as below

where $H(x, k)$ is twice continuously differentiable with respect to $x$.

Lemma 2.3. ^[34] System (2.2) admits a unique ESD if the following assumptions hold:

${\rm{(i)}}$ $\gamma_{ij} > 0$ for any $i\neq j$;

${\rm{(ii)}}$ For each $k\in\mathbb{S}$, the matrix $\mathcal{P}(x, k) = (\mathcal{P}_{ij}(x, k))$ is symmetric satisfying

with some constant $\varphi\in(0, 1]$;

${\rm{(iii) }}$ There exists a twice continuously differentiable function $H(\cdot, k):\mathcal{D}^c\times\mathbb{S}\rightarrow \mathbb{R}$ such that

where $\mathcal{D}^c$ is the complement of a bounded open subset $\mathcal{D}\in\mathbb{R}^n$ with a smooth boundary. Moreover, there exists a unique stationary density $\pi(\cdot, \cdot)$ for any Borel measurable function $\varsigma(\cdot, \cdot):\mathbb{R}^n\times\mathbb{S}\rightarrow \mathbb{R}$ such that

which means that

In what follows, two lemmas are exhibited to prove the existence and uniqueness of global positive solutions for models (1.3) and (1.4), respectively.

Lemma 2.4. For any initial data $X(0) = (S(0), E(0), I(0), V(0))\in\mathbb{R}^4_+$ and $t\in[0, +\infty)$, model (1.3) admits a unique global positive solution $X(t) = (S(t), E(t), I(t), V(t))$ a.s.

Proof. Define a $C^2$-function $\mathcal{W}_1$: $\mathbb{R}^4_+\rightarrow \mathbb{R}_+$ for $0 < a < 1$:

Utilizing Itô's formula, we get

where

Define

we can obtain

and

Let $f''(S) = 0$, we get

When $S < S_0$, we have $f''(S) > 0$; then, $f'(S)$ is monotonically increasing. While $S > S_0$, we have $f''(S) < 0$; then, $f'(S)$ is monotonically decreasing. Obviously, $f'(0) = \frac{\beta C_2}{e}+\sigma_{11}\sigma_{12} > 0,$ so we can obtain that $f'(S)$ has a maximum value $f'(S_0) > 0$.

Also, when $S\rightarrow +\infty$, $f'(S)\rightarrow -\infty$ which means that there must exist $S_1$ such that $f'(S_1) = 0$. Similarly, $f'(S) > 0$ when $S\in(0, S_1)$; then, $f(S)$ is monotonically increasing and $f'(S) < 0$ when $S\in(S_1, +\infty)$; then, $f(S)$ is monotonically decreasing. Since $f(S)$ is a continuous function, $f(S)$ has a guaranteed upper-bound which is positive. The same analysis can be applied to $E$, $I$ and $V$. Then we obtain

Thus $P_1$ is a positive constant. Hence we can obtain $L\mathcal{W}_1\leq P_1$. The remaining proof is analogous to that of Theorem 2.1 in ^[35]; hence, we skip the details.

Lemma 2.5. For any initial data $(X(0), r(0))\in\mathbb{R}^4_+\times\mathbb{S}$, model (1.4) has a unique positive solution $(X(t), r(t))$ on $t\geq0$ a.s.

Proof. Define $\mathcal{W}_2(S, E, I, V, k):\mathbb{R}_+^4\times\mathbb{S}\rightarrow \mathbb{R}$:

where $\varrho\in(0, 1)$ is a constant and $k\in\mathbb{S}$. Note that

where $\mathbb{G}_n = (\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)$. It is easy to see that $\mathcal{W}_2$ has a minimal value $\overline{\mathcal{W}}_0$ in $\mathbb{R}^4_+\times\mathbb{S}$; then, we consider a function $\overline{\mathcal{W}}_2(S, E, I, V, k)$, as follows:

By using Itô's formula, we have

where

Similar to the proof of the upper-bound $P_1$ of $y_1(S, E, I, V)$ (see (2.5)–(2.9)), we can obtain a positive constant $P_2$ such that

As a result, we can get $L\overline{\mathcal{W}}_2\leq P_2.$ The remaining proof is nearly identical to those in Theorem 2.1 of ^[36] and is omitted.

3.   Main results

In this section, by using the Has'minskii theory ^[34] and a Lyapunov functional, we shall prove the existence of the unique ESD for models (1.3) and (1.4), respectively.

Theorem 3.1. If

then for any initial data $X(0)\in\mathbb{R}_+^4$, model (1.3) admits a unique ESD.

Proof. We only need to verify that each assumption in Lemma 2.2 holds when Theorem 3.1 is valid. In what follows, we divide this proof into the following two steps.

Step 1 (Positive definiteness of diffusion matrix). According to model (1.3), we have

Obviously, $\Lambda(S, E, I, V)$ is positive definite. Therefore, the assumption $(\mathcal{H}_1)$ in Lemma 2.2 holds.

Step 2 (Construction of a non negative $C^2$-function). It is worth mentioning that the crucial part of this step is to construct a suitable non negative Lyapunov function so that the assumption $(\mathcal{H}_2)$ holds. For convenience, we first define the following $q$-stochastic critical value $R_0^s(q)$ corresponding to $R_0^s$:

where $q\in(0, 1)$ is a sufficiently small constant. Thus, it is easy to derive that

When $R_0^s > 1$, it follows from the continuity of $R_0^s(q)$ that there exists an adequately small $q$ such that $R_0^s(q) > 1$. To do so, we define

where

and $D_1$, $a_i$ and $b_i$ $(i = 1, 2)$ will be given later. According to model (1.3), we calculate

and

Let us choose

and then (3.3) can be simplified as

combining this with (3.1) and (3.2) we get

Similarly, let

where

$D_2$ and $d_i$ $(i = 1, 2, 3)$ will be given later. Applying Itô's formula, one has

and

Assign

then, (3.8) can be written as

which, together with (3.6) and (3.7), leads to

Assign

where

$D_3$ and $\theta_1$ will be given later. Applying Itô's formula, one can obtain that

Combining (3.12) and (3.13), one has

Let $\theta_1 = 1/(1-q)\sigma_{31}^q$; then,

Assign

Adding (3.5), (3.11) and (3.14), one can get

where

Define

utilizing Itô's formula, we get

Let

where the positive number $M$ is large enough to satisfy $-M(R_0^s(q)-1)+C_5\leq-2$ and $C_5$ will be established later. Since $\widetilde{U}(S, E, I, V)$ is a continuous function, it follows that

where $W_\vartheta = (\frac{1}{\vartheta}, \vartheta)\times(\frac{1}{\vartheta}, \vartheta)\times(\frac{1}{\vartheta}, \vartheta)\times(\frac{1}{\vartheta}, \vartheta)$ and $\vartheta > 1$ is a sufficiently large integer.

Now, we consider the following nonnegative $C^2$-function:

where $\widetilde{U}_0$ is the minimal value of $\widetilde{U}$. Consequently, we have

Assign a suitable compact subset as follows:

where $0 < \epsilon < 1$ is a suitably small constant to satisfy the following inequalities:

and the constants $J$ and $C_5$ are given explicitly in (3.30) and (3.31), respectively. We split $\mathbb{R}_+^4\setminus W_\epsilon$ into the following eight regions:

Obviously, $\mathbb{R}_+^4\setminus W_\epsilon = W_\epsilon^c = W_1\cup W_2\cup W_3\cup W_4\cup W_5\cup W_6\cup W_7\cup W_8$. Then, we need to verify that $LU(S, E, I, V)\leq-1$ on $(S, E, I, V)\in\mathbb{R}_+^4\setminus W_\epsilon$.

Case 1. If $(S, E, I, V)\in W_1$, then we obtain

where

Case 2. If $(S, E, I, V)\in W_2$, then combining (3.21) with (3.23), one gets that

Case 3. If $(S, E, I, V)\in W_3$, then by (3.21) and (3.24), we derive

where

Case 4. If $(S, E, I, V)\in W_4$, according to (3.21) and (3.25), we obtain

Case 5. If $(S, E, I, V)\in W_5$, in accordance with (3.21) and (3.26), one can see that

Case 6. If $(S, E, I, V)\in W_6$, it follows from (3.21) and (3.27) that

Case 7. If $(S, E, I, V)\in W_7$, by (3.21) and (3.28), one has

Case 8. If $(S, E, I, V)\in W_8$, in view of (3.21) and (3.29), we have

In short, there is a sufficiently small $\epsilon$ such that

Then we draw a conclusion from Steps 1 and 2 that model (1.3) has a unique ESD.

Theorem 3.2. The solution $(X(t), r(t))\in\mathbb{R}^4_+\times\mathbb{S}$ of model (1.4) admits a unique ESD when

Proof. We employ Lemma 2.3 to verify Theorem 3.2, which needs to satisfy the assumptions (ⅰ)–(ⅲ). Here, we shall divide the whole proof of Theorem 3.2 into three steps.

Step 1 (Transition coefficient). According to the basic property of the generator $\Gamma$ in Section 2, namely $\gamma_{ij} > 0$, assumption (ⅰ) in Lemma 2.3 is clearly valid.

Step 2 (Positive definiteness of diffusion matrix). The diffusion matrix of model (1.4) is

Evidently, the matrix $\Phi(S, E, I, V, k)$ is positive definite and assumption (ⅱ) is verified.

Step 3 (Construction of a non negative $C^2$-function). In this step, we will show that assumption (ⅲ) holds. Since the proof is similar to Step 2 in the proof of Theorem 3.1, we pay attention to the different parts. For an adequately small constant $\psi\in(0, 1)$, define

Evidently, $\liminf\limits_{\psi\to0^+}R_1^s(\psi) = R_1^s$.

What is more, we construct a suitable nonnegative $C^2$-function:

where $\widetilde{H} = G(H_1+H_2+H_3)-\ln S-\ln I-\ln V+H_4$ and $\widetilde{H}_0$ is the minimal value of $\widetilde{H}$. Among them,

where $a_3$, $a_4$, $b_3$, $b_4$, $d_4$, $d_5$, $d_6$, $\theta_2$, $T_1(k)$, $T_2(k)$, $T_3(k)$ and $G$ will be given later, and

Moreover, we assume $G > 0$ is a large enough constant that satisfies $-G(R_1^s(\psi))+\tilde{C}\leq-2$ and

Then, applying Itô's formula to $H_1$, a calculation, similar to (3.5), is

where

and

Since $\Gamma$ is irreducible, for $\overrightarrow{M_1} = (M_1(1), M_1(2), \cdots, M_1(N))^\top$, there is a vector $\overrightarrow{T_1} = (T_1(1), T_1(2)$, $\cdots, T_1(N))^\top$ satisfying the following Poisson system $\Gamma\overrightarrow{T_1} = \sum_{l = 1}^N\pi_lM_1(k)-\overrightarrow{M_1}$, that is,

which is substituted into (3.37); we have

Similarly, applying Itô's formula to $H_2$, one obtains

where

and

For $\overrightarrow{M_2} = (M_2(1), M_2(2), \cdots,$ $M_2(N))^\top$, we determine a vector $\overrightarrow{T_2} = (T_2(1), T_2(2),$ $\cdots, T_2(N))^\top$ satisfying the Poisson system $\Gamma\overrightarrow{T_2} = \sum_{l = 1}^N\pi_lM_2(k)-\overrightarrow{M_2}$, which implies

This yields that

By applying Itô's formula to $H_3$, one obtains

where

Define $\overrightarrow{T_3} = (T_3(1), T_3(2), \cdots, T_3(N))^\top$; then, $\overrightarrow{M_3} = (M_3(1), M_3(2), \cdots,$ $M_3(N))^\top$ satisfies the Poisson system $\Gamma\overrightarrow{T_3} = \sum_{l = 1}^N\pi_lM_3(k)-\overrightarrow{M_3}$, which indicates

We conclude that

The remaining proof is similar to that of Theorem 3.1. Thus, we omit it here.

To sum up, the above three steps assure that Model (1.4) has a unique ESD under $R_1^s > 1$.

4.   Numerical simulations

In the present section, we are interested in providing two numerical examples to confirm the obtained mathematical results for models (1.3) and (1.4). For the sake of analysis, we fix the initial condition $(S(0), E(0), I(0), V(0)) = (0.8, 0.2, 0.3, 0.3)$, and the parameter values involved in both models are given according to different actual needs. The details are below.

Example 4.1 Let us consider model (1.3) with the parameters $\mu = 1/77$, $\delta = 0.05$, $e = 10^{-6}$, $p = 0.2$, $\beta = 0.35$, $A = 0.5$, $\tau = 0.4$, $\varepsilon = 0.7$, $\eta = 0.2$ and the corresponding noise values $\sigma_{11} = 0.01$, $\sigma_{12} = 0.02$, $\sigma_{21} = 0.01$, $\sigma_{22} = 0.02$, $\sigma_{31} = 0.005$, $\sigma_{32} = 0.005$, $\sigma_{41} = 0.005$ and $\sigma_{42} = 0.005$. By a calculation, one could calculate that $R^S_0 = 2.1439 > 1$. It is obvious from Theorem 3.1 that model (1.3) has a unique ESD. We plot the sample trajectories and the probability density functions of $S(t)$, $E(t)$, $I(t)$ and $V(t)$, respectively in Figure 2. The obvious conclusion is that reducing the noise intensities $\sigma_{11}^2$, $\sigma_{12}^2$, $\sigma_{21}^2$, $\sigma_{22}^2$ and $\sigma_{31}^2$ and decreasing the fraction of the newborns vaccinated $p$ or increasing the disease transmission coefficient $\beta$ will lead to the long-term persistence of disease.

Example 4.2 For model (1.4), we only consider the Markov chain $r(t)$ switch among these two states $\mathbb{S} = \{1, 2\}$ with the generator $\Gamma = (\gamma_{ij})_{2\times 2}$, where $\gamma_{11} = -\frac{5}{9}$, $\gamma_{12} = \frac{5}{9}$, $\gamma_{21} = -\frac{4}{9}$ and $\gamma_{22} = \frac{4}{9}$, and the unique stationary of $r(t)$ is given by $\pi = (\pi_1, \pi_2) = (\frac{4}{9}, \frac{5}{9})$. Next, we consider two sets of different parameter values to characterize different environmental regimes.

When $r(t) = 1$, we take $\mu(1) = 1/77$, $\delta(1) = 0.05$, $e(1) = 10^{-6}$, $p(1) = 0.2$, $\beta(1) = 0.35$, $A(1) = 0.5$, $\tau(1) = 0.4$, $\varepsilon(1) = 0.7$, $\eta(1) = 0.2$, $\sigma_{11}(1) = 0.01$, $\sigma_{12}(1) = 0.02$, $\sigma_{21}(1) = 0.01$, $\sigma_{22}(1) = 0.02$, $\sigma_{31}(1) = 0.01$, $\sigma_{32}(1) = 0.02$, $\sigma_{41}(1) = 0.01$ and $\sigma_{42}(1) = 0.02$.

When $r(t) = 2$, choose $\mu(2) = 0.28$, $\delta(2) = 2.9$, $e(2) = 0.025$, $p(2) = 0.12$, $\beta(2) = 0.40$, $A(2) = 0.4$, $\tau(2) = 0.65$, $\varepsilon(2) = 0.02$, $\eta(2) = 0.02$, $\sigma_{11}(2) = 0.05$, $\sigma_{12}(2) = 0.04$, $\sigma_{21}(2) = 0.05$, $\sigma_{22}(2) = 0.04$, $\sigma_{31}(2) = 0.03$, $\sigma_{32}(2) = 0.04$, $\sigma_{41}(2) = 0.03$ and $\sigma_{42}(2) = 0.04$.

Combining the above two sets of parameter values, it is easy to check that the condition of Theorem 3.2 is satisfied and $R_1^s = 2.7782 > 1$, which means that model (1.4) owns a unique ESD; see Figure 3. In addition, all of the probability density functions of $r(t)$, $S(t)$, $E(t)$, $I(t)$ and $V(t)$ have two distribution curves, which are caused by two different switching regimes. After a further observation, if we reduce the noise intensities $\sigma_{11}^2(k)$, $\sigma_{12}^2(k)$, $\sigma_{21}^2(k)$, $\sigma_{22}^2(k)$ and $\sigma_{31}^2(k)$ and the fraction of the newborns vaccinated $p(k)$, or if we increase the disease transmission coefficient $\beta(k)$, the disease will be prevalent and persistent in the long term.

5.   Discussion

In this work, we constructed two novel SEIVS models with latency and temporary immunity under higher-order perturbation, incorporating both white noise and telephone noise; see models (1.3) and (1.4). And, the existence and uniqueness of the stationary distribution have been obtained by utilizing Has'minskii theory and a Lyapunov function approach. We found that the noise can influence the dynamic behaviors of disease, which reveals that reducing the noise intensity leads to the persistence of disease.

Meanwhile, sufficient criteria on the existence of the unique ESD of the above two models are established. When $R^s_0 > 1$, it follows from Theorem 3.1 that model (1.3) owns a unique ESD. Reviewing Theorem 3.2, model (1.4) owns a unique ESD when $R^s_1 > 1$. It can be concluded from the conditions of these two theorems that the noise intensity, disease transmission coefficient and vaccination rate play a crucial role in the prevalence of disease.

Looking to the future, there are some intriguing questions that warrant further consideration and investigation. For example, one can develop more elaborate models, such as focusing on model (1.1) and taking into account environmental noises, the dynamic properties of the stochastic SEIVS model with the nonlinear incidence rate ($\frac{\beta I^\rho S^\alpha}{1+eI^\kappa}$) are worth exploring. In addition, there exist time delays in the process of information dissemination, taking into account the fact that the impact of time delays can be introduced into model (1.2) ^[8,12]. Furthermore, there are also some interesting modeling ideas and approaches from different perspectives, such as the multi-group epidemic dynamic models ^[10,11] and fractional-order models ^[13]. These works will be shown in other articles.

Acknowledgments

The authors would like to thank the editor and referees for their valuable comments and suggestions, which have greatly improved the presentation of this paper. The work was supported by the NNSF of China (No. 11871201).

Conflict of interest

The authors declare that there is no conflict of interest.