Dynamical analysis of a stochastic SIRS epidemic model with saturating contact rate

1.   Introduction

Recently, a new type of pneumonia caused by the coronavirus, named COVID-19, is spreading around the world. The issue of infectious diseases has once again aroused people's great concern. How to prevent and control infectious diseases has been an important subject facing human beings^{[1,2,3,4,5,6,7,8]}. The SIR model assumes that the infected person can obtain permanent immunity after recovery. However, for smallpox, cholera, malaria and other diseases, individuals recovered from treatment can return to the susceptible category after temporary immunization, which can be described by SIRS model. Moreover, for some bacterial infectious diseases, such as meningitis and sexually transmitted diseases, some individuals can not produce effective antibodies after treatment and may be infected again. Others may gain temporary immunity, but then lose immunity and become susceptible ^[9,10,11,12]. Literature ^[10] established an SIRS model with a general population-size dependent contact rate $\lambda (N)$ and proportional transfer rate from the infective class to susceptible class. The authors studied the threshold conditions of disease extinction and discussed the stability of disease-free equilibrium and endemic equilibrium.

Infection rate is an important index to measure the intensity of disease transmission. Employing an appropriate infection rate based on a specific disease for the mathematical model plays a vital role in the disease prevention and control. In the literature ^[13], Thieme and Castillo-Chavez proposed the incidence $\frac{{\beta \varrho \left({N\left(t \right)} \right)S\left(t \right)I\left(t \right)}}{{N\left(t \right)}}$, where $N\left(t \right)$ represents the total population. On that basis, Heesterbeek et al. ^[14] gave the following saturating contact rate

Obviously, $\varrho\left({N\left(t \right)} \right)$ is a non-decreasing function of $N\left(t \right)$. $\frac{{\varrho\left({N\left(t \right)} \right)}}{{N\left(t \right)}}$ is a non-increasing function of $N\left(t \right)$. If $N\left(t \right)$ is sufficiently small, $\varrho\left({N\left(t \right)} \right) \sim bN$. Conversly, if $N\left(t \right)$ is fully large, $\varrho\left({N\left(t \right)} \right) \sim 1$. Compared with the bilinear incidence $\beta S\left(t \right)I\left(t \right)$ and the standard incidence $\frac{{\beta S\left(t \right)I\left(t \right)}}{{N\left(t \right)}}$, the saturating contact rate is more closer to the transmission of many diseases. The saturated contact rate is widely used in the study of infectious disease modeling. For example, Zhang et al.^[15] constructed an SEIS model with general saturated incidence rate, and proved the global asymptotic stability of the endemic equilibrium by using the autonomous convergence theorem. Lan et al.^[16] considered an SIS epidemic model with saturating contact rate, by using Itô's formula, the conditions for disease extinction and the existence of stationary solutions were obtained. In reference ^[11], Li et al. established an SIRS epidemic model with a general incidence, which considered both the transfer from the infected to the susceptible and the transfer from the recovered to the susceptible. Motivated by the above literature, we formulate a deterministic SIRS epidemic model with saturating contact rate and transfer from infectious to susceptible:

where $g\left({N\left(t \right)} \right) = 1 + bN\left(t \right) + \sqrt {1 + 2bN\left(t \right)}$, $N\left(t \right) = S\left(t \right) + I\left(t \right) + R\left(t \right)$ is the total population. $S\left(t \right)$, $I\left(t \right)$, $R\left(t \right)$ represent the number of susceptible individuals, infected individuals and recovered individuals, respectively. $\Lambda$ is the recruitment rate of susceptible individuals. $u$ denotes the natural mortality rate. $\alpha$ represents the mortality rate caused by diseases. ${\gamma _1}$ is the transfer rate from the infected individuals to the susceptible individuals. ${\gamma _2}$ is the transfer rate from the infected individuals to the recovered individuals, and $\delta$ denotes the immunity loss rate.

