Stochastic modeling for an SIVR coronavirus epidemic model with distributed delay: dynamic properties and optimal control

Abstract

This paper develops a stochastic SIVR epidemic model with distributed delay to analyze the dynamics of coronavirus spread. We first derive the basic reproduction number and equilibrium points for the deterministic version of the model. For the stochastic model, we establish sufficient conditions for disease persistence using a Lyapunov approach and define a critical threshold $\mathcal {R}_0^c$, proving the existence of a unique stationary distribution via Khasminskii theory. Moreover, we obtain the exact probability density function around the quasi-equilibrium by solving the associated Fokker–Planck equation. To support intervention planning, we formulate a stochastic optimal control problem with three targeted strategies for susceptible, vaccinated, and infected subpopulations and characterize optimal controls using the stochastic maximum principle. These analytical results clarify key mechanisms driving disease persistence from an epidemiological standpoint. Theoretical findings are validated through numerical simulations, where key parameters are estimated via the Markov chain Monte Carlo method using real-time data.

Appendix

1.1 Proof of Theorem 2.1

Proof

Observing the coefficient of Model (5), we find it satisfies the local Lipschitz condition. Hence, there exists a unique local solution for $t \in [0, \tau _e)$, with $\tau _e$ denoting the explosion time. Utilizing Itô’s formula, it can be demonstrated that the singular local solution of (5) maintains positivity. Next, we aim to demonstrate that this solution is global, implying $\tau _e = \infty$ almost certainly.

Let $n_0 > 0$ be sufficiently large such that the initial values of S(0), I(0), V(0), and U(0) lie within the interval $[\frac{1}{n_0}, n_0]$. For each integer $n \ge n_0$, we define a sequence of stopping times as follows:

$$\begin{aligned} \tau _n = \inf \left\{ t \in [0, \tau _e]: S(t) \notin (\frac{1}{n}, n ) \text { or } I(t) \notin (\frac{1}{n}, n ) \text { or } V(t) \notin (\frac{1}{n}, n ) \text { or } U(t) \notin (\frac{1}{n}, n )\right\} , \end{aligned}$$

where we define $\inf \emptyset = \infty$ (with $\emptyset$ representing the empty set). Since the sequence $\{\tau _n\}$ is nondecreasing as $n \rightarrow \infty$, there exists a limit given by

$$\begin{aligned} \tau _\infty = \lim _{n\rightarrow \infty }\tau _n. \end{aligned}$$

Next, our objective is to demonstrate that $\tau _\infty = \infty$ almost surely. If this assertion were false, it would imply that there exist $T > 0$ and $\varepsilon \in (0,1)$ such that

$$\begin{aligned} \mathcal {P} (\tau _\infty \le T ) > \varepsilon . \end{aligned}$$

Thus, there is an integer $n_1 \ge n_0$ that satisfies the inequality:

$$\begin{aligned} \mathcal {P} (\tau _n \le T ) \ge \varepsilon , \quad \forall n \ge n_1. \end{aligned}$$

(56)

Let a and b denote two positive constants, and let $\mathcal {V}: \mathbb {R}_+^4 \rightarrow \mathbb {R}_+$ be a nonnegative $\mathcal {C}^4$ function defined as follows:

$$\begin{aligned} \mathcal {V}(S, I, V, U)=\left( S-a-a \ln \frac{S}{a}\right) + ( I-1-\ln I)+(V-1-\ln V)+ b(U-1-\ln U). \end{aligned}$$

Since $-\ln x$ approaches $\infty$ as x tends to $0^+$, and $(-x-m-m) \ln \frac{x}{m}$ also tends to $\infty$ for any $m > 0$ as x approaches $0^+$, it follows that

$$\begin{aligned} \mathop {\textrm{lim inf}}\limits _{n\rightarrow \infty ,(S,I,V,U)\in \mathbb {R}_{+}^4\backslash \mathbb {D}_n}\mathcal {V} \left( S,I,V,U\right) =\infty , \end{aligned}$$

