Stationary distribution and density function expression for a stochastic SIQRS epidemic model with temporary immunity

Abstract

Recently, considering the temporary immunity of individuals who have recovered from certain infectious diseases, Liu et al. (Phys A Stat Mech Appl 551:124152, 2020) proposed and studied a stochastic susceptible-infected-recovered-susceptible model with logistic growth. For a more realistic situation, the effects of quarantine strategies and stochasticity should be taken into account. Hence, our paper focuses on a stochastic susceptible-infected-quarantined-recovered-susceptible epidemic model with temporary immunity. First, by means of the Khas’minskii theory and Lyapunov function approach, we construct a critical value ${\mathscr {R}}_0^S$ corresponding to the basic reproduction number ${\mathscr {R}}_0$ of the deterministic system. Moreover, we prove that there is a unique ergodic stationary distribution if ${\mathscr {R}}_0^S>1$. Focusing on the results of Zhou et al. (Chaos Soliton Fractals 137:109865, 2020), we develop some suitable solving theories for the general four-dimensional Fokker–Planck equation. The key aim of the present study is to obtain the explicit density function expression of the stationary distribution under ${\mathscr {R}}_0^S>1$. It should be noted that the existence of an ergodic stationary distribution together with the unique exact probability density function can reveal all the dynamical properties of disease persistence in both epidemiological and statistical aspects. Next, some numerical simulations together with parameter analyses are shown to support our theoretical results. Last, through comparison with other articles, results are discussed and the main conclusions are highlighted.

Appendices

Appendix A

(I)Proof of Lemma 2.3: Consider the algebraic equation $G_0^2+A_0\theta _0+\theta _0A_0^{\tau }=0$, where $\theta _0$ is a symmetric matrix. By direct calculation, we have

$$\begin{aligned} \theta _0=\left( \begin{array}{cccc} \sigma _{11}&{} 0 &{}\sigma _{13} &{} 0 \\ 0&{} \sigma _{22} &{} 0 &{}\sigma _{24} \\ \sigma _{13}&{} 0 &{} \sigma _{33} &{} 0 \\ 0&{} \sigma _{24} &{} 0 &{}\sigma _{44} \end{array} \right) , \end{aligned}$$

(6.1)

where $\sigma _{22}=\frac{a_3}{2[a_1(a_2a_3-a_1a_4)-a_3^2]},\ \sigma _{13}=-\sigma _{22},\ \sigma _{33}=\frac{a_1}{a_3}\sigma _{22},\ \sigma _{24}=-\frac{a_1}{a_3},\ \sigma _{11}=\frac{a_2a_3-a_1a_4}{a_3}\sigma _{22}, and\ \sigma _{44}=\frac{a_1a_2-a_3}{a_3a_4}\sigma _{22}$. Assume that $a_1>0,\ a_3>0,\ a_4>0,\ a_1(a_2a_3-a_1a_4)-a_3^2>0$. Then we can show that

$$\begin{aligned}&\sigma _{11}>0,\ \sigma _{11}\sigma _{22}>0,\ \sigma _{22}(\sigma _{11}\sigma _{33}-\sigma _{13}^2)>0,\\&(\sigma _{11}\sigma _{33}-\sigma _{13}^2)(\sigma _{22}\sigma _{44}-\sigma _{24}^2)>0. \end{aligned}$$

This means all the leading principal minors of matrix $\theta _0$ are positive. Consequently, $\theta _0$ is positive definite.

The proof is completed.

(II) Proof of Lemma 2.4: Consider the algebraic equation $G_0^2+B_0\theta _1+\theta _1B_0^{\tau }=0$, where $\theta _1$ is a symmetric matrix. We can get by direct computation that

$$\begin{aligned} \theta _1=\left( \begin{array}{cccc} \theta _{11}&{} 0 &{}\theta _{13} &{} 0 \\ 0&{} \theta _{22} &{} 0 &{}0 \\ \theta _{13}&{} 0 &{} \theta _{33} &{} 0 \\ 0&{} 0 &{} 0 &{}0, \end{array} \right) , \end{aligned}$$

(6.2)

where

$$\begin{aligned}&\theta _{22}=\frac{1}{2(b_1b_2-b_3)},\ \theta _{13}=-\theta _{22},\ \theta _{11}=b_2\theta _{22},\\&\theta _{33}=\frac{b_1}{b_3}\theta _{22}. \end{aligned}$$

