A comparative study of deterministic and stochastic computational modeling approaches for analyzing and optimizing COVID-19 control

Abstract

This paper presents a comparative analysis of deterministic and stochastic computational modeling approaches for the optimal control of COVID-19. We formulate a compartmental epidemic model with perturbation by white noise that incorporates various factors influencing disease transmission. By incorporating stochastic effects, the model accounts for uncertainties inherent in real-world epidemic data. We establish the mathematical properties of the model, such as well-posedness and the existence of stationary distributions, which are crucial for understanding long-term epidemic dynamics. Moreover, the study presents an optimal control strategies to mitigate the epidemic’s impact, both in deterministic and stochastic sceneries. Reported data from Algeria are used to parameterize the model, ensuring its relevance and applicability to practical satiation. Through numerical simulations, the study provides insights into the effectiveness of different control measures in managing COVID-19 outbreaks. This research contributes to advancing our understanding of epidemic dynamics and informs decision-making processes for epidemic controlling interventions.

Introduction

Mathematical modeling and simulation are essential tools for informed decision-making in epidemic control^1,2. Each epidemic exhibits a unique biological characteristics, necessitating the formulation of dynamical models to accurately capture the transmission processes and address real-world scenarios. The highly contagious Coronavirus infection which first arose in China in December 2019 and spread quickly to cover practically the whole world posing a serious threat to everyone on the globe³. The COVID-19 epidemic has created an unparalleled worldwide health crisis causing great human suffering and societal upheaval. The virus has been the target of multiple efforts to stop its spread, including non-pharmaceutical therapies, vaccination drives, and other public health measures. Mathematical models help to explore how a disease spreads through a population and suggest effective control interventions. These models can be applied to analyze public health measures’ efficacy and forecast the progression of epidemic. It’s a potent instrument that may be utilized to comprehend the spread of illness and formulate control methods^4,5.

Although, several deterministic epidemic models have been discussed in previous studies^6,7. Stochastic differential models are more effective than deterministic models for simulating biological phenomena, as they offer a higher level of realism^8,9,10. Running a stochastic model multiple times generates a distribution of possible outcomes, providing deeper and more practical insights. In contrast, deterministic models yields only a single predicted value^11,12. The transmission of infectious disease is inherently random, making it essential to incorporate environmental noise into models to capture environmental fluctuations. Different types of environmental noise, such as white noise¹³, Lévy noise¹⁴ and colored noise¹⁵, each contribute unique advantages to disease modeling. Whit noise in epidemic modeling captures random, uncorrelated fluctuations that address uncertainties in disease spread. It is frequently incorporated into deterministic models to account for variability in transmission processes, human behavior. It offers a practical approach to introducing randomness, improving the realism of epidemic predictions. In recent studies, many researchers have utilized white noise in epidemic modeling to better understand and simulate stochastic effects in disease dynamics^16,17,18. A recent study on computer virus propagation using stochastic modeling approach has been presented in¹⁹.

The optimal control theory in epidemics plays a significant role in providing effective preventative control techniques²⁰. Many researchers have applied optimal control theory to various epidemiological models for effetely mitigating the outbreaks. For example in²¹, the authors studied a deterministic SEIQR type optimal control model for mitigating COVID-19 in Nigeria. A novel stochastic optimal controlling model of hepatitis C was analyzed in²². In⁵ fractional control of coronavirus epidemic in Algeria was investigated. In⁶, the researchers analyzed COVID-19 spread in Indonesia, incorporating vaccination effects using real data.

This study aims to formulate a stochastic and deterministic optimal control frameworks for modeling and managing the dynamics of the COVID-19 epidemic. Our approach integrates mathematical modeling and optimization techniques to develop strategies that minimize the spread of the virus while considering the constraints imposed by limited resources and societal needs. The present study initiate by formulating a COVID-19 compartmental model that incorporates the effects of perturbations induced by white noise. The resulting stochastic problem helps to provide a more realistic scenario of the inherent uncertainties and random fluctuations observed in various real-world problems including epidemic dynamics. We aim to gain a deeper insights into the dynamics of the pandemic and identify optimal control measures. The impact of various control measures, such as vaccination campaigns, social distancing, and treatment strategies on reducing the incidence and new outbreaks of the epidemics are analyzed graphically. To ensure the practical applicability of our model, we parameterize it using real-world data reported in Algeria for a specific time period. We estimate the model parameters via a statistical tool that best fit the observed data, thereby enhancing the accuracy and reliability of our mathematical framework. To analyze the effectiveness of our proposed control measures, we employ numerical simulations to depict the analytical results obtained from our model. The numerical results demonstrate the efficacy of different control strategies in managing the COVID-19 pandemic. This knowledge can inform policymakers and public health authorities in making informed decisions regarding the implementation of control measures, resource allocation, and the optimization of public health interventions.

This study is mainly divided in six sections. Mathematical formulation of the stochastic model is provided in “Mathematical formulation of the model”. Some qualitative aspects are presented in “Qualitative analysis of our stochastic model”. Formulation and analysis of optimal control problem is presented in “Formulation of stochastic optimal control”. Numerical simulation and parameter estimation are presented in “Numerical simulations” to support our results. Finally, the work is concluded in “Conclusion”.

Mathematical formulation of the model

Motivated by the above, as well as by works^6,7, we construct a stochastic mathematical model for the dynamics of the transmission of the COVID-19 epidemic. The classification of the totla population N(t) take place into four classes: S(t) susceptible, V(t) vaccinated, I(t) infected, R(t) recovered, at time $t>0$ with $N=S+V+I+R$.

According to the characteristics of COVID-19, the following presumptions are considered in our model: there are no negative parameters, the susceptible population gets vaccinated, infected people is able to transmit the infection to susceptible or vaccinated person, recovered persons have temporary immunity.

The assumption that white noise is proportional to the compartments is a widely used and practical approach for modeling random fluctuations in population dynamics^16,17,18,23. This approach captures the notion that larger population undergo greater random variations.

Considering that white noise is in direct proportional to S,  V,  I,  R. Thus, we get the following stochastic model:

$$\begin{aligned} \begin{aligned} dS(t)&=\Big [\Delta - \beta \dfrac{S(t)I(t)}{N} + \mu R(t) - (k+\delta )S(t)\Big ]dt+\rho _{1}S(t)dW_{1}(t)\\ dV(t)&= \Big [kS(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t)\Big ]dt+\rho _{2}V(t)dW_{2}(t)\\ dI(t)&= \Big [\beta \dfrac{S(t)I(t)}{N}+(1-\tau )\beta \dfrac{V(t)I(t)}{N}-(\alpha +\delta +\delta _{0} )I(t)\Big ]dt+\rho _{3}I(t)dW_{3}(t) \\ dR(t)&= \Big [\alpha I(t)-(\delta +\mu )R(t)\Big ]dt+\rho _{4}R(t)dW_{4}(t), \end{aligned} \end{aligned}$$

(1)

where $W_{j}$ $\{j=1,..,4\}$ are independently Brownian motion. $\rho _{j}$ $\{j=1,..,4\}$ denotes the thew white noise intensity.

If we put $\rho _{j}=0, \ j=1,2,3,4.$, the following deterministic version of model (1) will result:

$$\begin{aligned} \begin{aligned} dS(t)&=\Big [\Delta - \beta \dfrac{I(t)S(t)}{N} + \mu R(t) - (\delta +k )S(t)\Big ]dt, \\ dV(t)&= \Big [kS(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t)\Big ]dt,\\ dI(t)&= \Big [\beta \dfrac{I(t)S(t)}{N}+\beta (1-\tau )\dfrac{V(t)I(t)}{N}-(\alpha +\delta +\delta _{0} )I(t)\Big ]dt, \\ dR(t)&= \Big [\alpha I(t)-(\mu +\delta )R(t)\Big ]dt. \end{aligned} \end{aligned}$$

(2)

Detail of the model embedded parameters are listed in the below table 1.

Table 1 Parameter’s biological description.

The dynamical characteristics of the COVID-19 model (2) can be thoroughly analyzed based on the threshold number, which, as referenced in²⁴, and is defined as follows:

$$\begin{aligned} R_{0}^{d}=\dfrac{\beta \delta +(1-\tau )\beta k}{(k+\delta )(\alpha +\delta +\delta _{0})}. \end{aligned}$$

More precisely, if $R_{0}^{d} < 1$, the disease-free equilibrium $P_{0}=(\frac{\Delta }{k+\delta },\frac{k\Delta }{\delta (k+\delta )},0,0)$ is globally asymptotically stable (GAS). If $R_{0}^{d} > 1$ the model posses an endemic equilibrium denoted $P^{*} = (S^{*},V^{*}, I^{*}, R^{*})$, and it is further GAS.

In the subsequent section, we present a comprehensive qualitative study and dynamic analysis of compartmental model (1). In this study, we presume that $(\Omega ,{\mathcal {F}}, {\mathcal {P}})$ is a complete probability space coupled with a filtration $({\mathcal {F}}_{t})_{t\ge 0}$ that meets the standard condition.

Qualitative analysis of our stochastic model

Existence of unique and global positive solution

In the following, we will demonstrate that the system (1) has a unique positive solution.

Theorem 1

For all $(S(0),V(0),I(0),R(0)) \in {\mathbb {R}}^{4}_{+}$, a unique solution of the system (1) can be determined. Moreover, the solution will remain in ${\mathbb {R}}^{4}_{+}$$\forall$ $t\ge 0$ almost surely.

Proof

As the coefficients of the model are locally Lipschitz continues for all initial values in $\in {\mathbb {R}}^{4}_{+}$, there exists a unique solution (S(t), V(t), I(t), R(t)) for $t \in [0,\tau _{e})$, where $\tau _{e}$ denotes the final time. To show that the aforementioned solution is global, it is necessary to prove that $\tau _{e}=\infty$ a.s. Assuming that m is a sufficiently large non-negative number and S(0), V(0), I(0) and R(0) all lie within $[\frac{1}{m}, m]$. The stopping time is given by