where $\mathbb {D}_n=(\frac{1}{n},n)\times (\frac{1}{n},n)\times (\frac{1}{n},n)\times (\frac{1}{n},n)$ for each integer $n \ge n_0$.

By applying Itô’s formula to $\mathcal {V}$, we derive:

$$\begin{aligned} \textrm{d}\mathcal {V}&=\bigg ((1-\frac{a}{S})(\Lambda -\beta C_0 S U - (\omega _1+\mu )S)+(1-\frac{1}{I})(\beta C_0 S U +\gamma _0 V U -(\delta +\mu +r_1) I) \\&\quad +(1-\frac{1}{V})(\omega _1 S - \gamma _0 V U -\mu V)+ \alpha b(1-\frac{1}{U})( I- U)+\frac{a}{2}\sigma _1^2 +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2}\bigg )\textrm{d}t \\&\quad +\sigma _1 S(1-\dfrac{ a }{ S })\textrm{d} B_1(t)+\sigma _2 I (1-\dfrac{1}{ I })\textrm{d} B_2(t)+\sigma _3 V (1-\dfrac{1}{ V })\textrm{d} B_3(t) \\&=\bigg (\Lambda +a(\omega _1+\mu )+(\delta +\mu +r_1)+\mu +\alpha b +\frac{a \sigma _1^2}{2} +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2}+(\alpha \beta C_0 +\gamma _0- \alpha b)U \\&\quad +(\alpha b -(\delta +\mu +r_1))I-(\mu (S+V) +a\frac{\Lambda }{S} +\beta C_0 \frac{SU}{I}+ \gamma _0 \frac{VU}{I} +\omega _1 \frac{S}{V} + \alpha b\frac{I}{U})\bigg ) \textrm{d}t \\&\quad +\sigma _1 S(1-\dfrac{ a }{ S })\textrm{d} B_1(t)+\sigma _2 I(1-\dfrac{1}{ I })\textrm{d} B_2(t)+\sigma _3 V(1-\dfrac{1}{ V })\textrm{d} B_3(t)\\&=\mathscr {L}\mathcal {V}\textrm{d}t + \sigma _1 S(1-\dfrac{ a }{ S })\textrm{d} B_1(t)+\sigma _2 I(1-\dfrac{1}{ I })\textrm{d} B_2(t)+\sigma _3 V(1-\dfrac{1}{ V })\textrm{d} B_3(t), \end{aligned}$$

where

$$\begin{aligned} \mathscr {L}\mathcal {V}&=(1-\frac{a}{S})(\Lambda -\beta C_0 S U - (\omega _1+\mu )S)+(1-\frac{1}{I})(\beta C_0 S U +\gamma _0 V U -(\delta +\mu +r_1) I) \\&\quad +(1-\frac{1}{V})(\omega _1 S - \gamma _0 V U -\mu V)+ \alpha b(1-\frac{1}{U})( I- U)+\frac{a}{2}\sigma _1^2 +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2} \\&=(\Lambda +a(\omega _1+\mu )+(\delta +\mu +r_1)+\mu +b \alpha +\frac{a\sigma _1^2}{2} +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2})+(\alpha \beta C_0 +\gamma _0- \alpha b)U \\&\quad +(\alpha b-(\delta +\mu +r_1))I-(\mu (S+V) +a\frac{\Lambda }{S} +\beta C_0 \frac{SU}{I}+ \gamma _0 \frac{VU}{I} +\omega _1 \frac{S}{V} + \alpha b\frac{I}{U}). \end{aligned}$$

Choosing $a=\frac{\delta +\mu +r_1-\gamma _0}{\beta C_0}$ and $b=\frac{\delta +\mu +r_1}{\alpha }$ yields