If $b_1>0,\ b_3>0,\ b_1b_2-b_3>0$, noting that

$$\begin{aligned} \theta _{11}>0,\ \theta _{11}\theta _{22}>0,\ \theta _{22}(\theta _{11}\theta _{33}-\theta _{13}^2)>0, \end{aligned}$$

which means three leading principal minors of matrix $\theta _1$ are positive. Hence, $\theta _1$ is semi-positive definite. The proof is confirmed.

(III) Proof of Lemma 2.5: For the algebraic equation $G_0^2+C_0\theta _2+\theta _2C_0^{\tau }=0$, since $\theta _2$ is a symmetric matrix, we obtain

$$\begin{aligned} \theta _2=\left( \begin{array}{cccc} \vartheta _{11}&{} 0 &{}0 &{} 0 \\ 0&{} \vartheta _{22} &{} 0 &{}0 \\ 0&{} 0 &{} 0 &{} 0 \\ 0&{} 0 &{} 0 &{}0 \end{array} \right) , \end{aligned}$$

(6.3)

where $\vartheta _{11}=\frac{1}{2c_1},\ \vartheta _{22}=\frac{1}{2c_1c_2}$.

If $c_1>0$ and $c_2>0$, then $\theta _2$ is a semi-positive definite matrix. This completes the proof.

Appendix B (Theory in obtaining standardized transformation matrix)

By means of the invertible linear transformations, we will derive the corresponding standardized transformation matrices of standard $R_1, R_2$, and $R_3$ matrices.

(I) The theory of obtaining standard $R_1$ matrix: For the algebraic equation $G^2+A\Sigma +\Sigma A^{\tau }=0$, where $G=diag(\sigma ,0,0,0)$, and

$$\begin{aligned} A=\left( \begin{array}{cccc} a_{11}&{} a_{12} &{}a_{13} &{} a_{14} \\ a_{21}&{} a_{22} &{} a_{23} &{}a_{24} \\ 0&{} a_{32} &{} a_{33} &{} a_{34} \\ 0&{} 0 &{} a_{43} &{}a_{44} \end{array} \right) . \end{aligned}$$

(6.4)

First, we assume that

$$\begin{aligned} a_{21}\ne 0\text {, } a_{32}\ne 0\text {, } a_{43}\ne 0. \end{aligned}$$

Define $X=(x_1,x_2,x_3,x_4)^{\tau }$ which follows $dX=AXdt$. Considering the following vector $Y=(y_1,y_2,y_3,y_4)^{\tau }$,

$$\begin{aligned} \begin{aligned} y_4&=x_4,\ y_3=y_4'=a_{43}x_3+a_{44}x_4,\\ y_2&=y_3'=a_{43}dx_3+a_{44}dx_4=a_{32}a_{43}x_2+(a_{33}\\&\quad +a_{44})a_{43}x_3+(a_{44}^2+a_{34}a_{43})x_4,\\ y_1&=y_2'=a_{21}a_{32}a_{43}x_1+[(a_{22}\\&\quad +a_{33}+a_{44})a_{32}a_{43}]x_2\\&\quad +[a_{43}(a_{23}a_{32}+a_{34}a_{43}\\&\quad +a_{33}a_{44}+a_{33}^2+a_{44}^2)]x_3\\&\quad + [a_{24}a_{32}a_{43}+(a_{33}+a_{44})a_{34}a_{43} \\&\quad +(a_{34}a_{43}+a_{44}^2)a_{44}]x_4:=m_1x_1\\&\quad +m_2x_2+m_3x_3+m_4x_4. \end{aligned} \end{aligned}$$

Then the corresponding standardized transformation matrix is given by

$$\begin{aligned} M=\left( \begin{array}{cccc} m_1&{} m_2 &{}m_3 &{} m_4 \\ 0&{} a_{32}a_{43} &{} (a_{33}+a_{44})a_{43} &{}a_{44}^2+a_{34}a_{43} \\ 0&{} 0 &{} a_{43} &{} a_{44} \\ 0&{} 0 &{}0 &{}1 \end{array} \right) . \end{aligned}$$

(6.5)

Given the above, we derive that $Y=MX$, which implies that $dY=MdX=MAXdt=(MAM^{-1})Ydt$. Meanwhile, based on the relationship of the vector Y’s components, one has