$$\begin{aligned} \tau _{l}=\inf \biggl \{t\in [0, \tau _{e}): \min \Bigl \{S,V,I,R\Bigr \} \le \frac{1}{l} \ \text {Or} \ \max \Bigl \{ S,V,I,R\Bigr \} \ge l \biggr \}, \ \ \ \forall l \ge m, \end{aligned}$$

where in the subsequent results, we set $\inf \varnothing =\infty$ where, $\varnothing$ is an empty set. It is clear that $\tau _{l}$ has an increasing nature as $l\rightarrow \infty$. Further, setting $\tau _{\infty }=\lim _{l \rightarrow \infty }\tau _{l}$, with $\tau _{e} \ge \tau _{\infty } \ a.s$.

Under the confirmation that $\tau _{\infty }=\infty \ a.s.$, then, we have $\tau _{e}=\infty$ and in a result $(S(t),V(t),I(t),R(t))\in {\mathbb {R}}^{4}_{+} \ a.s$$\forall$ $t\ge 0$. Further, we proceed to prove that $\tau _{\infty }=\infty$. Based on the assumption that $\tau _{\infty }< \infty$, then there exists T greater than 0 and $\varepsilon \in (0,1)$ satisfying ${\mathbb {P}}\{\tau _{\infty }\le T\}\ge \varepsilon$. Therefore, there exists an integer $l_{1} \ge m$ such that

$$\begin{aligned} {\mathbb {P}}\{\tau _{l} \le T \} \ge \varepsilon , ~~ \forall ~~~ l \ge l_{1}. \end{aligned}$$

(3)

We define a $C^{2}$-function $\Sigma : {\mathbb {R}}^{4}_{+}$ $\rightarrow$ ${\mathbb {R}}_{+}$ as follows :

$$\begin{aligned} \Sigma (S,V,I,R)= S+V+I+R-4-(\ln S+ \ln V + \ln I+\ln R). \end{aligned}$$

$\Sigma$ is non-negativity as can verified via $x^*-\ln x^*-1 \ge 0$, $\forall x^*>0.$ Moreover, the implementation of Itô formula gives

$$\begin{aligned} d\Sigma (S,V,I,R)=&\biggl (1-\dfrac{1}{S}\biggr )dS +\dfrac{1}{2S^{2}}(dS)^{2}+\biggl (1-\dfrac{1}{V}\biggr )dV\\&+ \dfrac{1}{2V^{2}}(dV)^{2}+\biggl (1-\dfrac{1}{I}\biggr )dI +\dfrac{1}{2I^{2}}(dI)^{2}\\&+\biggl (1-\dfrac{1}{R}\biggr )dR +\dfrac{1}{2R^{2}}(dR)^{2} \\ =&\biggl (1-\dfrac{1}{S}\biggr )\biggl (\Big [\Delta - \beta \dfrac{S(t)I(t)}{N} + \mu R(t) - (k+\delta )S(t)\Big ]dt+\rho _{1}S(t)dW_{1}(t) \biggr ) \\&+\dfrac{1}{2S^{2}}\biggl ( \Big [\Delta - \beta \dfrac{S(t)I(t)}{N} + \mu R(t) - (k+\delta )S(t)\Big ]dt+\rho _{1}S(t)dW_{1}(t)\biggr )^{2}\\&+\biggl (1-\dfrac{1}{V}\biggr )\biggl (\Big [kS(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t)\Big ]dt+\rho _{2}V(t)dW_{2}(t) \biggr )\\ &+ \dfrac{1}{2V^{2}}\biggl ( \Big [kS(t)- \beta (1-\tau )\dfrac{V(t)I(t)}{N} -\delta V(t)\Big ]dt+\rho _{2}V(t)dW_{2}(t)\biggr )^{2} \\ &+\biggl (1-\dfrac{1}{I}\biggr )\biggl ( \Big [\beta \dfrac{S(t)I(t)}{N}+\beta (1-\tau )\dfrac{V(t)I(t)}{N}-(\alpha +\delta +\delta _{0} )I(t)\Big ]dt+\rho _{3}I(t)dW_{3}(t)\biggr )\\&+\dfrac{1}{2I^{2}}\biggl (\Big [\beta \dfrac{S(t)I(t)}{N}+\beta (1-\tau )\dfrac{V(t)I(t)}{N}-(\alpha +\delta +\delta _{0} )I(t)\Big ]dt+\rho _{3}I(t)dW_{3}(t) \biggr )^{2} \\ &+\biggl (1-\dfrac{1}{R}\biggr )\biggl (\Big [\alpha I(t)-(\delta +\mu )R(t)\Big ]dt+\rho _{4}R(t)dW_{4}(t) \biggr ) \\ &+\dfrac{1}{2R^{2}}\biggl ( \Big [\alpha I(t)-(\delta +\mu )R(t)\Big ]dt+\rho _{4}R(t)dW_{4}(t) \biggr )^{2} \\ =&L\Sigma (S,V,I,R) dt +\rho _{1} (S(t)-1)dW_{1}(t)+ \rho _{2} (V(t)-1)dW_{2}(t)+ \rho _{3} (I(t)-1) dW_{3}(t)\\&+ \rho _{4} (R(t)-1) dW_{4}(t), \end{aligned}$$

where $L\Sigma : \ {\mathbb {R}}^{4}_{+} \rightarrow {\mathbb {R}}_{+}$ is given by

$$\begin{aligned} L\Sigma (S,V,I,R)=&\Delta -\delta (S+V+R+I) -\dfrac{\Delta }{S}+\beta \dfrac{I}{N} - \mu \dfrac{R}{S}\\&- k\dfrac{S}{V}+\beta (1-\tau )\dfrac{I}{N}-\beta \dfrac{S}{N}-\delta _{0}I\\&-\beta (1-\tau )\frac{V}{N}-\alpha \dfrac{I}{R}+\mu +\alpha +4\delta +\delta _{0}+k+ \dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{3}^{2}+\rho _{4}^{2}}{2}\\ \le&\Delta +\beta +(1-\tau )\beta +4\mu +\alpha +\delta _{0} +k+\delta + \dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{3}^{2}+\rho _{4}^{2}}{2} : = M. \end{aligned}$$

In the given context, M denotes a constant that is independent of t, S, I, V, and R.

$$\begin{aligned} \Rightarrow \ \ \ d\Sigma (S,V,I,R)\le M dt + \rho _{1} (S-1)dW_{1}(t)+ \rho _{2} (V-1)dW_{2}(t)+ \rho _{3} (I-1) dW_{3}(t)+\rho _{4} (I-1) dW_{4}(t).\end{aligned}$$

Integration from 0 to $\tau _{l} \wedge T = \min \{ \tau _{l},T\}$ gives

$$\begin{aligned} \int _{0}^{\tau _{l} \wedge T} d\Sigma (S,V,I,R)\le & \int _{0}^{\tau _{l} \wedge T} M dr +\int _{0}^{\tau _{l} \wedge T} \rho _{1} (S(r)-1)dW_{1}(r)\\ & + \int _{0}^{\tau _{l} \wedge T} \rho _{2} (V(r)-1)dW_{2}(r) \\ & +\int _{0}^{\tau _{l} \wedge T} \rho _{3} (I(r)-1) dW_{3}(r) +\int _{0}^{\tau _{l} \wedge T} \rho _{4} (R(r)-1) dW_{4}(r).\end{aligned}$$

Taking expectation on both sides, We obtain

$$\begin{aligned} & {\mathbb {E}} \Big [\Sigma \big (S(\tau _{l} \wedge T ),V(\tau _{l} \wedge T ),I(\tau _{l} \wedge T ),R(\tau _{l} \wedge T )\big )\Big ]\\ & \quad \le \Sigma \big (S(0),V(0),I(0),R(0)\big )+ {\mathbb {E}} \bigg [ \int _{0}^{\tau _{l} \wedge T} M dr\bigg ] \\ & {\mathbb {E}} \Big [\Sigma \big (S(\tau _{l} \wedge T ),V(\tau _{l} \wedge T ),I(\tau _{l} \wedge T ),R(\tau _{l} \wedge T )\big )\Big ]\\ & \quad \le \Sigma \big (S(0),V(0),I(0),R(0)\big )+ MT. \end{aligned}$$

Let $\Gamma _{l}=\{\tau _{l}\le T\}$ $\forall \ l\ge l_{1}$ and by the result stated in (3) we have ${\mathbb {P}}(\Gamma _{l})\ge \varepsilon$. It is worth noticing that for every $\omega$ in $\Gamma _{l}$, we can find at least one $S(\tau _{l},\omega )$, $V(\tau _{l},\omega )$, $I(\tau _{l},\omega )$, and $R(\tau _{l} \wedge T)$ that are equal to l or $\frac{1}{l}$, consequently

$$\begin{aligned}\Sigma \big (S(\tau _{l},\omega ),V(\tau _{l},\omega ),I(\tau _{l},\omega ),R(\tau _{l},\omega )\big )\ge (l-\ln l-1) \wedge \big (\frac{1}{l}-1+\ln l \big ). \end{aligned}$$

Therefore,

$$\begin{aligned} \Sigma \big (S(0),V(0),I(0),R(0)\big )+ MT&\ge {\mathbb {E}}\Big [ {\mathbb {I}}_{\Gamma _{l}}(\omega ) \Sigma \big (S(\tau _{l} \wedge T ),V(\tau _{l} \wedge T ),I(\tau _{l} \wedge T ),R(\tau _{l} \wedge T )\big )\Big ] \\&= {\mathbb {E}}\Big [ {\mathbb {I}}_{\Gamma _{l}}(\omega ) \Sigma \big ( S(\tau _{l},\omega ),V(\tau _{l},\omega ),I(\tau _{l},\omega ),R(\tau _{l},\omega )\big )\Big ]\\&\ge {\mathbb {E}}\Big [ {\mathbb {I}}_{\Gamma _{l}}(\omega ) (l-\ln l-1) \wedge \big (\frac{1}{l}+\ln l-1 \big ) \Big ]\\&\ge {\mathbb {E}}\Big [ {\mathbb {I}}_{\Gamma _{l}}(\omega ) \Big ] (l-\ln l-1) \wedge \big (\frac{1}{l}+\ln l-1 \big ) \\&\ge \varepsilon (l-\ln l-1) \wedge \big (\frac{1}{l}+\ln l-1 \big ). \end{aligned}$$