Due to the influence of environmental noise, the prevalence and transmission of diseases is often random. For example, the change of temperature and the influence of climate will lead to the fluctuation of mortality, morbidity and so on. In recent years, mathematical models of infectious diseases described by stochastic differential equations have been widely concerned ^{[17,18,19,20,21,22]}. There are many ways to construct a stochastic differential equation model, such as adding random perturbations to the parameters of deterministic system ^{[23,24,25,26]}, or introducing proportional perturbations to state variables ^{[27,28,29,30,31]}. Recently, considering the effect of two different white noises on the model parameters, reference ^[32] established a stochastic SIS model with two correlated Brownian motions, in which the threshold of disease extinction as well as the variance and mean of the stationary distribution were investigated. In this paper, we consider that the incidence coefficient $\beta b$ is disturbed by white noise, that is $\beta b \to \beta b + \sigma dB\left(t \right)$, where ${\sigma ^2}$ is the intensity of white noise, $B\left(t \right)$ is defined as the standard Brownian motion in the complete probability space $(\Omega, \mathcal{F}, \{\mathcal{F}\}_{t\geq0}, P)$. Thus, the above model (1.1) is transformed into the following SIRS stochastic epidemic model:

As far as we know, there have been a lot of studies on the epidemic model disturbed by environmental noise, but there are few stochastic models considering saturating contact rate and transfer from infectious to susceptible, especially the existence of stationary solution. The paper is organized as follows: In section 2, we give some notations and related lemmas. The thresholds of deterministic system and stochastic system are established in sections 3 and 4, respectively. Sufficient condition for the existence of stationary solution in the stochastic system (1.2) is provided in section 5. In section 6, we verify the results of theoretical derivation by numerical simulations.

2.   Preliminaries

Throughout this paper, we let $R_ + ^3 = \left\{ {{x_i} > 0, i = 1, 2, 3} \right\}$. For an integrable function $h$ on $[0, +\infty)$, we define $\langle h(t)\rangle = \frac{1}{t}\int^{t}_{0}h(\pi) \rm{d}\pi$. By using the methods from Liu et al. ^[33], we can prove that the region

is a positively invariant set of system (1.2).

Lemma 2.1. For any given initial value $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right)} \right) \in R_ + ^3$, then the model (1.2) has a unique positive solution $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right)$ on $t \ge 0$, and the solution will remain in $R_ + ^3$ with probability 1.

Next, we will introduce some contents of stationary Markov process. The $n$-dimensional stochastic differential equation can be expressed by the following formula

with the initial value $x\left({{t_0}} \right) = {x_0} \in \mathbb{R}{^n}$. Integrating from $0$ to $t$ for both sides of the Eq (2.1), one can obtain that

Assume that the vectors $f\left({x, t} \right)$, ${g_1}\left({x, t} \right)$, ..., ${g_l}\left({x, t} \right)$$\left({t \ge {t_0}, x \in \mathbb{R}{^n}} \right)$ are continuous functions of $\left({x, t} \right)$, satisfying the following conditions for some constant $D$,

According to the literature ^[34] Theorem 3.7, we have the following Lemma.

Lemma 2.2. Suppose that the coefficients of (2.2) are independent of $t$, and the conditions (2.3) hold in ${U_G}\left({\forall G > 0} \right)$. There exists a function $V\left(x \right) \in {C^2}$ with the following properties in $\mathbb{R}{^n}$

where ${C^2}$ represents a class of functions that are twice continuously differentiable relative to $x$ in $\mathbb{R}{^n}$. Further, we assume that there is at least one $x \in \mathbb{R}{^n}$, such that the process ${X^x}\left(t \right)$ is regular. Then there exists a solution of system (2.2) which is a stationary Markov process.

For Lemma 2.2, we need to note the following two points.

Remark 2.1 (ⅰ) Condition (2.3) can be replaced by the global existence of solution of system (2.2) (see ^[35] Remark 5);

(ⅱ) Condition (2.4) can be replaced by the weaker condition $LV\left(x \right) \le - 1$ (see ^[34] Chapter 4).

3.   Dynamics of system (1.1)

In this section, we concentrate on the stability of the equilibria. Firstly, using $N\left(t \right)$ as a variable instead of the variable $S\left(t \right)$, we convert system (1.1) into the following form

Obviously, the system (3.1) exists a boundary equilibrium ${P_0} = \left({\frac{\Lambda }{u}, 0, 0} \right)$. Define

By direct calculation, if ${R_0} > 1$, we get system (3.1) has a positive equilibrium ${P^*} = \left({{N^ * }, {I^ * }, {R^ * }} \right)$ with

where ${{N^ * }}$ is the unique positive root of the following function

In fact, we have

and

Let $\frac{{d\phi \left(N \right)}}{{dN}} = 0$, it can be got

