Stationary distribution and extinction of a stochastic SEIQ epidemic model with a general incidence function and temporary immunity

1.   Introduction

Infectious diseases, such as AIDS, Ebola, Zika, influenza, Hepatitis B virus, COVID-19 and so on, are all public health concerns of the world. It has an important impact on public health, social and economic development and even human life^[5,6,9,10]. By the reference of the World Health Organization (WHO) 1996: Nearly 50000 children, men and women die from infectious diseases every day^[41,42]. For many diseases, there are no drugs or vaccines for treatment or cure. We are on the verge of a global epidemic crisis. No country can do without them. No country can ignore their threat any more, the Director-General of WHO says in the WHO report^[41]. Thus, in order to better understand the transmission mechanism of infectious diseases, and obtain good methods to control the infectious diseases propagation, mathematicians have proposed many mathematical models to discuss the dynamic behavior of infectious disease transmission^{[2,16,17,22,35]}. For instance, the SIS (Susceptible-Infectious-Susceptible) epidemic model^[12,14,27], the SIR (Susceptible-Infected-Recovered) epidemic model^[3,20,21], the SIRI (Susceptible-Infected-Recovered-Infected) epidemic model^[33,37,44], the SIRS (Susceptible-Infected-Recovered-Susceptible) epidemic model^[25,32], SEIR (Susceptible-Expose-Infected-Recovered) epidemic model^[30,39,43] and the infectious diseases model with vaccination^[4,28].

In diseases such as influenza and COVID-19, some people can be isolated once they are found to have been infected with infectious disease in the exposed or infectious state. To portray the character of the quarantined state and combine with the characteristics of classical epidemic models, Chen et al.^[8] proposed the infectious disease model with quarantined compartment as follows:

where $S(t)$, $E(t)$, $I(t)$ and $Q(t)$ represent the number of the susceptible, exposed, infected and quarantined individuals at any time $t$, respectively. The meanings of the parameters are shown in Table 1. In ^[8], many important results, such as the local and global asymptotic stability of the disease-free equilibrium and the endemic equilibrium of model (1.1), were obtained under certain conditions, and estimated the domain of attraction for the endemic equilibrium.

As is known to all, the incidence rate of infectious diseases plays a key role in the investigation of mathematical epidemiology. Most classical infectious disease models employ bilinear incidence rates, which is based on the homogeneous mixing assumption^[2,29,45]. However, experiments ^[11] have strongly suggested that the incidence is an increasing function of the parasite dose, and is usually sigmoidal in shape (see, for example, ^[36]). As a result, Anderson and May ^[1] proposed a saturated incidence $\frac{\alpha SI}{1+\beta S}$; Lv et al.^[31] introduced a Beddington-DeAngelis incidence rate $\frac{\beta xv}{1+mx+nv}$ into the HIV model; Caraballo et al.^[7] considered an incidence function $(\beta_{1}-\frac{\beta_{2}I}{m+I})SI$ to discuss the influence of media coverage on propagation dynamics. Motivated by the above discussion, we present the following SEIQ epidemic model with general incidence function:

where the meanings of all parameters in model (1.2) are consistent with that in model (1.1). By a simple calculation, we obtain a basic reproduction number

Similar to the proof method in ^[8], we also show that if $\mathcal{R}_{0} < 1$, then system (1.2) has a unique free disease-equilibrium, which is globally asymptotically stable. The disease equilibrium is unique and also globally asymptotically stable when $\mathcal{R}_{0} > 1$.

On the other hand, due to the existence of environmental noise, the parameters involved in system (1.2) are not constants, and they always fluctuate near the average value. Therefore, numerous researchers considered the influence of random factors in the epidemic models and studied the dynamic behavior of stochastic epidemic models with general incidence rate^[19,25,26]. Fan et al.^[15] presented a stochastic delayed SIR epidemic model with generalized incidence rate and temporary immunity. On the basis of the strong law of large numbers, they established sufficient conditions for extinction and permanence in the mean. Fatini et al.^[13] adopted parameter perturbation method to investigate the stationary distribution of a generalized SIRS epidemic model. They considered the parameter $\beta$ fluctuates near the average value, and showed the extinction of the disease when $\mathcal{R}_{0}^{s} < 1$ while the persistence in mean of the disease when $\mathcal{R}_{0}^{s} > 1$. Liu et al.^[24] considered a stochastic SIRS epidemic model with logistic growth and general nonlinear incidence rate. By constructing a suitable stochastic Lyapunov function, they obtain sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of the positive solutions to the stochastic SIRS epidemic model.