${\mathbb {I}}_{\Gamma _{l}(\omega )}$ denotes the corresponding indicator function for $\Gamma _{l}$.

If $l \rightarrow \infty$, then $\infty > \Sigma \big (S(0),V(0),I(0),R(0)\big )+ MT = \infty$ is a contradiction. Hence, we conclude that $\tau _{\infty }=\infty$ a.s.. $\square$

Extinction criterions of the infection

In this part of the study, we focus on investigating the sufficient conditions that lead to infection extension in the system (1). This remains a fundamental problem in epidemiology.

Define a parameter

$$\begin{aligned} R_{0}=\dfrac{\beta +\beta (1-\tau )}{(\delta +\alpha +\delta _{0}+\frac{\rho _{3}^{2}}{2})}. \end{aligned}$$

Theorem 2

(S, V, I, R) be a solution of the model (1) with $(S(0), V(0), I (0), R(0))\in {\mathbb {R}}^{4}$, if$R_{0}<1$ then $\lim \limits _{t\rightarrow \infty }\sup \dfrac{\ln I(t)}{t}<0 \ a.s.$, namely, $I(t) \longrightarrow 0$ exponentially a.s.

Proof

Make use of Itô formula we have

$$\begin{aligned} d\ln I(t)=\Big [ \beta \dfrac{S(t)}{N}+(1-\tau )\beta \dfrac{V(t)}{N} -(\alpha +\delta +\delta _{0}+\frac{\rho _{3}^{2}}{2})\Big ]dt +\rho _{3}dW_{3}(t). \end{aligned}$$

(4)

Integrating over [0, t] and then with the division with t yields

$$\begin{aligned} \dfrac{\ln I(t)}{t}= & \dfrac{\beta }{t}\int _{0}^{t} \dfrac{S(r)}{N}dr +\dfrac{\beta (1-\tau )}{t}\int _{0}^{t} \dfrac{V(r)}{N} dr - (\alpha +\delta +\delta _{0}\nonumber \\ & +\frac{\rho _{3}^{2}}{2}) +\dfrac{\rho _{3}}{t}\int _{0}^{t}dW_{3}(r)+\dfrac{\ln I(0)}{t}\nonumber \\\le & \beta (1+(1-\tau )) - (\alpha +\delta +\delta _{0}+\frac{\varrho _{3}^{2}}{2}) +\dfrac{\rho _{3}}{t}\int _{0}^{t}dW_{3}(r)+\dfrac{\ln I(0)}{t}\nonumber \\\Rightarrow & \ \ \dfrac{\ln I(t)}{t}\le (\alpha +\delta +\delta _{0}+\frac{\rho _{3}^{2}}{2})(R_{0} - 1) +\dfrac{\rho _{3}}{t}\int _{0}^{t}dW_{3}(r)+\dfrac{\ln I(0)}{t}. \end{aligned}$$

(5)

Moreover $\frac{1}{t}\int _{0}^{t}dW_{3}(r)$ is a continuous local martingale. By lemma (Strong Laws of a Large Numbers), we obtain

$$\begin{aligned}\lim _{t\rightarrow \infty } \dfrac{1}{t}\int _{0}^{t}dW_{3}(r)=0 \ \ a.s.\end{aligned}$$

Taking $\lim \limits _{t\rightarrow \infty }\sup$ of (5) and if $R_{0}<1$ we get

$$\begin{aligned} \lim \limits _{t\rightarrow \infty }\sup \dfrac{\ln I(t)}{t}\le (\alpha +\delta +\delta _{0}+\frac{\rho _{3}^{2}}{2})(R_{0} - 1)<0, \end{aligned}$$

which further simplifies $\lim \limits _{t\rightarrow \infty } I(t)=0 \hspace{0.3cm} a.s.$ $\square$

Stationary distribution and ergodicity

Studying how diseases stay or die out in epidemiology is a crucial problem. In this part we will demonstrate that system described in (1) posses a stationary distribution with some specific conditions.

Let ${\mathcal {Z}}(t)$ follows the homogeneous Markov process in ${\mathbb {R}}^{n}$ which ca be shown by the SDE

$$\begin{aligned}d{\mathcal {Z}}(t)=b({\mathcal {Z}})dt+\sum _{r=1}^{\kappa }\sigma _{r}({\mathcal {Z}})dB_{r}(t).\end{aligned}$$

The diffusion matrix is descried as

$$\begin{aligned}{\mathcal {D}}(y)=\Big (d_{ij}(y)\Big ); \ \ \ d_{ij}(y)=\sum _{r=1}^{k}\sigma _{r}^{i}(y)\sigma _{r}^{j}(y).\end{aligned}$$

Lemma 1

²⁵ The Markov process ${\mathcal {Z}}(t)$ possesses a unique stationary distribution $\pi (\cdot )$ if there we can find a bounded domain $U \subset {\mathbb {R}}^{n}$ with a regular boundary $\Gamma$.

• (a) There exists a positive constant ${\mathcal {M}}$ such that

  $$\begin{aligned}\sum _{i,j=1}^{n} d_{ij}(y)\Theta _{i}\Theta _{j}\ge {\mathcal {M}}\vert \Theta \vert ^{2} \ y\in U, \ \Theta \in {\mathbb {R}}^{n}. \ \ \end{aligned}$$

• (b) A non-negative $C^{2}$-function V and a neighborhood U exists, where LV is negative for any ${\mathbb {R}}^{n}\setminus U$.

Theorem 3

Under the assumption

$$\begin{aligned}\bar{R}_{0}:= \dfrac{\delta k\beta (1-\tau )}{ (\delta +k+\frac{\rho _{1}^{2}}{2})(\delta +\frac{\rho _{2}^{2}}{2})(\delta _{0}+\alpha +\delta +\frac{\rho _{3}^{2}}{2})}>1,\end{aligned}$$

the solution of (1) posses a unique stationary distribution $\pi (.)$, and it is ergodic.

Proof

To prove Theorem 3, it is sufficient to verify the conditions of Lemma 1.

First we investigate a condition (b), We will constructing a non-negative $C^{2}$-function $\chi : {\mathbb {R}}^{4}_{+}\longrightarrow {\mathbb {R}}_{+}$.

We define

$$\begin{aligned}\Upsilon _{1}(S,V,I,R)=S+I+V+R-w_{1}\ln S-w_{2}\ln V - w_{3}\ln I,\end{aligned}$$

where $\omega _{1}, \omega _{2}, \omega _{3}$ are positive constants we will define later.

Utilizing the Itô formula we obtain

$$\begin{aligned} L(S+I+V+R)&=\Delta -\delta (S+I+V+R)-\delta _{0}I, \\ L(-\ln S)&=-\dfrac{\Delta }{S}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+\delta , \\ L(-\ln V)&= -\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}, \\ L(-\ln I)&= -\dfrac{\beta S}{N}-\dfrac{(1-\tau )\beta V}{N}+\alpha +\delta _{0}+\delta +\dfrac{\rho _{3}^{2}}{2}. \end{aligned}$$

Consequently,

$$\begin{aligned} L\Upsilon _{1}=&\Delta -\delta (S+V+I+R)-\delta _{0}I-\dfrac{w_{1}\Delta }{S}+\dfrac{w_{1}\beta I}{N} - \dfrac{w_{1}\mu R}{S}+w_{1}(\dfrac{\rho _{1}^{2}}{2}+k+\delta ) \\&-\dfrac{w_{2}k S}{V} +\dfrac{w_{2}(1-\tau )\beta I}{N}+w_{2}(\delta +\dfrac{\rho _{2}^{2}}{2}) -\dfrac{w_{3}\beta S}{N}-\dfrac{w_{3}(1-\tau )\beta V}{N}+w_{3}(\alpha +\delta _{0}+\delta +\dfrac{\rho _{3}^{2}}{2}), \end{aligned}$$

using the inequality $\root 4 \of {a_{1}a_{2}a_{3}a_{4}}\le \frac{1}{4}(a_{1}+a_{3}+a_{2}+a_{4}) \ \ \, a_{1}, a_{3} a_{2}, a_{4} >0$.

$$\begin{aligned} L\Upsilon _{1} \le&-4\root 4 \of {\delta (S+V+I+R)\frac{w_{1}\Delta }{S}\frac{w_{2}k S}{V}\frac{w_{3}(1-\tau )\beta V}{N}}+\Delta -\delta _{0}I+\dfrac{w_{1}\beta I}{N} - \dfrac{w_{1}\mu R}{S} \\&+w_{1}(\dfrac{\rho _{1}^{2}}{2}+k+\delta ) +\dfrac{w_{2}(1-\tau )\beta I}{N}+w_{2}(\delta +\dfrac{\rho _{2}^{2}}{2}) -\dfrac{w_{3}\beta S}{N}+w_{3}(\alpha +\delta _{0}+\delta +\dfrac{\rho _{3}^{2}}{2})\\ =&-4\root 4 \of {\delta w_{1}\Delta w_{2} k w_{3}(1-\tau )\beta }+ \Delta -\delta _{0}I+\dfrac{w_{1}\beta I}{N} - \dfrac{w_{1}\mu R}{S}+w_{1}(\dfrac{\rho _{1}^{2}}{2}+k+\delta ) +\dfrac{w_{2}(1-\tau )\beta I}{N}\\&+w_{2}(\delta +\dfrac{\rho _{2}^{2}}{2}) -\dfrac{w_{3}\beta S}{N}+w_{3}(\alpha +\delta _{0}+\delta +\dfrac{\rho _{3}^{2}}{2}). \end{aligned}$$