which implies that $\phi \left(N \right)$ is increasing if $N \ge {N_{ * * }}$, and $\phi \left(N \right)$ is decreasing if $N < {N_{ * * }}$. Therefore, function $\phi \left(N \right)$ has a unique positive root, then the system (3.1) exists a unique positive equilibrium ${P^ * } = \left({{N^ * }, {I^ * }, {R^ * }} \right)$.

Theorem 3.1. For system (3.1), we have

(i) If ${R_0} < 1$, then ${P_0} = \left({\frac{\Lambda }{u}, 0, 0} \right)$ is a unique stable equilibrium, which implies the disease of system (3.1) goes extinct.

(ii) If ${R_0} > 1$, then ${P^ * } = \left({{N^ * }, {I^ * }, {R^ * }} \right)$ is a stable positive equilibrium, which implies the disease of system (3.1) is permanent.

Proof. (ⅰ) The Jacobian matrix of system (3.1) evaluated at ${P_0} = \left({\frac{\Lambda }{u}, 0, 0} \right)$ is

which has three eigenvalues:

Therefore, according to stability theory, ${P_0}$ is stable if ${R_0} < 1$.

(ⅱ) The Jacobian matrix at ${P^ * } = \left({{N^ * }, {I^ * }, {R^ * }} \right)$ is

where

Hence, the characteristic equation of ${J^ * }$ is

where

and

Then,

Therefore, the equilibrium ${P^ * }$ is stable when it exists. This completes the proof of Theorem 3.1.

4.   The threshold of the system (1.2)

In the previous section, we have obtained the threshold for ordinary differential equation (ODE) system (1.1). Similarly, the threshold of stochastic differential equation (SDE) system (1.2) is also crucial, which determines the extinction and persistence of the disease. {Define the following parameter

where $g\left(x \right) = 1 + bx + \sqrt {1 + 2bx}$. In this section, we will prove that $R_0^s$ is the threshold of system (1.2). Now, we give the following definition.

Definition 4.1.

(i) The disease $I(t)$ is said to be extinct if $\mathop {\lim }\limits_{{\rm{t}} \to + \infty } I(t) = 0$;

(ii) The disease $I(t)$ is said to be permanent in mean if there is a positive constant $\phi$ such that $\mathop {\lim }\limits_{{\rm{t}} \to \infty } \inf \langle {I(t)} \rangle \ge \phi$.

4.1.   Extinction

Theorem 4.1. Set $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right)$ be a solution of system (1.2) with any given initial value $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right)} \right) \in \Gamma$.

$(i)$ If ${\sigma ^2} > \max \left\{ {\frac{{{\beta ^2}{b^2}}}{{2\left({u + {\gamma _1} + {\gamma _2} + \alpha } \right)}}, \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda }} \right\}$ holds, then

$(ii)$ If $R_0^s < 1$ and ${\sigma ^2} < \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda }$ holds, then

which implies that the disease dies out with probability 1. In addition, we have

Proof. Set $\mathop V\limits^ \sim = \ln I\left(t \right)$, by the Itô's formula, one can get

where $\Phi \left(x \right) = - \frac{{{\sigma ^2}}}{2}{x^2} + \beta bx - \left({u + {\gamma _1} + {\gamma _2} + \alpha } \right)$.

Case (ⅰ): If condition (ⅰ) is satisfied, then

Substituting (4.2) into (4.1), one can obtain that

By using the strong law of large numbers for martingales, we obtain

Integrating from $0$ to $t$, dividing by $t$, and taking superior limit on both sides of Eq (4.3) yields

Case (ⅱ): If condition (ⅱ) $R_0^s < 1$ and ${\sigma ^2} < \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda }$ are satisfied, then

From Eq (4.1), one can get

The inequalities (4.4) and (4.6) imply $\mathop {\lim }\limits_{t \to \infty } I\left(t \right) = 0$ and the disease goes to extinction.

Next, adding up the three equations of system (1.2), integrating both sides from $0$ to $t$ and dividing by t, one can see that

So, from Eq (4.7) we have

where $\Theta \left(t \right) = \frac{1}{u}\left[ {\frac{{S\left(t \right) - S\left(0 \right)}}{t} + \frac{{I\left(t \right) - I\left(0 \right)}}{t} + \frac{\delta }{{u + \delta }}\frac{{R\left(t \right) - R\left(0 \right)}}{t}} \right]$. Obviously, $\mathop {\lim }\limits_{t \to \infty } \Theta \left(t \right) = 0$ and we have

Similarly, from the third equation of the system (1.2) yields

Let $t \to \infty$, we get that