$$\begin{aligned} \mathscr {L}\mathcal {V}&=(\Lambda +a(\omega _1+\mu )+(\delta +\mu +r_1)+\mu + \alpha b+\frac{a\sigma _1^2}{2} +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2}) \\&\quad -(\mu (S+V) +a\frac{\Lambda }{S} +\beta C_0 \frac{SU}{I}+ \gamma _0 \frac{VU}{I} +\omega _1 \frac{S}{V} + \alpha b\frac{I}{U}) \\&< \Lambda +2\mu +\delta +r_1 +a(\omega _1+\mu ) +\alpha b+\frac{a \sigma _1^2}{2} +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2} \\&=\Lambda +3\mu +2(\delta +r_1) +\frac{(\omega _1+\mu )(\delta +\mu +r_1-\gamma _0)}{\beta C_0}+\frac{(\delta +\mu +r_1-\gamma _0)}{2\beta C_0}\sigma _1^2 +\frac{\sigma _2^2}{2}+\frac{\sigma _3^2}{2} \\&:= K, \end{aligned}$$

where $K$ represents a positive constant. In this case, we observe that

$$\begin{aligned} \textrm{d}\mathcal {V}(S,I,V,U)\le K \textrm{d}t+\sigma _1 S(1-\dfrac{ a }{ S })\textrm{d} B_1(t)+\sigma _2 I(1-\dfrac{1}{ I })\textrm{d} B_2(t)+\sigma _3 V(1-\dfrac{1}{ V })\textrm{d} B_3(t). \end{aligned}$$

Integrating both sides of the inequality above over the interval from 0 to $\tau _n\wedge T$, we obtain:

$$\begin{aligned} \int _0^{\tau _n\wedge {T}}\textrm{d}\mathcal {V}&\le \int _0^{\tau _n\wedge {T}}K\textrm{d}t+\int _0^{\tau _n\wedge {T}} \sigma _1 S(1-\frac{ a }{ S })\textrm{d} B_1(s)+\int _0^{\tau _n\wedge {T}}\sigma _2 I(1-\frac{1}{ I })\textrm{d} B_2(s)\\&+\int _0^{\tau _n\wedge {T}} \sigma _3 V(1-\frac{1}{ V })\textrm{d} B_3(s), \end{aligned}$$

where $\tau _n\wedge T=\min \{ \tau _n, T \}$. Taking the expectations of both sides, we get:

$$\begin{aligned} \textbf{E}\mathcal {V} (S (\tau _n\wedge T), I(\tau _n\wedge T ), V(\tau _n\wedge T ) , U(\tau _n\wedge T )) \le \mathcal {V}(S(0), I(0), V(0), U(0))+KT. \end{aligned}$$

Let $\Omega _n=\left\{ \tau _n\le T \right\}$ for $n \ge n_1$. From Eq. (56), we have $\mathcal {P}(\Omega )\ge \varepsilon$. For every $v\in \Omega _n$, there exists some i such that $x_i(\tau _n,v)$ equals either n or $\dfrac{1}{n}$ for $i=1,2,3,4$. Therefore, $\mathcal {V} (S (\tau _n, v ),I (\tau _n, v ),V (\tau _n, v ), U (\tau _n, v ) )$ is no less than $\min \{ n-1-\ln n,\frac{1}{n}-1-\ln \frac{1}{n}\}$. Thus, we obtain

$$\begin{aligned} \mathcal {V} (S(0), I(0), V(0), U(0))+KT \ge \textbf{E} (\textbf{I}_{\Omega _{n}(v)}\mathcal {V} (S(\tau _{n}),I(\tau _{n}),V(\tau _{n}),U(\tau _{n}) )) \ge \varepsilon \theta _n, \end{aligned}$$

where

$$\begin{aligned} \theta _{n}=\min \{n-a- \ln \frac{n}{a}, \frac{1}{n}-a+a \ln (na), n-1- \ln n, \frac{1}{n} -1 +\ln n,b(n-1-\ln n), b(\frac{1}{n} -1 +\ln n) \}, \end{aligned}$$

where $\textbf{I}_{\Omega _{n}(v)}$ is the indicator function of $\Omega _{n}$. Letting $n\rightarrow \infty$ results in the contradiction $\infty =\mathcal {V}(S(0),I(0),V(0),U(0))+KT<\infty$. This completes the proof. $\square$