Putting

$$\begin{aligned} \Delta =w_{1}(\delta +k+\frac{\rho _{1}^{2}}{2})=w_{2}(\delta +\frac{\rho _{2}^{2}}{2})=w_{3}(\delta _{0}+\alpha +\delta +\frac{\rho _{3}^{2}}{2}), \end{aligned}$$

then

$$\begin{aligned}w_{1}= \dfrac{\Delta }{(\delta +k+\frac{\rho _{1}^{2}}{2})}, \ \ \ \ w_{2}=\dfrac{\Delta }{(\delta +\frac{\rho _{2}^{2}}{2})},\ \ \ \ w_{3}=\dfrac{\Delta }{(\delta _{0}+\alpha +\delta +\frac{\rho _{3}^{2}}{2})}.\end{aligned}$$

As a result

$$\begin{aligned} L\Upsilon _{1}&\le -4\root 4 \of {\dfrac{\delta k\beta (1-\tau )\Delta ^{4}}{ (\delta +k+\frac{\rho _{1}^{2}}{2})(\delta +\frac{\rho _{2}^{2}}{2})(\delta _{0}+\alpha +\delta +\frac{\rho _{3}^{2}}{2})}}+4\Delta -\delta _{0}I+\dfrac{w_{1}\beta I}{N} - \dfrac{w_{1}\mu R}{S}\\&\quad +\dfrac{w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{3}\beta S}{N} \\&\le -4\Delta \Big (\root 4 \of {\dfrac{\delta k\beta (1-\tau )}{ (\delta +k+\frac{\rho _{1}^{2}}{2})(\delta +\frac{\rho _{2}^{2}}{2})(\delta _{0}+\alpha +\delta +\frac{\rho _{3}^{2}}{2})}}-1\Big )+\dfrac{w_{1}\beta I}{N} +\dfrac{w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{3}\beta S}{N}\\&=-4\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{1}\beta I}{N} +\dfrac{w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{3}\beta S}{N}. \end{aligned}$$

We define

$$\begin{aligned} \Upsilon _{2}(S,V,I,R)=w_{4}\big (S+V+I+R-w_{1}\ln S-w_{2}\ln V - w_{3}\ln I\big )+S+V+I+R-\ln S-\ln R -\ln V.\end{aligned}$$

Obviously

$$\begin{aligned} \liminf \limits _{q\rightarrow +\infty ,(S,V,I,R)\in {\mathbb {R}}_{+}^{4}\setminus U_{q}}\Upsilon _{2}(S,V,I,R)= +\infty , \end{aligned}$$

where $U_{q}=]\frac{1}{q},q[\times ]\frac{1}{q},q[\times ]\frac{1}{q},q[\times ]\frac{1}{q},q[$.

On the other hand

$$\begin{aligned} & \dfrac{\partial \Upsilon _{2}(S,V,I,R)}{\partial S}=w_{4}+1-\dfrac{w_{4}w_{1}+1}{S}, \ \ \ \ \ \hspace{0.5cm} \dfrac{\partial \Upsilon _{2}(S,V,I,R)}{\partial V}=w_{4}+1-\dfrac{w_{4}w_{2}+1}{V},\\ & \dfrac{\partial \Upsilon _{2}(S,V,I,R)}{\partial I}= w_{4}+1-\dfrac{w_{4}w_{3}}{I}\ \ \ \ \ \hspace{0.5cm} \dfrac{\partial \Upsilon _{2}(S,V,I,R)}{\partial R}= w_{4}+1-\dfrac{1}{R},\end{aligned}$$

the function $\Upsilon _{2}$ has a unique stagnation point that is

$$\begin{aligned} \big (S_{*},V_{*},I_{*},R_{*}\big )=\Big ( \dfrac{w_{4}w_{1}+1}{w_{4}+1}, \dfrac{w_{4}w_{2}+1}{1+w_{4}}, \dfrac{w_{4}w_{3}}{1+w_{4}},\dfrac{1}{1+w_{4}} \Big ).\end{aligned}$$

The Hessian matrix of $\Upsilon _{2}$ at $(S_{*},V_{*},I_{*},R_{*})$ is given by

$$\begin{aligned} \begin{bmatrix} \frac{w_{4}w_{1}+1}{S_{*}^{2}}& 0& 0 \\ 0& \frac{w_{4}w_{2}+1}{V_{*}^{2}}& 0 \\ 0& 0& \frac{w_{4}w_{3}}{I_{*}^{2}} \\ 0& 0& 0& \frac{1}{R_{*}^{2}} \end{bmatrix} \end{aligned}$$

The Hessian matrix is positive definite and $\Upsilon _{2}$ is continuous function then we can conclude that $\Upsilon _{2}$ has a unique minimum $(S_{*},V_{*},I_{*},R_{*}) \in {\mathbb {R}}^{4}_{+}$.

Let

$$\begin{aligned}\chi (S,V,I,R)=\Upsilon _{2}(S,V,I,R)-\Upsilon _{2}(S_{*},V_{*},I_{*},R_{*}).\end{aligned}$$

Applying Itô formula we obtain

$$\begin{aligned} L\chi (S,V,I,R)&=L\Upsilon _{2}(S,V,I,R)\\&\le w_{4}\Big [-4\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{1}\beta I}{N} +\dfrac{w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{3}\beta S}{N}\Big ]+\Delta -\delta (S+V+I+R) \\&\quad -\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\&= -4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&\quad -\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}. \end{aligned}$$

Let $\varepsilon _{l}>0, l=1,2,...,8$, we state a bounded closed set

$$\begin{aligned}U=\Bigl \{(S,V,I,R)\in {\mathbb {R}}_{+}^{4}: \varepsilon _{1}\le S\le \frac{1}{\varepsilon _{2}}, \varepsilon _{3}\le V\le \frac{1}{\varepsilon _{4}}, \varepsilon _{5}\le I\le \frac{1}{\varepsilon _{6}},\varepsilon _{7}\le R\le \frac{1}{\varepsilon _{8}}\Bigr \}.\end{aligned}$$

We can choose $\varepsilon _{l}$ and $w_{4}$ positive constant such that the subsequent conditions fulfilled

• $w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +\dfrac{\rho _{1}^{2}}{2}+k+3\delta +\mu +\dfrac{\rho _{4}^{2}}{2}+\dfrac{\rho _{2}^{2}}{2}-\dfrac{\Delta }{\varepsilon _{1}}<0,$

• $w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-\dfrac{k}{\varepsilon _{1}}<0,$

• $w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\delta \varepsilon _{1}+\dfrac{\beta }{\varepsilon _{1}^{2}}<0,$

• $w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\alpha }{\varepsilon _{5}}<0,$

• $w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\delta }{\varepsilon _{2}} <0$,

• $w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\delta }{\varepsilon _{4}}<0,$

• $w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-\dfrac{\delta }{\varepsilon _{6}} <0,$

• $w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\delta }{\varepsilon _{8}} <0.$

Next step we show that $L\chi <0$ on ${\mathbb {R}}^{4}_{+}\setminus U$, ${\mathbb {R}}^{4}_{+}\setminus U=\bigcup _{j=1}^{8}U_{j}$ where

$$\begin{aligned} & U_{1}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},0<S<\varepsilon _{1}\Bigr \}, \ \ \ \ U_{2}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},0<V<\varepsilon _{3}, S>\varepsilon _{1}\Bigr \}, \\ & U_{3}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},0<I<\varepsilon _{5},S>\varepsilon _{1}\Bigr \},\ \ \ \ U_{4}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},0<R<\varepsilon _{7}, I>\varepsilon _{5}\Bigr \}, \\ & U_{5}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},S>\dfrac{1}{\varepsilon _{2}}\Bigr \},\ \ \ \ U_{6}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},V>\dfrac{1}{\varepsilon _{4}}\Bigr \}, \\ & U_{7}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},I>\dfrac{1}{\varepsilon _{6}}\Bigr \}, \ \ \ \ U_{8}=\Bigl \{(S,V,I,R)\in {\mathbb {R}}^{4}_{+},R>\dfrac{1}{\varepsilon _{8}}\Bigr \}. \end{aligned}$$

Case 1. For $(S,V,I,R)\in U_{1}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&{w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +\dfrac{\rho _{1}^{2}}{2}+k+3\delta +\mu +\dfrac{\rho _{4}^{2}}{2}+\dfrac{\rho _{2}^{2}}{2}-\dfrac{\Delta }{S} } \\ \le&{w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +\dfrac{\rho _{1}^{2}}{2}+k+3\delta +\mu +\dfrac{\rho _{4}^{2}}{2}+\dfrac{\rho _{2}^{2}}{2}-\dfrac{\Delta }{\varepsilon _{1}}<0.} \end{aligned}$$

Case 2. For $(S,V,I,R)\in U_{2}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-k\dfrac{S}{V} \\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\beta +\Delta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-k\dfrac{\varepsilon _{1}}{\varepsilon _{3}}.\\ \text {Putting} \ \ \ \varepsilon _{3}=\varepsilon _{1}^{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\beta +\Delta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-\dfrac{k}{\varepsilon _{1}}<0. \end{aligned}$$

Case 3. Let $(S,V,I,R)\in U_{3}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +(w_{4}w_{2}+1)(1-\tau )\beta +\Delta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-\delta S+\beta \dfrac{I}{N}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +k+\Delta +3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-\delta \varepsilon _{1}+\beta \dfrac{\varepsilon _{5}}{\varepsilon _{1}}.\\ \text {Let } \varepsilon _{5}=\frac{1}{\varepsilon _{1}^{2}}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\delta \varepsilon _{1}+\dfrac{\beta }{\varepsilon _{1}^{2}}<0. \end{aligned}$$

Case 4. Considering $(S,V,I,R)\in U_{4}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\alpha \dfrac{I}{R}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\alpha \dfrac{\varepsilon _{7}}{\varepsilon _{5}}\\ \text {let } \varepsilon _{7}=\varepsilon _{5}^{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}(1-\tau )\beta +\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\alpha }{\varepsilon _{5}}<0. \end{aligned}$$