This finishes the proof.

4.2.   Permanence in mean

Theorem 4.2. If $R_0^s > 1$, the disease is persistence in mean. Moreover, we have

where

Proof. In view of Eq (4.1), we have

Integrating from $0$ to $t$ and dividing by $t$ on both sides of Eq (4.10), one can get that

Substituting Eq (4.8) into Eq (4.11), we can get

Hence,

where $F\left(t \right) = \ln I\left(0 \right) + \int_0^t {\frac{{\sigma S\left(\xi \right)}}{{g\left({N\left(\xi \right)} \right)}}d\xi } - \frac{{\beta bt\Theta \left(t \right)}}{{ug\left({\frac{\Lambda }{u}} \right)}}$. Obviously, $\mathop {\lim }\limits_{t \to \infty } \frac{{F\left(t \right)}}{t} = 0\quad a.s.$. From Lemma 1 in ^[12], we obtain

According to Eqs (4.8) and (4.9), we can see that

and

The proof of Theorem 4.2 is completed.

Remark 4.1. According to Theorem 4.1 (i), if the intensity of white noise is large enough that the condition ${\sigma ^2} > \max \left\{ {\frac{{{\beta ^2}{b^2}}}{{2\left({u + {\gamma _1} + {\gamma _2} + \alpha } \right)}}, \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda }} \right\}$ holds, then the disease goes to extinction. Therefore, a large environmental nose intensity can suppress the spread of disease. In addition, by comparing the thresholds of systems (1.2) and (1.1), it can be found that if the intensity of environmental noise ${\sigma ^2}{\rm{ = }}0$, then $R_0^s = {R_0}$; if ${\sigma ^2} \ne 0$, then $R_0^s < {R_0}$. When $R_0^s < 1 < {R_0}$, the deterministic system (1.1) has a stable positive equilibrium, while the disease of the stochastic system (1.2) dies out with probability 1. This means that the presence of environmental noise is conducive to disease control.

5.   The existence of stationary solution

Theorem 5.1. If $R_0^s > 1$, there exists a solution of system (1.2) which is a stationary Markov process.

Proof. Let $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right)$ be a solution of system (1.2) with any given initial value $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right)} \right) \in \Gamma$. Then we construct a ${C^2}$-function $G$ as following:

where ${V_1} = - \ln I\left(t \right) - {c_1}\ln S\left(t \right) - {c_2}{g^2}\left({N\left(t \right)} \right) + \frac{{4b{c_2}g\left({\frac{\Lambda }{u}} \right)}}{{u + \delta }}R\left(t \right)$, ${V_2} = - \ln S\left(t \right)$, ${V_3} = - \ln R\left(t \right)$, ${V_4} = - \ln \left({N\left(t \right) - \frac{\Lambda }{{u + \alpha }}} \right)$, ${V_5} = - \ln \left({\frac{\Lambda }{u} - N\left(t \right)} \right)$, ${c_1} = \frac{{u + {\gamma _1} + {\gamma _2} + \alpha + \frac{{{\Lambda ^2}{\sigma ^2}}}{{2{u^2}{g^2}\left({\frac{\Lambda }{u}} \right)}}}}{u}$ and ${c_2} = \frac{{u + {\gamma _1} + {\gamma _2} + \alpha + \frac{{{\Lambda ^2}{\sigma ^2}}}{{2{u^2}{g^2}\left({\frac{\Lambda }{u}} \right)}}}}{{2b\Lambda g\left({\frac{\Lambda }{u}} \right)}}$. The $M$ is a positive constant and satisfies the following condition

Furthermore, $G\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right)$ is a continuous function, which exists a minimum point $\left({\mathop {{S_0}}\limits^ -, \mathop {{I_0}}\limits^ -, \mathop {{R_0}}\limits^ - } \right)$. Next, we define a nonnegative ${C^2}$-function $V$

Applying Itô's formula for ${V_1}$, we can see

Similarly, we have

Therefore, in view of Eqs (5.2)–(5.6), we get that

where $\lambda = \frac{{{c_1}\beta b}}{{g\left({\frac{\Lambda }{{u + \alpha }}} \right)}} + 4b{c_2}g\left({\frac{\Lambda }{u}} \right)\left({\alpha + \frac{{{\gamma _2}}}{{u + \delta }}} \right)$.

Define the following bounded closed set

where $\varepsilon$ is a sufficiently small constant satisfying the following inequalities (5.7)–(5.11)