Based on the above concerns, in this paper, we try to do some work in this field, and discuss the existence of a stationary distribution of solutions of stochastic SEIQ epidemic model. Our method is similar to that of Wang et al.^[40]. Now, we assume the environmental fluctuations are white noise type, which are directly proportional to $S$, $E$, $I$ and $Q$, and influenced on $\frac{dS}{dt}$, $\frac{dE}{dt}$, $\frac{dI}{dt}$ and $\frac{dQ}{dt}$ in model (1.1). Then we obtain the following stochastic SEIQ epidemic model corresponding to model (1.2):

where $B_{i}(t)$ are mutually independent standard one-dimensional motion of Brownian motions with $B_{i}(0) = 0$, $\sigma_{i}^{2} > 0$ mean the intensities of the white noise, $i = 1, 2, 3, 4$. Furthermore, we assume that the contact between a susceptible individual and a infected individual is given by an incidence function $f(\cdot)$, which satisfies the following conditions:

${\rm(A1)}$ $f(\cdot)\in C^1(\mathbb{R}_{+}, \mathbb{R}_{+})$ such that $f'(I) > 0$ and $f(0) = 0$ $\forall$ $I\geq0$;

${\rm(A2)}$ $\frac{f(I)}{I}$ is monotonically non-increasing on $(0, \infty)$ and $f(I) > If'(I)$ for all $u\in \mathbb{R}_+$.

From (A2), one can see that $f(I) < If'(0)$ for all $I > 0$. These assumptions are biologically motivated, see ^[13] in detail. Obviously, it includes the common incidence functions such as the linear function $f(I) = \beta I$ and the saturation incidence function $f(I) = \frac{\beta I}{1+aI}$, where $\beta$, $a > 0$.

As is known to all, in the deterministic epidemic model, the basic reproduction number and endemic equilibrium are two very important factors. Because the basic reproduction number represents the threshold of the epidemic spreading, and the endemic equilibrium represents the prevalence of the disease. However, the influence of environmental noises is taken into account, it is difficult to determine the threshold of stochastic epidemic model, and most stochastic models have no the endemic equilibrium. Therefore, the stationary distribution of stochastic models has been widely concerned. On the other hand, from the perspective of biomathematics, the existence and ergodicity of stationary distribution illustrates that infectious diseases will exist for a long time and continue to develop. Therefore, we concentrate on two points: (i) develop a stochastic SEIQ epidemic model with generalized incidence function $Sf(I)$, and on the basis of the properties of $f(I)$, prove the existence and uniqueness of the global positive solution. (ii) try to give the extinction threshold of the disease and investigate the stationary distribution of stochastic epidemic model with generalized incidence function.

The paper is organized as follows. In Section 2, we introduce some basic definitions and related theories in this paper. In Section 3, we show the existence and uniqueness of the global positive solution of system (1.3) with any initial value. In Section 4, we give some sufficient conditions for extinction of the disease. In Section 5, we show the existence of a stationary distribution of solutions to model (1.3) by using the method of Khasminskii and constructing a suitable Lyapunov function. In the Section 6, some numerical simulations are presented to verify the theoretical results. The main results of this paper ends with conclusion in Section 7.

2.   Preliminaries

Throughout the paper, let $(\Omega, \mathscr{F}, \mathbb{P})$ be a complete probability space with filtration $\{\mathscr{F}_{t}\}_{t\geq0}$ satisfying the usual conditions (i.e., it is right continuous and increasing while $\mathscr{F}_{0}$ contains all $\mathbb{P}$-avoid set) and $B_{i}(t)$ are denoted on the complete probability space. Denote $\mathbb{R}_{+}^{4} = \{x\in\mathbb{R}_{+}^{4}: x_{i} > 0, i = 1, 2, 3, 4\}$, $a\vee b = \max\{a, b\}$ and $a\wedge b = \min\{a, b\}$, for all $a, b\in\mathbb{R}$.

For the four-dimensional stochastic differential equation

with the initial value $x(0) = x_{0}\in\mathbb{R}^{4}$. Let $C^{2, 1}(\mathbb{R}^{4}\times(0, \infty); \mathbb{R}_{+})$ denote the family of all nonnegative functions $V(x, t)$ which defined on $\mathbb{R}^{4}\times(0, \infty)$ such that they are continuously twice differentiable in $x$ and once in $t$. Define the differential operator $L$ of (2.1) as follows^[34]