Case 5. For $(S,V,I,R)\in U_{5}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\delta S\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\delta }{\varepsilon _{2}} <0. \end{aligned}$$

Case 6. For $(S,V,I,R)\in U_{6}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\delta V\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\delta }{\varepsilon _{4}}<0. \end{aligned}$$

Case 7. For $(S,V,I,R)\in U_{7}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\delta I\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{4}^{2}+\rho _{2}^{2}+\rho _{1}^{2}}{2}-\dfrac{\delta }{\varepsilon _{6}} <0. \end{aligned}$$

Case 8. For $(S,V,I,R)\in U_{8}$

$$\begin{aligned} L\chi (S,V,I,R)\le&-4w_{4}\Delta \big [(R^{s}_{0})^{\frac{1}{4}}-1\big ] +\dfrac{w_{4}w_{1}\beta I}{N} +\dfrac{w_{4}w_{2}(1-\tau )\beta I}{N} -\dfrac{w_{4}w_{3}\beta S}{N}+\Delta -\delta (S+V+I+R)\\&-\dfrac{\Delta }{S}-\dfrac{\alpha I}{R}+\dfrac{\beta I}{N} - \dfrac{\mu R}{S}+\dfrac{\rho _{1}^{2}}{2}+k+2\delta +\mu +\dfrac{\rho _{4}^{2}}{2}-\dfrac{k S}{V} +\dfrac{(1-\tau )\beta I}{N}+\delta +\dfrac{\rho _{2}^{2}}{2}\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\delta R\\ \le&w_{4}w_{1}\beta +w_{4}w_{2}\beta (1-\tau )+\Delta +\beta +k+3\delta +\mu +\dfrac{\rho _{1}^{2}+\rho _{2}^{2}+\rho _{4}^{2}}{2}-\dfrac{\delta }{\varepsilon _{8}} <0. \end{aligned}$$

According to the above $L\chi <0$ for all $(S,V,I,R)\in {\mathbb {R}}^{4}_{+}\setminus U$ therefor, the condition (b) of 1 satisfied.

Next, we show a condition (a). The diffusion matrix of the problem (1) is determined by

$$\begin{aligned} \begin{bmatrix} \rho _{1}^{2}S^{2}& 0& 0& 0\\ 0& \rho _{2}^{2}V^{2}& 0& 0 \\ 0& 0& \rho _{3}^{2}I^{2}& 0\\ 0& 0& 0& \rho _{4}^{2}R^{2}\\ \end{bmatrix}, \end{aligned}$$

and

$$\begin{aligned} \sum _{i,j=1}^{4}a_{ij}(S,V,I,R)\Theta _{i}\Theta _{j}&=\rho _{1}^{2}S^{2} \Theta _{1}^{2}+\rho _{2}^{2}V^{2}\Theta _{2}^{2}+\rho _{3}^{2}I^{2}\Theta _{3}^{2}+\rho _{4}^{2}R^{2}\Theta _{4}^{2} \\&\ge m \vert \Theta \vert ^{2} \ \ \text {for all} \ \ \ (S,V,I,R) \in U, \ \Theta \in {\mathbb {R}}^{4}_{+}, \end{aligned}$$

where $m=\min \{\rho _{1}^{2}S^{2}, \rho _{2}^{2}V^{2}, \rho _{3}^{2}I^{2}, \rho _{4}^{2}R^{2}\}$, thus, the condition (a) is verified. Therefore, the transmission model (1) admits a uniquely determined ergodic stationary distribution, denoted as $\pi (\cdot )$. $\square$

Formulation of stochastic optimal control

This section aims to determine the optimally intervention measures for minimizing the transmission of illness via a mathematical optimization techniques. The main goal is to identify the controls that, while taking into account a variety of constraints, reduce the overall impact of the epidemic. In order to do this, the epidemic model is constructed and then control variables are included to reflect the various interventions that can be employed to stop the disease’s spread. Measures such as social isolation, immunization, quarantine and treatment are examples of these control variables.

The following two control variables $c_{1}(t),$ and $c_{2}(t)$ will be applied in this work to analyze our suggested models (1)and (2), where

• The control variable $c_{1}(t)$ is the COVID-19 time-dependent vaccine.

• The control variable $c_{2}(t)$ is the time-dependent treatment effects.

The objective in this study is to reduce susceptibility and infection burden and to maximize the recovered individuals at the lowest cost of control variants.

Analysis of deterministic optimal control problem

In this part, we introduced a control scheme that is considered to be optimal if it reduce the cost function

$$\begin{aligned} J_{1}(c_{1},c_{2})=\int _{0}^{T_{e}}\Big [a_{1}S+a_{2}I+\frac{1}{2}(p_{1}c^{2}_{1}(t)+p_{2}c^{2}_{2}(t))\Big ]dt, \end{aligned}$$

(6)

subject to

$$\begin{aligned} \begin{aligned} dS(t)&=\Big [\Delta - \beta \dfrac{I(t)S(t)}{N} + \mu R(t) - (kc_{1}(t)+\delta )S(t)\Big ]dt, \\ dV(t)&= \Big [kS(t)c_{1}(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t)\Big ]dt,\\ dI(t)&= \Big [\beta \dfrac{I(t)S(t)}{N}+\beta (1-\tau )\dfrac{I(t)V(t)}{N}-(\alpha c_{2} +\delta +\delta _{0} )I(t)\Big ]dt, \\ dR(t)&= \Big [\alpha c_{2}(t) I(t)-(\delta +\mu )R(t)\Big ]dt, \end{aligned} \end{aligned}$$

(7)

along with the following initial values of the state variables

$$\begin{aligned} S(0)=S_{0}>0, I(0)=I_{0}\ge 0, V(0)=V_{0}\ge 0, R(0)=R_{0}>0. \end{aligned}$$

(8)

In cost function (6), the positive constants $a_{1},a_{2}$ represent the balancing factors to maintain a balance of the size of susceptible, infected individuals respectively. The positive coefficients $p_{1}$ and $p_{2}$ are the cost on vaccination and treatment respectively and $T_{e}$ is the end time. From the expression of the cost function (6), our main goal is to reduce the number of the susceptible and the infected while reducing the cost of controls $c_{1}$ and $c_{2}$. Hence, we are intendment to evaluate optimal variables $(c^{*}_{1},c^{*}_{2})\in \Omega$ such that

$$\begin{aligned} J_{1}(c^{*}_{1},c^{*}_{2})=\min _{\Omega }J_{1}(c_{1},c_{2}), \end{aligned}$$

(9)

subject to the system (7) with (8), where $\Omega$ is control admissible set

$$\begin{aligned}\Omega :=\Bigl \{(c_{1},c_{2}) \vert \ \ c_{j}:[0,T_{e}]\rightarrow [0,1] \ \text {Lebesgue measurable},\ j=1,2.\Bigr \}.\end{aligned}$$

Existence of optimal control problem

This section present the necessary condition for the existence of a solution of problem (6)-(8), and then we prove the existence of optimal control variables. The control variables $c_{1}$ and $c_{2}$ are positive, bounded, and Lebesgue measurable. With a positive initial value, the solution remains positive and bounded.

Let

$$\begin{aligned}Z_{t}=AZ+G(Z)\end{aligned}$$

, where

$$\begin{aligned} Z= & \begin{bmatrix} S(t)\\ V(t)\\ I(t)\\ R(t),\\ \end{bmatrix}, \hspace{0.25 cm} A= \begin{bmatrix} -(kc_{1}(t)+\delta )& 0& 0& \mu \\ kc_{1}(t)& -\delta & 0& 0\\ 0& 0& -(\alpha c_{2}(t)+\delta +\delta _{0})& 0\\ 0& 0& \alpha c_{2}(t) & -(\delta +\mu )\\ \end{bmatrix}, \\ G(Z)= & \begin{bmatrix} \Delta - \frac{\beta I(t)S(t)}{N}\\ -\frac{\beta (1-\tau )I(t)V(t)}{N} \\ \frac{\beta I(t)S(t)}{N}+\frac{\beta (1-\tau )I(t)V(t)}{N}\\ 0 \end{bmatrix}. \end{aligned}$$

$$\begin{aligned} G(Z_{1})-G(Z_{2})&\le b_{1}\vert S_{1}-S_{2}\vert +b_{2}\vert V_{1}-V_{2}\vert +b_{3}\vert I_{1}-I_{2}\vert + b_{4}\vert R_{1}-R_{2}\vert \\&\le b\bigg [\vert S_{1}-S_{2}\vert +\vert V_{1}-V_{2}\vert +\vert I_{1}-I_{2}\vert +\vert R_{1}-R_{2}\vert \bigg ], \end{aligned}$$

where $b = \max \{b_{1}, b_{2}, b_{3}, b_{4}\}$ denotes a positive constant that remains independent of the model’s state variables.

$$\begin{aligned} \vert F(Z_{1})-F(Z_{2})\vert \le 2M\vert Z_{1}-Z_{2}\vert ,\end{aligned}$$

and $M=\max \{b,\Vert A\Vert \}<\infty$n, then F is uniformly Lipshitz continuous. The control and state parameters are clearly positive hence, we conclude that the solution of model (7) exist.

Now, we show the result that ensures that we have the optimal control that minimizes the cost function

Theorem 4

There is a control variable $c^{*}=(c^{*}_{1},c^{*}_{2})$ $\in \Omega$ satisfying

$$\begin{aligned} J_{1}(c^{*}_{1},c^{*}_{2})=\min _{\Omega }J_{1}(c_{1},c_{2}),\end{aligned}$$