$$\begin{aligned} dY= d\left( \begin{array}{c} y_1 \\ y_2\\ y_3\\ y_4\end{array}\right) = \left( \begin{array}{cccc} -a_1&{} -a_2 &{}-a_3 &{} -a_4 \\ 1&{} 0 &{} 0 &{}0 \\ 0&{} 1 &{} 0 &{} 0 \\ 0&{} 0 &{} 1 &{} 0 \end{array} \right) \left( \begin{array}{c} y_1 \\ y_2\\ y_3\\ y_4\end{array}\right) dt. \end{aligned}$$

Obviously, we obtain the corresponding standard $R_1$ matrix $MAM^{-1}:=A_0$, which refers to (2.3). Let $\rho _1=a_{21}a_{32}a_{43}\sigma$ and $\theta _0=\rho _1^{-2}M\Sigma M^{\tau }$. Then the above equation can be equivalently transformed into the following equation:

$$\begin{aligned} G_0^2+A_0\theta _0+\theta _0 A_0^{\tau }=0. \end{aligned}$$

(6.6)

(II) The method of transforming standard $R_2$ matrix: For the algebraic equation $G^2+B\Sigma +\Sigma B^{\tau }=0$, where $G=diag(\sigma ,0,0,0)$, and

$$\begin{aligned} B=\left( \begin{array}{cccc} b_{11}&{} b_{12} &{}b_{13} &{} b_{14} \\ b_{21}&{} b_{22} &{} b_{23} &{}b_{24} \\ 0&{} b_{32} &{} b_{33} &{} b_{34} \\ 0&{} 0 &{} 0 &{}b_{44} \end{array} \right) . \end{aligned}$$

(6.7)

Similarly, we stipulate that

$$ b_{21}\ne 0\text {, } b_{32}\ne 0. $$

Let the vector $X=(x_1,x_2,x_3,x_4)^{\tau }$ follow $dX=BXdt$. For the following vector $Y=(y_1,y_2,y_3,y_4)^{\tau }$,

$$\begin{aligned} \begin{aligned} y_4&=x_4,\ \ y_3=x_3,\ \ y_2=y_3'=b_{32}x_2+b_{33}x_3+b_{34}x_4,\\ y_1&=y_2'=b_{32}dx_2+b_{33}dx_3+b_{34}dx_4\\&=b_{21}b_{32}x_1+(b_{22}+b_{33})b_{32}x_2\\&+(b_{33}^2+b_{23}b_{32})x_3+[(b_{33}+b_{44})b_{34}+b_{24}b_{32}]x_4. \end{aligned} \end{aligned}$$

Then the relevant standardized transformation matrix M is described by

$$\begin{aligned} M=\left( \begin{array}{cccc} b_{21}b_{32}&{} (b_{22}+b_{33})b_{32} &{} b_{33}^2+b_{23}b_{32} &{} (b_{33}+b_{44})b_{34}+b_{24}b_{32} \\ 0&{} b_{32} &{} b_{33} &{} b_{34} \\ 0&{} 0 &{} 1 &{} 0 \\ 0&{} 0 &{}0 &{}1 \end{array} \right) . \end{aligned}$$

(6.8)

Moreover, we get that $Y=MX$, which means $dY=MdX=MBXdt=(MBM^{-1})Ydt$. Similarly, according to the relationship of the vector Y’s components, we obtain a standard $R_2$ matrix $MBM^{-1}:=B_0$, which refers to (2.4). Denote $\rho _2=b_{21}b_{32}\sigma$, $\theta _1=\rho _2^{-2}M\Sigma M^{\tau }$. Then the above equation is equivalent to

$$\begin{aligned} G_0^2+B_0\theta _1+\theta _1 B_0^{\tau }=0. \end{aligned}$$

(6.9)

(III) The method of transforming standard $R_3$ matrix: Consider the algebraic equation $G^2+C\Sigma +\Sigma C^{\tau }=0$, where $G=diag(\sigma ,0,0,0)$, and

$$\begin{aligned} C=\left( \begin{array}{cccc} c_{11}&{} c_{12} &{}c_{13} &{} c_{14} \\ c_{21}&{} c_{22} &{} c_{23} &{}c_{24} \\ 0&{} 0 &{} c_{33} &{} c_{34} \\ 0&{} 0 &{} c_{43} &{}c_{44} \end{array} \right) . \end{aligned}$$