For convenience, we divide $\Gamma \backslash {D_\varepsilon }$ into five domains

So, we only need to prove $LV \le - 1$ on the above five domains.

Case ⅰ: If $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right) \in {D_1}$, from Eq (5.7), one can obtain that

Case ⅱ: If $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right) \in {D_2}$, in view of Eq (5.8), we have

Case ⅲ: If $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right) \in {D_3}$, according to Eq (5.9), we obtain

Case ⅳ: If $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right) \in {D_4}$, by Eq (5.10), one can see that

Case ⅴ: If $\left({S\left(t \right), I\left(t \right), R\left(t \right)} \right) \in {D_5}$, in view Eq (5.11), we derive

Consequently, for a sufficiently small $\varepsilon$, we have

In view of Lemma 2.2, there exists a solution of system (1.2) which is a stationary Markov process. This completes the proof.

Remark 5.1. Theorem 5.1 shows that if the intensity of white noise is small enough to make $R_0^s > 1$, then the system (1.2) has a stationary solution. That is to say, disease will exist for a long time and form endemic disease.

6.   Conclusions and numerical simulations

In this paper, we study the dynamics of a stochastic SIRS epidemic model with saturating contact rate. Firstly, the threshold of disease extinction for deterministic system (1.1) is obtained by using the Jacobian matrix. If ${R_0} < 1$, system (1.1) has a unique stable equilibrium and the disease goes to extinct. If ${R_0} > 1$, system (1.1) forms endemic disease after a sufficiently long time. Secondly, the threshold parameter $R_0^s$ of stochastic system (1.2) is established. If the intensity of the white noise is small enough to satisfy the condition $R_0^s > 1$, the disease is persistent. Otherwise, the disease will die out. Finally, we prove that there exists a stationary solution under condition $R_0^s > 1$ in system (1.2).

In order to demonstrate the above theoretical derivation, we use MATLAB software to carry out some numerical simulations. Next, we choose the relevant parameters of system (1.1) as follows:

Firstly, we simulate the deterministic system and set $\beta = 0.66$. By an ordinary computation, ${R_0}\; = \; 0.9335 < 1$. From the first condition of the Theorem 3.1, one can obtain that the system (1.1) has a unique stable equilibrium point ${E_0} = (10, 0, 0)$ (Figure 1). Furthermore, if we increase $\beta$ to $0.8$, in this case, we have ${R_0} = 1.1315 > 1$. The condition (ⅱ) in Theorem 3.1 is established, then by Theorem 3.1, the disease of system (1.1) is persistent (Figure 2).

Next, we perform numerical simulations on stochastic system. We keep the parameters of deterministic system (1.1) unchanged, and only select different intensities of white noise $\sigma$ in system (1.2).

Case ⅰ: In order to verify the conclusion (ⅰ) of the Theorem 4.1, we let $\sigma = 0.8$. By computing, we can obtain ${\sigma ^2} = 0.64$, $g\left({\frac{\Lambda }{u}} \right) = 10.6056$ and $\max \left\{ {\frac{{{\beta ^2}{b^2}}}{{2\left({u + {\gamma _1} + {\gamma _2} + \alpha } \right)}}, \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda }} \right\} = 0.5091$, which implies the parameters satisfy the condition (ⅰ) ${\sigma ^2} > \max \left\{ {\frac{{{\beta ^2}{b^2}}}{{2\left({u + {\gamma _1} + {\gamma _2} + \alpha } \right)}}, \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda }} \right\}$. This shows that the disease in system (1.2) dies out with probability 1 (Figure 3).

Case ⅱ: In Figure 4, we assume that $\sigma = 0.45$. It is not difficult to obtain that $\sigma$ satisfies the second condition of Theorem 4.1. {Then, we have $R_0^s = 0.9236 < 1$ and $0.2025 = {\sigma ^2} < \frac{{\beta bug\left({\frac{\Lambda }{u}} \right)}}{\Lambda } = 0.5091$.} As can be seen from Figure 4, the disease is going extinct.

Case ⅲ: If $\sigma = 0.1$, by calculation, we get $R_0^s = 1.1190 > 1$. According to Theorem 4.2, the disease in system (1.2) is permanent in the time mean (Figure 5). Furthermore, by Theorem 5.1, the system (1.2) has a stationary solution (Figure 6).

Acknowledgments

This work is supported by Shandong Provincial Natural Science Foundation of China (No. ZR2019MA003).

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this paper.