for all $(S(0),V(0),I(0),R(0)) \in {\mathbb {R}}^{4}_{+}$.

Proof

We need to verify the conditions of Corollary (4.1) in²⁶ for the necessary result

• The set of admissible control with the respective state variables is nonempty, because the solution of system exist.

• The set $\Omega$ is closed and convex and therefore, by definition $\Omega$ is closed. Let $u_{1}, u_{2} \in \Omega$ and $0\le \lambda \le 1$ such that $u_{1}=(c_{1},c_{2})$ and $u_{2}=(c_{1}^{\prime },c_{2}^{\prime })$. Hence, $\lambda u_{1}+(1-\lambda )u_{2}=(\lambda c_{1}+(1-\lambda )c_{1}^{\prime },\lambda c_{2}+(1-\lambda )c_{2}^{\prime })$ $\in \Omega$. Thus $\Omega$ is convex.

• The optimal control system is bounded by a linear function in the state and control variables. Our system is linear in $c_{1}$ and $c_{2}$, in addition, the solution is absolutely continuous and from it we conclude the boundness of the solution, hence, the condition is fulfilled.

• The integrand ${\mathcal {L}}_{1}$ of the cost function is strictly convex on $\Omega$, because the Hessian matrix of ${\mathcal {L}}_{1}$ on $\Omega$ is given by

  $$\begin{aligned}\begin{bmatrix} p_{1}& 0\\ 0& p_{2} \\ \end{bmatrix} \end{aligned}$$

  is definite positive.