Acting $L$ on $V\in C^{2, 1}(\mathbb{R}^{4}\times(0, \infty); \mathbb{R}_{+})$, we have

where $V_{t} = \frac{\partial V}{\partial t}$, $V_{x} = (\frac{\partial V}{\partial x_{1}}, \frac{\partial V}{\partial x_{2}}, \frac{\partial V}{\partial x_{3}}, \frac{\partial V}{\partial x_{4}})$, $V_{xx} = \left(\frac{\partial^{2}V}{\partial x_{i}\partial x_{j}}\right)_{4\times4}$. Thus, by the Itô's formula, if $x(t)\in\mathbb{R}^{4}$, then

To obtain the stationary distribution of system (1.3), we propose some definitions and known results in the sequel. Let $X(t)$ be homogeneous Markov process in the $d$-dimensional Euclidean space $\mathbb{R}^{d}$, which satisfies the autonomous equation

It is well-known that the diffusion matrix corresponding to Eq (2.2) has the form

Lemma 2.1. ^[23] Let $U\in\mathbb{R}^{d}$ be a bounded open domain with regular boundary $\Gamma$, which has the following properties:

${\rm (\textbf{H1})}$ There exists a constant $M > 0$; $\sum\limits_{i, j = 1}^{l}a_{ij}(x)\xi_{i}\xi_{j}\geq M|\xi|^{2}$, for all $x\in U$ and $\xi\in\mathbb{R}^{d}$;

${\rm (\textbf{H2})}$ There is a $C^{2}-$function $V\geq0$ satisfying $LV < 0$ for any $x\in\mathbb{R}^{d}-U$.

Then the Markov process $X(t)$ admits a unique ergodic stationary distribution $\pi(\cdot)$. In this situation,

where $f(\cdot)$ is an integrable function with respect to the measure $\pi$.

3.   Existence and uniqueness of the global positive solution

In order to investigate dynamics of an SEIQ model, the first concern is whether the solution is global and positive. In this section, we will construct a suitable Lyapunov function to verify the existence and uniqueness of global positive solutions of the model (1.3).

Theorem 3.1. For any initial value $(S(0), E(0), I(0), Q(0))\in\mathbb{R}_{+}^{4}$, a unique positive solution $(S(t), E(t), I(t), Q(t))$ of model (1.3) exists on $t\geq0$ and the solution remains in $\in\mathbb{R}_{+}^{4}$ with probability one, namely, the solution $(S(t), E(t), I(t), Q(t))\in\mathbb{R}_{+}^{4}$ for all $t\geq0$ almost surely.(a.s.)

Proof. Since the coefficients of model (1.3) satisfy the local Lipschitz condition, for the initial value $(S(0), E(0), I(0), Q(0)) \in\mathbb{R}_{+}^{4}$, model (1.3) exists a unique local solution $(S(t), E(t), I(t), Q(t))$ on $t\in[0, \tau_{e}]$, where $\tau_{e}$ is the explosion time^[34]. To prove the solution is global, we need to show $\tau_{e} = \infty$ a.s.. Let $n_{0}$ be large enough such that $S(0)$, $E(0)$, $I(0)$ and $Q(0)$ all lie within the interval $[\frac{1}{n_{0}}, n_{0}]$. For every $n\geq n_{0}$, define the following sequence of stopping time

where set $\inf\emptyset = \infty$(in general, $\emptyset$ is as the empty set). Obviously, $\tau_{n}$ is increasing as $n\rightarrow \infty$. Let $\tau_{\infty} = \lim\limits_{n\rightarrow \infty}\tau_{n}$, then $\tau_{\infty}\leq\tau_{n}$ a.s.. If we can verify $\tau_{\infty} = \infty$ a.s., then $\tau_{e} = \infty$ a.s. and $(S(t), E(t), I(t), Q(t))\in\mathbb{R}_{+}^{4}$ a.s. for all $t\geq0$. Namely, to complete the proof, we only need to prove $\tau_{\infty} = \infty$ a.s.. If the statement is violated, then there exists a constant $T > 0$ and $\epsilon\in(0, 1)$ such that

Define a non-negative $C^{2}-$function $V$: $\mathbb{R}_{+}^{4}\rightarrow \mathbb{R}_{+}$

where $a$ is a positive constant to be determined later. Applying Itô's formula^[34], we get