(6.10)

First, we assume that $c_{21}\ne 0$. For a vector $X=(x_1,x_2,x_3,x_4)^{\tau }$ determined by $dX=CXdt$, the following vector $Y=(y_1,y_2,y_3,y_4)^{\tau }$ satisfies

$$\begin{aligned} \begin{aligned} y_4&=x_4\text {, }y_3=x_3\text {, } y_2=c_{21}x_1\\&\quad +c_{22}x_2+c_{23}x_3+c_{24}x_4,\\ y_1&=y_2'=c_{21}(c_{11}+c_{22})x_1+(c_{12}c_{21}+c_{22}^2)x_2\\&\quad +[c_{13}c_{21}+c_{23}(c_{22}+c_{33})+c_{24}c_{43}]x_3\\&\quad +[c_{14}c_{21}+c_{24}(c_{22}+c_{44})+c_{23}c_{34}]x_4. \end{aligned} \end{aligned}$$

By defining the corresponding standardized transformation matrix

$$\begin{aligned} M= & {} \left( \begin{array}{cccc} c_{21}(c_{11}+c_{22})&{} c_{12}c_{21}+c_{22}^2 &{} c_{13}c_{21}+c_{23}(c_{22}+c_{33})+c_{24}c_{43} &{} c_{14}c_{21}+c_{24}(c_{22}+c_{44})+c_{23}c_{34} \\ c_{21}&{} c_{22} &{} c_{23} &{} c_{24} \\ 0&{} 0 &{} 1 &{} 0 \\ 0&{} 0 &{}0 &{}1 \end{array} \right) , \end{aligned}$$

(6.11)

we can obtain $Y=MX$. That is to say, $dY=MdX=MCXdt=(MCM^{-1})Ydt$. Hence, the standard $R_3$ matrix $MCM^{-1}:=C_0$ is obtained, which refers to (2.5). Let $\rho _3=c_{21}\sigma$ and $\theta _2=\rho _3^{-2}M\Sigma M^{\tau }$. Then it can be equivalently transformed into the following equation:

$$\begin{aligned} G_0^2+C_0\theta _2+\theta _2 C_0^{\tau }=0. \end{aligned}$$

(6.12)

Appendix C

Consider the following k-dimensional stochastic differential equation

$$\begin{aligned} dX(t)=f(X(t),t)dt+g(X(t),t)dB(t)\qquad \text {for } t\ge t_0, \end{aligned}$$

with the initial value $X(0)=X_0 \in {\mathbb {R}}^k$, where B(t) depicts a k-dimensional standard Brownian motion defined in the above complete probability space. The common differential operator ${\mathscr {L}}$ is described by

$$\begin{aligned} {\mathscr {L}}= & {} \frac{\partial }{\partial t}+\sum _{i=1}^{k}f_i(X(t),t)\frac{\partial }{\partial X_i}\\&+\frac{1}{2}\sum _{i,j=1}^{k}\bigl [g^{\tau }(X(t),t)g(X(t),t)\bigr ]_{ij}\frac{\partial ^2}{\partial X_i\partial X_j}. \end{aligned}$$

Let the operator ${\mathscr {L}}$ act on a function $V\in C^{2,1}({\mathbb {R}}^k\times [t_0,\infty ];{\mathbb {R}}_+^1)$. Then one can determine that

$$\begin{aligned}&{\mathscr {L}}V(X,t)=V_t(X(t),t)+V_X(X(t),t)f(X(t),t)\\&\quad +\frac{1}{2}trace\bigl [g^{\tau }(X(t),t)V_{XX}(X(t),t)g(X(t),t)\bigr ], \end{aligned}$$

where $V_t=\frac{\partial V}{\partial t}$, $V_X=(\frac{\partial V}{\partial x_1},\ldots ,\frac{\partial V}{\partial x_k})$, and $V_{XX}=(\frac{\partial ^2V}{\partial x_i \partial x_j})_{k\times k}$. If $X(t) \in {\mathbb {R}}^k$, we have

$$\begin{aligned}&dV(X(t),t)={\mathscr {L}}V(X(t),t)dt\\&\quad +V_X(X(t),t)g(X(t),t)dB(t). \end{aligned}$$