• There exist $m_{1}>0,m_{2}>0$ and $m_{3}>1$ such that

  $$\begin{aligned} {\mathcal {L}}_{1}(c_{1},c_{2})\ge m_{1}\vert (c_{1},c_{2})\vert ^{m_{3}}-m_{2}, \end{aligned}$$
  $$\begin{aligned} {\mathcal {L}}_{1}(c_{1},c_{2})&= a_{1}S+a_{2}I+\frac{1}{2}(p_{1}c^{2}_{1}(t)+p_{2}c^{2}_{2}(t)\\ &\ge \frac{1}{2}min\{p_{1},p_{2}\} (c^{2}_{1}(t)+c^{2}_{2}(t))\\&\ge \frac{1}{2}min\{p_{1},p_{2}\} (c^{2}_{1}(t)+c^{2}_{2}(t))-1\\&=\frac{1}{2}min\{p_{1},p_{2}\} \vert (c_{1},c_{2})\vert ^{2}-1. \end{aligned}$$

  Let $m_{1}=\frac{1}{2}min\{p_{1},p_{2}\}$ , $m_{2}=1$ and $m_{3}=2$.

Thus, the above conditions guarantee the existence of the optimal set $(c^{*}_{1},c^{*}_{2})$. $\square$

Solution of the optimal control system

This section investigate the necessary optimality criterion and the characterization of control optimal for our control problem, by using Pontryagin’s principle²⁶. This principle shifts the problem (6)–(8) into a problem of minimizing point-wise a Hamiltonian ${\mathcal {H}}_{1}$ corresponds to the controls $c_{1},c_{2}$. We define the Lagrangian ${\mathcal {L}}_{1}$ associated to the control problem (7)

$$\begin{aligned}{\mathcal {L}}_{1}(S,I,c_{1},c_{2})=a_{1}S+a_{2}I+\frac{1}{2}(b_{1}c^{2}_{1}(t)+b_{2}c^{2}_{2}(t)),\end{aligned}$$

and the Hamiltonian function ${\mathcal {H}}_{1}$ is given by

$$\begin{aligned}{\mathcal {H}}_{1}(z,c,\gamma )=\gamma .\varphi (z,c)+{\mathcal {L}}_{1}(z,c),\end{aligned}$$

where

$$\begin{aligned}z=(S,V,I,R), \ \gamma =(\gamma _{1}, \ \gamma _{2},\gamma _{3},\gamma _{4}),\ c=(c_{1},c_{2}), \ \varphi =(\varphi _{1},\varphi _{2},\varphi _{3},\varphi _{4}).\end{aligned}$$

$$\begin{aligned} \varphi _{1}(z,c)&=\Delta - \beta \dfrac{S(t)I(t)}{N} + \mu R(t) - (kc_{1}(t)+\delta )S(t),\\ \varphi _{2}(z,c)&=kc_{1}(t)S(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t),\\ \varphi _{3}(z,c)&=\dfrac{S(t)I(t)}{N}+(1-\tau )\beta \dfrac{V(t)I(t)}{N}-(\alpha c_{2} +\delta +\delta _{0} )I(t),\\ \varphi _{4}(z,c)&=\alpha c_{2}(t) I(t)-(\delta +\mu )R(t). \end{aligned}$$

So the Hamiltonian will have this forme

$$\begin{aligned} {\mathcal {H}}_{1}(z,c,\gamma )=&a_{1}S+a_{2}I+\frac{1}{2}(b_{1}c^{2}_{1}(t)+b_{2}c^{2}_{2}(t)) \\&+\gamma _{1}\big [\Delta - \beta \dfrac{S(t)I(t)}{N} + \mu R(t) - (kc_{1}(t)+\delta )S(t)\big ]\\&+ \gamma _{2}\big [kc_{1}(t)S(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t) \big ] \\&+ \gamma _{3}\big [\dfrac{S(t)I(t)}{N}+(1-\tau )\beta \dfrac{V(t)I(t)}{N}-(\alpha c_{2} +\delta +\delta _{0} )I(t) \big ] \\&+ \gamma _{4}\big [ \alpha c_{2}(t) I(t)-(\delta +\mu )R(t)\big ]. \end{aligned}$$

To find the characterization of $(c^{*}_{1}(t),c^{*}_{2}(t))$, we apply Pontryagin minimum principle which states that if $(z^{\star },c^{\star })$ be an optimality solution for the control system, we can determine a non-trivial vector say $\gamma$ such that

$$\begin{aligned} \begin{aligned} \dfrac{d\gamma (t)}{dt}&=-\dfrac{\partial {\mathcal {H}}_{1}}{\partial z}(z^{\star }(t),c^{\star }(t),\gamma (t))\\ 0&= \dfrac{\partial {\mathcal {H}}_{1}}{dc}(z^{\star }(t),c^{\star }(t),\gamma (t)), \end{aligned} \end{aligned}$$

(10)

with transversality condition

$$\begin{aligned} \gamma (T_{e})=0, \end{aligned}$$

(11)

and the maximality condition

$$\begin{aligned} {\mathcal {H}}_{1}(z^{\star }(t),c^{\star }(t),\gamma (t))=\min _{c \in \Omega }{\mathcal {H}}_{1}(z^{\star }(t),c(t),\gamma (t)). \end{aligned}$$

(12)

Theorem 5

Let $(S^{\star }V^{\star }I^{\star }R^{\star })$ be the optimal solution associated to the optimal control $c^{\star }_{1}$ and $c^{\star }_{2}$ of the proposed problem, then $\exists$ adjoint function $\gamma _{j}(t)\ne 0, \ \text {for} \ j=1...4$ that satisfies the following

$$\begin{aligned} \begin{aligned} \gamma ^{\prime }_{1}&=\beta \dfrac{I}{N}\gamma _{1}+kc_{1}\gamma _{1}-a_{1},\\ \gamma ^{\prime }_{2}&=(1-\tau )\beta \dfrac{I}{N}\gamma _{2}+\delta \gamma _{2},\\ \gamma ^{\prime }_{3}&=-\beta \dfrac{S}{N}\gamma _{3}+(\tau -1)\beta \dfrac{V}{N}\gamma _{3}+(\alpha c_{2}+\delta +\delta _{0})\gamma _{3}-a_{2},\\ \gamma ^{\prime }_{4}&=(\delta +\mu )\gamma _{4}, \end{aligned} \end{aligned}$$

(13)

with terminal conditions

$$\begin{aligned} \gamma _{j}(T_{e})=0, \ \ j=1,...4, \end{aligned}$$

(14)

and

$$\begin{aligned} \begin{aligned} c^{*}_{1}&=\min \Bigl \{\max \Bigl \{\dfrac{kS(\gamma _{1}-\gamma _{2})}{a_{3}},0\Bigr \},1\Bigr \},\\ c^{*}_{2}&=\min \Bigl \{\max \Bigl \{\dfrac{\alpha I(\gamma _{3}-\gamma _{4})}{a_{4}},0\Bigr \},1\Bigr \}. \end{aligned} \end{aligned}$$

(15)

Proof

Applying the adjoint system (10) of the Pontryagin Minimum Principle, it gives us the system (13), and the condition (14) is a direct result of terminal condition (11).

In order to get the formula of the optimal control $c^{\star }_{1},c^{\star }_{2}$ which is given by (15), we solve the following system

$$\begin{aligned}&\dfrac{d{\mathcal {H}}_{1}}{dc_{1}}=kS(\gamma _{2}-\gamma _{1})+a_{3}c_{1}=0,\\&\dfrac{d{\mathcal {H}}_{1}}{dc_{2}}=\alpha I(\gamma _{4}-\gamma _{3})+a_{4}c_{2} =0, \end{aligned}$$

by considering the upper and lower limits of the control. Thus, we obtain the result (15) $\square$

After that, we will discuss the control analysis of stochastic case of problem (2) with the application of the same measures mentioned above.

Stochastic optimal control strategy

In this part, we formulate and study the stochastic control version of model (1) by considering the same control function. The proposed stochastic control system as following

$$\begin{aligned} \begin{aligned} dS(t)&=\Big [\Delta - \beta \dfrac{I(t)S(t) }{N} + \mu R(t) - (kc_{1}(t)+\delta )S(t)\Big ]dt+\rho _{1}S(t)dW_{1}(t), \\ dV(t)&= \Big [kc_{1}(t)S(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t)\Big ]dt+\rho _{2}V(t)dW_{2}(t),\\ dI(t)&= \Big [\beta \dfrac{I(t)S(t)}{N}+(1-\tau )\beta \dfrac{V(t)I(t)}{N}-(\alpha c_{2} +\delta +\delta _{0} )I(t)\Big ]dt+\rho _{3}I(t)dW_{3}(t), \\ dR(t)&= \Big [\alpha c_{2}(t) I(t)-(\mu +\delta )R(t)\Big ]dt+\rho _{4}R(t)dW_{4}(t). \end{aligned} \end{aligned}$$

(16)

To simplify writhen, we put

$$\begin{aligned}z=(S,V,I,R), \ \ c=(c_{1},c_{2}), \ \ h= (h_{1},h_{2},h_{3},h_{4}), \ \ \sigma =(\sigma _{1},\sigma _{2},\sigma _{3},\sigma _{4}),\end{aligned}$$

then system (16) become

$$\begin{aligned}dz(t)=h(z,c)dt+\sigma (z)dW(t),\end{aligned}$$

where

$$\begin{aligned} h_{1}(z(t),c(t))&=\Delta - \beta \dfrac{I(t)S(t)}{N} + \mu R(t) - (kc_{1}(t)+\delta )S(t), \\ h_{2}(z(t),c(t))&=kc_{1}(t)S(t)- (1-\tau )\beta \dfrac{V(t)I(t)}{N} -\delta V(t), \\ h_{3}(z(t),c(t))&=\beta \dfrac{S(t)I(t)}{N}+(1-\tau )\beta \dfrac{V(t)I(t)}{N}-(\alpha c_{2} +\delta +\delta _{0} )I(t), \\ h_{4}(z(t),c(t))&= \alpha c_{2}(t) I(t)-(\mu +\delta )R(t), \end{aligned}$$

and

$$\begin{aligned}\sigma _{1}=\rho _{1}S, \ \ \sigma _{2}=\rho _{2}V,\ \ \sigma _{3}=\rho _{3}I, \ \ \sigma _{4}=\rho _{4}R. \end{aligned}$$

The Cost function is given by

$$\begin{aligned}j_{2}(c)={\mathbb {E}}\Big [\int _{0}^{T_{e}}\Big [a_{1}S+a_{2}I+\frac{1}{2}(a_{3}c^{2}_{1}(t)+a_{4}c^{2}_{2}(t))\Big ]dt, +\frac{p_{1}}{2}S^{2}(T_{e})+\frac{p_{2}}{2}V^{2}(T_{e})+\frac{p_{3}}{2}I^{2}(T_{e})+\frac{p_{4}}{2}R^{2}(T_{e})\Big ],\end{aligned}$$

where $a_{i},p_{i}\ \ i=1,...,4$ are non-negative real numbers.

Our objective to find optimal controls $c^{\star }=(c^{\star }_{1},c^{\star }_{2})$, such that

$$\begin{aligned}j_{2}(c^{\star })=\min _{c \in \Omega ^{*}}j_{2}(c), \end{aligned}$$

where

$$\begin{aligned}\Omega ^{*}=\bigl \{c_{j}(t): \ \ c_{j}(t)\in [0,c^{max}_{j}],\ \forall c_{j}\in L^{2}[0,T_{e}], \ \ j=1,2\bigr \}.\end{aligned}$$

We define the Hamiltonian function

$$\begin{aligned}{\mathcal {H}}_{2}(z,c,v,\eta )=\langle h(z,c),v\rangle +\langle \sigma (z),\eta \rangle +{\mathcal {L}}_{2}(z,c), \end{aligned}$$

where $\langle .,.\rangle$ represent the inner product in Euclidean space and $v = (v_{1},v_{2},v_{3},v_{4})$ and $\eta =(\eta _{1},\eta _{2},\eta _{3},\eta _{4})$ represent the vectors adjoint respectively.

Then

$$\begin{aligned} {\mathcal {H}}_{2}=&a_{1}S+a_{2}I+\frac{1}{2}(a_{3}c^{2}_{1}+a_{4}c^{2}_{2})+v_{1}(\Delta - \beta \dfrac{SI}{N} + \mu R - (kc_{1}+\delta )S)\\&+v_{2}(kc_{1}S- (1-\tau )\dfrac{\beta IV}{N} -\delta V)\\&+v_{3}(\beta \dfrac{SI}{N}+(1-\tau )\dfrac{\beta IV}{N}-(\alpha c_{2} +\delta +\delta _{0} )I)\\&+v_{4}(\alpha c_{2} I-(\delta +\mu )R)+\eta _{1}\rho _{1}S+\eta _{2}\rho _{2}V +\eta _{3}\rho _{3}I+\eta _{4}\rho _{4}R. \end{aligned}$$

The stochastic maximum principle confirmed that the following relation hold:

$$\begin{aligned} & dz^{\star }(t)=\dfrac{\partial {\mathcal {H}}_{2}(z^{\star }(t),c^{\star }(t),v,\eta )}{\partial v}dt+\sigma (z^{\star })dW(t),\\ & dv(t)=-\dfrac{\partial {\mathcal {H}}_{2}(z^{\star }(t),c^{\star }(t),v,\eta )}{\partial z}dt+\eta dW(t).\end{aligned}$$

$$\begin{aligned} {\mathcal {H}}_{2}(z^{\star }(t),c^{\star }(t),v,\eta )=\min _{c \in \Omega ^{*}}{\mathcal {H}}_{2}(z^{\star }(t),c(t),v,\eta ), \end{aligned}$$

(17)

with initial and terminal condition

$$\begin{aligned} z^{\star }(0)=(S_{0},V_{0},I_{0},R_{0}), \ \ v(T_{e})=-\dfrac{\partial \theta (z(T_{e}))}{\partial z},\end{aligned}$$

where

$$\begin{aligned}\theta =\frac{p_{1}}{2}S^{2}(T_{e})+\frac{p_{2}}{2}V^{2}(T_{e})+\frac{p_{3}}{2}I^{2}(T_{e})+\frac{p_{4}}{2}R^{2}(T_{e}),\end{aligned}$$

and $z^{*}(t)$ represent the optimal trajectory of z(t).

Theorem 6

Our control problem has an optimal control as the following form

$$\begin{aligned}(c^{\star }_{1},c^{\star }_{2})=\Big (\min \Bigl \{\max \Bigl \{\dfrac{kS(v_{1}-v_{2})}{a_{3}},0\Bigr \},c^{max}_{1}\Bigr \},\min \Bigl \{\max \Bigl \{\dfrac{\alpha I(v_{3}-v_{4})}{a_{4}},0\Bigr \},c^{max}_{2}\Bigr \}\Big ).\end{aligned}$$

Proof

Applying the well-know minimum principle in stochastic case, we get the adjoint system

$$\begin{aligned} dv_{1}(t)&=(-a_{1}+(v_{1}-v_{3})\beta \dfrac{I(t)}{N}+(v_{1}-v_{2})kc_{1}(t)+v_{1}\delta +\eta _{1}\rho _{1})dt+\eta _{1}dW_{1},\\ dv_{2}(t)&=((v_{2}-v_{3})(1-\tau )\beta \dfrac{I(t)}{N}+v_{2}\delta +\eta _{2}\rho _{2})dt+\eta _{2}dW_{2},\\ dv_{3}(t)&=(-a_{2}+(v_{1}-v_{3})\beta \dfrac{s(t)}{N}+(v_{3}-v_{4})(\alpha c_{2})+v_{3}(\delta +\delta _{0}))dt+\eta _{3}dW_{3},\\ dv_{4}(t)&=(-v_{1}\mu +v_{4}(\delta +\mu )+\eta _{4}\rho _{4})dt+\eta _{4}dW_{4}, \end{aligned}$$

with terminal condition

$$\begin{aligned} v_{1}(T_{e})=-p_{1}S(T_{e}), v_{2}(T_{e})=-p_{2}V(T_{e}),v_{3}(T_{e})=-p_{3}I(T_{e}),v_{4}(T_{e})=-p_{4}R(T_{e}),\end{aligned}$$

and we have

$$\begin{aligned}&\dfrac{d{\mathcal {H}}_{2}}{dc_{1}}=kS(v_{2}-v_{1})+a_{3}c_{1} =0,\\&\dfrac{d{\mathcal {H}}_{2}}{dc_{2}}=\alpha I(v_{4}-v_{3})+a_{4}c_{2} =0. \end{aligned}$$

Then

$$\begin{aligned} c^{\star }_{1}(t)&=\dfrac{kS^{\star }(t)(v_{1}(t)-v_{2}(t))}{a_{3}},\\ c^{\star }_{2}(t)&=\dfrac{\alpha I^{\star }(t)(v_{3}(t)-v_{4}(t))}{a_{4}}. \end{aligned}$$

From the lower and upper bounds of $c_{1},c_{2}$, we derived

$$\begin{aligned} c^{\star }_{1}(t)&=\min \Bigl \{\max \Bigl \{\dfrac{kS^{\star }(t)(v_{1}(t)-v_{2}(t))}{a_{3}},0\Bigr \},1\Bigr \},\\ c^{\star }_{2}(t)&=\min \Bigl \{\max \Bigl \{\dfrac{\alpha I^{\star }(t)(v_{3}(t)-v_{4}(t))}{a_{4}},0\Bigr \},1\Bigr \}. \end{aligned}$$

$\square$

Numerical simulations

Parameter estimation

In this section, we discuss the process used to estimate the parameters of our model, employing the nonlinear least-squares method. This approach is a widely recognized statistical technique that is commonly used for estimating parameters in mathematical models. It operates by minimizing the sum of the squared differences between observed data and the corresponding model predictions, ensuring the best possible fit of the model to the given data. Nonlinear least squares is particularly advantageous in regression analysis and curve fitting. To obtain accurate parameter values, we utilized a combination of previously established data and real-world COVID-19 case reports. Some of the model parameters were directly obtained from the reference⁵. The remaining parameters were determined by fitting our model to actual COVID-19 case data reported in Algeria during the period from October 22, to March 14, 2022. These data were sourced from the website²⁷.

The parameters values are listed in Table 2 and the fitted plot is depicted in figure 1.

Table 2 Values of the model embedded parameters.

Fig. 1

Observed corona virus cases in Algeria and the deterministic model predicted cumulative infected.

Sensitivity analysis

Sensitivity analysis plays a crucial role in epidemic modeling, helping to assess the impact of specific parameters on variations in the model’s variables. It quantifies the extent to which $R_{0}$ responds to changes in the model’s inherent parameters. We use the normalized sensitivity procedure to evaluate the relative indices of crucial parameters. This approach is expressed by the following formula²⁸:

$$\begin{aligned}{\mathcal {S}}_{p}=\dfrac{\partial R_{0}}{\partial p}.\dfrac{p}{R_{0}}.\end{aligned}$$

Mathematical expressions for the normalized sensitivity indexes can be expressed as follows:

$$\begin{aligned} & {\mathcal {S}}_{\beta }=1,\ \ {\mathcal {S}}_{\alpha }=\dfrac{-\alpha }{\alpha +\delta +\delta _{0}+\frac{\eta _{3}^{2}}{2}}, \ \ {\mathcal {S}}_{\delta }=\dfrac{-\delta }{\alpha +\delta +\delta _{0}+\frac{\eta _{3}^{2}}{2}}, \ \ {\mathcal {S}}_{\delta _{0}}= \dfrac{-\delta _{0}}{\alpha +\delta +\delta _{0}+\frac{\eta _{3}^{2}}{2}},\\ & {\mathcal {S}}_{\tau }=\dfrac{-\tau }{2-\tau },\ \ {\mathcal {S}}_{\rho _{3}}= \dfrac{-\rho _{3}^{2}}{\alpha +\delta +\delta _{0}+\frac{\eta _{3}^{2}}{2}}.\end{aligned}$$

As shown in Table 3 and Figure 2, model parameters with positive index suggest a direct relation with $R_{0}$. In other words raising the values of these parameters while maintaining the others unchanged will lead to an increase in $R_{0}$ indicating that COVID-19 would spread more widely within the population. Conversely, parameters with negative index imply a reduction in the infection rate. Based on the sensitivity indixes presented in Table 3, the most critical factors influencing the spread of COVID-19 include the transmission rate $\beta$, which accelerates the virus’s spread. Also, one of the most important parameters contributing to the reduction of prevalence is $\alpha$. To more effectively combat the virus, it is essential to implement an optimal strategy that focuses on improving the recovery rate $\alpha$ while simultaneously reducing the transmission rate $\beta$.

Table 3 The sensitivity index of $R_{0}$ to the parameters.

Fig. 2

Sensitivity analysis results of $R_{0}$ for model (1).

Numerical simulation examples

To show the theoretical findings, we perform simulations of the problem in this part. We applied stochastic Range-Kutta method to discretize system (1) as follows:

$$\begin{aligned} \begin{aligned} S_{j+1}&=S_{j}+\Big [\Delta - \beta \dfrac{S_{j}I_{j}}{N_{j}} + \mu R_{j} - (k+\delta )S_{j}\Big ]dt + \rho _{1}S_{j}\sqrt{ dt} B_{1_{j}}+\frac{\rho _{1}^{2}}{2}S_{j}(B_{1_{j}}-1) dt, \\ V_{j+1}&= V_{j}+\Big [kS_{j}- (1-\tau )\beta \dfrac{V_{j}I_{j}}{N_{j}} -\delta V_{j}\Big ]dt+ \rho _{2}V_{j}\sqrt{ dt} B_{2_{j}}+\frac{\rho _{2}^{2}}{2}V_{j}(B_{2_{j}}-1) dt,\\ I_{j+1}&=I_{j}+ \Big [\beta \dfrac{S_{j}I_{j}}{N_{j}}+(1-\tau )\beta \dfrac{V_{j}I_{j}}{N_{j}}-(\alpha +\delta +\delta _{0} )I_{j}\Big ]dt+ \rho _{3}I_{j}\sqrt{ dt} B_{3_{j}}+\frac{\rho _{3}^{2}}{2}I_{j}(B_{3_{j}}-1) dt, \\ R_{j+1}&= R_{j}+\Big [\alpha I_{j}-(\delta +\mu )R_{j}\Big ]dt+ \rho _{4}R_{j}\sqrt{ dt} B_{4_{j}}+\frac{\rho _{4}^{2}}{2}R_{j}(B_{4_{j}}-1) dt, \end{aligned} \end{aligned}$$

(18)

where $B_{j} \ (j=1,..,4)$ are normally-distributed ${\mathcal {N}}(0,1)$ independent random variables and dt is step size.

The following aspects will be highlighted through a number of empirical examples:

• The dynamical behavior of model (1) if $R_{0}<1$.

• The existence of stationary distributin of (1) if $\bar{R}_{0}>1$.

• The impact of control strategies on COVID-19 dynamics.

Numerical examples for extinction and stationary distribution

We assumed the values of the parameters shown in the two examples to verify the extinction and stationary distribution. The values in example 1 are used to demonstrate Theorem 2 numerically. According to this theorem, the disease will die out with probability one if $R_{0}<1$ . Figure 3 confirms this, as the parameters in this case satisfy the theorem’s condition.

We use another set of parameter values in example 2 to show the results implied by theorem 3. This set of parameters ensures that ${\bar{R}}_{0}>1$ as can be seen in Figure 4.

Example 1

We assume that: $\Delta =1530$, $\beta =0.451$, $k=0.664$ , $\mu =0.01$, $\tau =0.576$, $\delta _{0}=0.795$, $\delta =0.65$, $\alpha =0.791$, $\rho _{1}=0.2,\rho _{2}=0.1, \rho _{3}=0.3, \rho _{4}=0.2$, and initial values $S(0)=500 , V(0)=50 , I(0)=30 , R(0)=0$, we get $R_{0}<1$, then the epidemic will be extinct from the population.

Example 2

Taking another parameter’s values : $\Delta =1530$, $\beta =0.451$, $k=0.664$, $\tau =0.057$, $\delta _{0}=0.127$, $\mu =0.01$, $\delta =0.65$, $\alpha =0.179$ and we reduce the noise intensity values: $\rho _{1}=0.1;\rho _{2}=0.1;\rho _{3}=0.1;\rho _{4}=0.1$, and the same initial values of example 1. Based on these parameters we find ${\bar{R}}_{0}>1$, by theorem 3 model (1) has a unique stationary distribution, biologically this means that the disease will persist in the population and will not be naturally eliminated under the present condition.

Fig. 3

A comparative deterministic and stochastic trajectories of solution of model 2 and model 1 respectively with parameter values of Example 1.

Fig. 4

A comparative deterministic and stochastic trajectories of solution of model 2 and model 1 respectively with parameter values of Example 2.

Through the integration of theoretical analysis and numerical simulations, we discovered that smaller white noise enables model (1) to exhibit a unique ergodic stationary distribution when ${\bar{R}}_{0}>1$. Conversely, larger white noise can lead to disease extinction in model (1) when $R_{0}<1$. In comparison to the deterministic model, the inclusion of random white noise in the epidemic model significantly influences the persistence and extinction of the desease, thereby enriching the dynamic behavior of the epidemic model.

Numerical example on control strategies of COVID-19

In this part, we will present the simulation results of the optimal control model. The Rung-Kutta method will be used to provide the numerical solution for the control problem. In paticular, we will discretize the models, the adjoint systems and the control models while adding the control caracterisation and saubsidiary conditions.

The parameters in Table 1 were used for the simulation in the presence and absence of control. The results are as shown in Figs. 5 and 6.

Fig. 5

Simulations of the deterministic and stochastic models describing dynamics of susceptible and vaccination classes in both with and without optimal controls.

Fig. 6

A comparative deterministic and stochastic simulation of infections and recovered individuals under constant and variable controls.

Figures 5 and 6 show simulation of deterministic and stochastic models describing the dynamics of S, V and I, R classes respectively, both with and without optimal control strategies ($c^{*}_{1}$ and $c^{*}_{2}$). We note that the implementation of vaccine control and treatment control leads to a decline in the S and I populations while increasing the V and R groups. We can interpret this as the vaccination reduces the number of susceptible by moving them into the vaccinated class. Simultaneously, the optimal treatment helps I individuals recover more quickly, decreasing the I population and increasing the number of R.

Through these results we concluded that time dependent treatment control and vaccine control play crucial roles in COVID-19 epidemic models, shaping efficient public health strategies to mitigate the virus’s impact. Treatment control, which includes antiviral drugs and supportive care, reduces disease severity and shortens recovery time, thereby alleviating pressure on healthcare systems. On the other hand, vaccine control decreases the number of S, lowers transmission rates, and prevents severe cases, contributing to herd immunity and long-term epidemic control. Public health policies must integrate both approaches, ensuring widespread vaccination while maintaining effective treatments to manage emerging variants and protect high-risk populations. A balanced strategy combining vaccination and treatment is essential for sustainable disease control and preparedness for future pandemic.

Conclusion

This paper aimed to presents a comprehensive mathematical framework for optimal control of the COVID-19 in stochastic case using a robust computational mathematical modeling approach. By integrating stochastic elements, implementation with empirical data, and sophisticated optimization techniques, our approach offers a potent tool for modeling, analyzing, and effectively minimizing the spread of COVID-19 outbreaks. The model in stochastic case is rigorously analyzed in theoretical point of view. The theoretical results are validated graphically using the estimated parameters. Moreover, a stochastic optimal control model was developed by incorporating time-dependent control variables for the necessary mitigation of future outbreaks. The necessary optimality conditions and solution of the stochastic control system was presented in detail. The insights gleaned from this study hold immense value for policymakers and public health authorities, empowering them to devise evidence-driven strategies aimed at curbing the transmission of the virus, safeguarding public health, and fostering societal well-being.