Where

Choose $a = \frac{\mu+\alpha+\gamma_{2}}{f'(0)}$ such that $af'(0)- (\mu+\alpha+\gamma_{2}) = 0$, then we obtain

where $K$ is a suitable constant. Consequently,

For any $n\geq n_{0}$, integrating the both sides of (3.5) from 0 to $\tau_{n}\wedge T$, and taking the expectation, we have

Moreover, we set $\Omega_{n} = \{\tau_{n}\leq T\}$ for any $n\geq n_{0}$, then (3.2) leads to $\mathbb{P}(\Omega_{n})\geq\epsilon$. For every $\omega\in\Omega_{n}$, at least one of $S(\tau_{n}, \omega)$, $E(\tau_{n}, \omega)$, $I(\tau_{n}, \omega)$ or $Q(\tau_{n}, \omega)$ equals to $n$ or $\frac{1}{n}$. Thus, $V(S(\tau_{n}, \omega), E(\tau_{n}, \omega), I(\tau_{n}, \omega), Q(\tau_{n}, \omega))$ is no less than either

It follows that

Therefore,

where $1_{\Omega_{n}}$ means the indicator function of $\Omega_{n}$. Let $n\rightarrow \infty$, we have

which leads to a contradiction. Therefore, the claim $\tau_{\infty} = \infty$ a.s. hold. This completes the proof.

Remark 3.1. Theorem 3.1 illustrates that for any positive initial value, model (1.3) has unique positive solution with probability one. Here, let $N(t) = S(t)+E(t)+I(t)+Q(t)$, then

Consider the following system

By Theorem 3.9 in ^[34], we have $\lim\limits_{t\rightarrow \infty}X(t) < \infty$ a.s. Then according to the comparison theorem for stochastic differential equations^[38], this yields $\limsup\limits_{t\rightarrow \infty}N(t) < \infty$ a.s.

4.   Extinction of the diseases

One of the main concerns of epidemiology is how we regulate the dynamics of the disease in order to eliminate it in a long term. In this section, we give some sufficient conditions for the extinction of the disease in model (1.3). For the convenience, we define

Lemma 4.1. ^[46] Let $(S(t), E(t), I(t), Q(t))$ be the solution of model (1.3) with initial value $(S(0), E(0), I(0), Q(0))\in\mathbb{R}_{+}^{4}$. Then

Furthermore, if $\mu > \frac{\sigma_{1}^{2}\vee\sigma_{2}^{2}\vee\sigma_{3}^{2}\vee\sigma_{4}^{2}}{2}$, then

Define the parameter as follows:

Theorem 4.1. Suppose $(S(t), E(t), I(t), Q(t))$ is the solution of model (1.3) with initial value $(S(0), E(0), I(0), Q(0))\in\mathbb{R}_{+}^{4}$. If $\widehat{\mathcal{R}}_{0}^{*} < 1$ and $\mu > \frac{\sigma_{1}^{2}\vee\sigma_{2}^{2}\vee\sigma_{3}^{2}\vee\sigma_{4}^{2}}{2}$, then the solution $(S(t), E(t), I(t), Q(t))$ of model (1.3) has

and

Proof. For model (1.3), we have

Integrating the both sides of (4.3) from 0 to $t$ and combining with Lemma 2, we obtain

Let $P(t) = \varepsilon E(t)+(\mu+\varepsilon+\delta_{1}+\gamma_{1})I(t)$, by the Ito's formula^[34], we have

where

Obviously,

Therefore, we have

Integrating the both sides of (4.5) from 0 to $t$, together with Lemma 2 and $\widehat{\mathcal{R}}_{0}^{*} < 1$, we obtain

which implies that

Therefore, from the last equation of model (1.3), we can easily obtain $\lim\limits_{t\rightarrow \infty}Q(t) = 0$ a.s.

On the other hand, for any small $\xi > 0$, there exists $T > 0$ such that $I\leq\frac{\xi}{f'(0)}$ for $t\geq T$ a.s. Substituting $I\leq\frac{\xi}{f'(0)}$ into the first equation of model (1.3), we obtain

Integrating the both sides of (4.6) from 0 to $t$ and together with Lemma 2, we get

The arbitrariness of $\xi$, one has

Thus, (4.4), (4.6) and (4.8) imply that

This completes the proof.

Remark 4.1. Theorem 4.1 illustrates that if $\widehat{\mathcal{R}}_{0}^{*} < 1$, the disease will go to extinct. Note the express of $\widehat{\mathcal{R}}_{0}^{*}$, we notice that the larger the intensity of white noises are, the easier the extinction of the disease is. Thus, we can control the outbreak of disease by adjusting the intensity of environmental noises.

5.   Stationary distribution

In studying epidemiological dynamics, we are not only concerned about when the disease will be extinct in the population, but also about when the disease will persist in the population. For model (1.3), there is no endemic equilibrium. Therefore, in this section, we use the method of Khasminskii^[23] to focus on the existence of a stationary distribution, which reflects whether the disease is prevalent.

Theorem 5.1. Suppose

then for any initial value $(E(0), S(0), Q(0), I(0))\in\mathbb{R}_{+}^{4}$, model (1.3) has a ergodic stationary distribution $\pi(\cdot)$.

Proof. From (2.3), we know the diffusion matrix of model (1.3) has the form

According to Lemma 1, we only need to verify the conditions (H1) and (H2). Let $M: = \min\limits_{(S, E, I, Q)\in\overline{U}_{n}\subset\mathbb{R}_{+}^{4}}\big\{\sigma_{1}^{2}S^{2}, \sigma_{2}^{2}E^{2}, \sigma_{3}^{2}I^{2}, \sigma_{4}^{2}Q^{2}\big\} > 0$, then we have

for all $(S, E, I, Q)\in\overline{U}_{n}: = [{1}/{n}, n]^{4}$ provided that $n\in\mathbb{Z}_+$ is sufficiently large. Note that $\xi = (\xi_{1}, \xi_{2}, \xi_{3}, \xi_{4})\in\mathbb{R}^{4}$. Hence, the condition (H1) holds.

To verify (H2), we denote

Since $\mu+\frac{\sigma_{1}^{2}}{2} > 0$ and $\hat{\mathcal{R}}_{0} > 1$, we have $b > 0$. Now, define the $C^{2}$-function $V(\cdot): \mathbb{R}_{+}^{4}\rightarrow \mathbb{R}_{+}$ by

with

where

$p$ and $\rho$ are positive constants satisfying the conditions:

and

Therefore, we can easily check that

where $\mathbb{U}_{k} = (\frac{1}{k}, k)\times(\frac{1}{k}, k)\times(\frac{1}{k}, k)\times(\frac{1}{k}, k)$. Hence, $V(S, E, I, Q)$ is a continuous function and has a minimum point($\overline{E}_{0}$, $\overline{S}_{0}$, $\overline{Q}_{0}$, $\overline{I}_{0}$) in the interior of $\mathbb{R}_{+}^{4}$. Then, we define the non-negative $C^{2}$-function $\tilde{V}: \mathbb{R}_{+}^{4}\rightarrow \mathbb{R}_{+}$ as follows:

Using Itô's formula^[34], we obtain

Similarly, we have

Moreover, we get

and

Therefore, it follows that

Construct the following compact subset

where $\varepsilon_{1}$ is a sufficiently small constant such that

here $D$ $F$, $G$, $H$, $J$ are positive constants, which satisfying

and

For sufficiently small $\varepsilon_{1}$, condition (5.3) holds due to $-pb+B\leq-2$. For convenience of calculation, we divide the set $\mathbb{R}_{+}^{4}\backslash D_{\varepsilon_{1}}$ into eight domains,

then we will illustrate $LV\leq-1$ on the set $\mathbb{R}_{+}^{4}\backslash D_{\varepsilon_{1}}$.

Case 1. If $(S, E, I, Q)\in D_{1}$, we obtain

By inequality (5.2), we obtain $LV\leq-1$ for all $(S, E, I, Q)\in D_{1}$.

Making use of inequalities (5.3), (5.4) and (5.5), and being similar to the proof of Case 1, we make the conclusion that $LV\leq-1$ for all $(S, E, I, Q)\in D_{2}$, $(S, E, I, Q)\in D_{3}$ and $(S, E, I, Q)\in D_{4}$.

Case 2. If $(S, E, I, Q)\in D_{5}$, we have that

According to inequality (5.6), we achieve that $LV\leq-1$ for all $(S, E, I, Q)\in D_{5}$. Furthermore, being similar to the proof of Case 2 and combining with (5.7), (5.8) and (5.9), the conclusion $LV\leq-1$ can be obtained for all $(S, E, I, Q)\in D_{6}$, $(S, E, I, Q)\in D_{7}$ and $(S, E, I, Q)\in D_{8}$.

On the basis of the above discussion, it follows that

which illustrates the condition (H2) of Lemma 1 is satisfied, and model (1.3) has a unique stationary distribution and the ergodicity holds. The proof is completed.

Remark 5.1. Theorem 5.1 reveals that if $\hat{\mathcal{R}}_{0} > 1$, model (1.3) has a unique stationary distribution, that is to say, the disease will persist. Note the express of $\hat{\mathcal{R}}_{0}$, we can see that the smaller the intensity of white noises are, the easier the persistence of the disease is. Especially, when $\sigma_{i} = 0, i = 1, 2, 3, 4$, $\hat{\mathcal{R}}_{0}$ will degenerate into the basic reproduction number $\mathcal{R}_{0}$ of the deterministic model (1.2).

6.   Numerical simulations

In this section, we use two concrete examples to illustrate the obtained the theoretical results. For this purpose, we choose $f(I) = \frac{\beta I}{1+nI}$ with $n = 0.15$ and use the Milstein method mentioned in ^[18] with time $\Delta t = 0.01$, the following discretization system is obtained

where $\tau_{i, k}(i = 1, 2, 3, 4)$ are four independent the Gaussian random variables with $N(0, 1)$.

Example 6.1. Let $\Lambda = 0.4$, $\beta = 0.07$, $\mu = 0.35$, $c = 0.25$, $\varepsilon = 0.5$, $\delta_{1} = 0.05$, $\gamma_{1} = 0.05$, $\alpha = 0.1$, $\delta_{2} = 0.2$, $\gamma_{2} = 0.1$, $\gamma_{3} = 0.03$ and the initial values $(S(0), E(0), I(0), Q(0)) = \left(\frac{5}{7}, 0, \frac{3}{7}, 0\right)$. By a simple calculation, we obtain $\hat{\mathcal{R}}_{0}^{*} = 0.95 < 1$ and $\mu = 0.35 > \frac{\sigma_{1}^{2}\vee\sigma_{2}^{2}\vee\sigma_{3}^{2}\vee\sigma_{4}^{2}}{2} = 0.32$, so the condition of Theorem 4.1 holds. In other word, the disease will go to extinction exponentially, which is shown in Figure 1.

Example 6.2. Taking $\Lambda = 2$, $\beta = 0.5$, $\mu = 0.25$, $c = 0.25$, $\varepsilon = 0.5$, $\delta_{1} = 0.125$, $\gamma_{1} = 0.125$, $\alpha = 0.5$, $\delta_{2} = 0.25$, $\gamma_{2} = 0.25$, $\gamma_{3} = 0.5$ and the initial values $(S(0), E(0), I(0), Q(0)) = (6, 0, 2, 0)$. By the condition of Theorem 5.1, we have $\hat{\mathcal{R}}_{0} = 1.3221 > 1$. Thus, model (1.3) has a unique ergodic stationary distribution as shown in Figure 2 and Figure 3, which means that the disease is persistent.

7.   Conclusions

In this paper, we have investigated the dynamic behavior of a novel stochastic SEIQ disease model with general incidence rate and temporary immunity, which include the bilinear incidence rate ($f(I) = \beta I$) and the saturation incidence rate $\left(f(I) = \frac{\beta I}{1+aI}\right)$. We set up sufficient conditions for extinction of the disease. By using the method of Khasminskii and constructing a suitable Lyapunov function, we prove the existence of a stationary distribution to the model (1.3). From the numerical examples, we observe the stationary distribution of model (1.3) is governed by $\hat{\mathcal{R}_{0}}$, which leads to the persistence of infectious diseases. In addition, the persistence of infectious diseases is also affected by the intensity of environmental noises.

Some interesting and open topics deserve further consideration. For one thing, if the system is affected by other factors, such as time delay or pulse interference, how will the dynamic properties of the system change. For another, we can propose some realistic but complex models, such as considering the influence of Markov switching or lévy jump. We leave these problems as our future work.

Acknowledgements

We will say thanks to reviewers for the valuable ideas and comments on improving the article. The said article is supported by Funding for Outstanding Doctoral Dissertation in NUAA (No. BCXJ18-09), Postgraduate Research & Practice Innovation Program of Jiangsu Province (No. KYCX18_0234), and the Fundamental Research Funds for the central Universities (No. 2021qntd21).

Conflict of interest

The authors declare that there are no conflicts of interest.