Stationary distribution of stochastic COVID-19 epidemic model with control strategies

1.   Introduction

Coronavirus SARS-CoV-2 (CoV-2) or the acute respiratory syndrome (SARS) is a respiratory virus. It's known as COVID-19 because of the virus that causes this disease (coronavirus disease 2019). COVID-19 was found in humans for the first time. During the middle of December, 2019, the epidemic was identified in the Chinese city of Wuhan in the province of Hubei. On January 30, 2020, the World Health Organization (WHO) stated that the SARS-CoV-2 outbreak had reached pandemic proportions and had been declared a global health emergency. Quarantine was put in place on January 23, 2020, by the Chinese government in Wuhan, China, to restrict the spread of the pandemic ^[1]. Mathematical modeling is a powerful technique for describing and dealing with diverse phenomena in nature. Recently, there has been a lot of focus on the development of mathematical models for understanding infectious diseases ^[2,3]. Many authors have established different epidemic models for the development and control of transmissible disease in communities. When compared to cardiovascular disease, infection-related diseases rank as the second leading cause of death on the globe. Mathematical models have become an important part of infectious disease epidemiology. In the twenty-first century, infectious diseases will become more prominent in both developed and developing countries. Global health challenges associated to infectious diseases are at an all-time high on the priority lists of world leaders, public health officials, and philanthropists, and they are expected to continue to rise in importance. Many researchers have spent the last few years using various ways to investigate infectious diseases and their mechanisms ^[4,5,6,7]. Aside from controlling the transmission of infectious diseases, this helps to prevent them in daily life. Epidemiological models are used to consider the advancement of contagious illnesses in populations. Many researchers have studied epidemic models to investigate and analyze diseases including avian influenza, hepatitis B, tuberculosis, and leishmaniasis ^[8,9,10]. We cannot use a single model to describe the entire disease system around the world since the existence or eradication of COVID-19 is dependent on so many different characteristics of the affected system. Since COVID-19, the human population has been quarantined, as the spread of the disease has been directly linked to this action. In most cases, quarantine is divided into two categories: One is infected whereas the other is susceptible quarantine. In this research work, we use the term "infected quarantine, " which refers to the practice of isolating those who have been exposed to an infectious disease and placing them in isolation.

Mathematical models have become an important part of infectious disease epidemiology. It has become popular to use mathematical modeling to study infectious illnesses that are communicable (see, for example, ^{[11,12,13,14]}). Mathematical modeling has been observed in the exploration of preventive mechanisms and propagations dynamics in recent years. Epidemiological models are used to take the spread of infectious diseases within the population into account ^[15,16]. One major area of concern for mathematical models in epidemiology is qualitative analysis. Deterministic and stochastic models distinguish between two categories of epidemic models. In comparison to its deterministic counterparts, the stochastic system is the most effective epidemic system for these forms and may even have a higher level of realism ^[17,18]. Because a stochastic model can be run multiple times to build up a distribution of the predicted results, as opposed to a deterministic model which can only yield a single predicted value, stochastic models yield more valuable results ^[19,20,21]. The stochastic nature of infectious disease transmission arises from the unpredictable nature of human contact.

Biological phenomena are often influenced by environmental noise in the real world. Beyond epidemiology, stochastic theory and techniques are widely used in other fields of nonlinear applied sciences. The stochastic approach has gained popularity due to its ability to accurately represent randomness and uncertainties in real-world scenarios. As a result, stochastic epidemic systems have been extensively studied (see ^[22]). Stochastic variations have long intrigued mathematicians and researchers because they can lead to significant changes that better align with real-world problems. For example, Rihan and Alsakaji ^[23] introduced a stochastic susceptible, infected, asymptomatic, quarantined, recovered $(SIAQR)$ epidemic model for COVID-19 with time delays. In addition to qualitative analysis, they used reported data from the United Arab Emirates to predict the disease's behavior. Similarly, the authors in ^[24] proposed a stochastic susceptible-infected-recover-cross immune $(SIRC)$ model for COVID-19, establishing sufficient conditions for the existence of an ergodic stationary distribution. They also showed that while diseases may go extinct under high levels of white noise, recurrence and periodic outbreaks are still possible due to the time-delayed feedback in transmission dynamics. The epidemic models can display qualitatively different dynamical behaviors when the conventional presumption of mass-action law is discarded. Numerous infectious illnesses like COVID-19 exhibit periodic changes in their prevalence. Such periodicities may be driven by extrinsic factors, as reflected in periodic transmission rates, e.g., seasonality ^[25], or may be caused by time delays (e.g., ^[26,27,28]) or nonlinearity of incidence rates. To describe the dynamics of COVID-19, the authors in ^[28] used the stochastic epidemic model with bilinear incidence rate and time-delay. In this study, we investigate the effects of weakening the bi-linearity assumption. As the disease COVID-19 is spreading at very fast rate with an increasing infected population, it is more realistic to consider the incidence of the form $\frac{\beta SI(\delta I+1)}{N(t)}$. Here, the notion $N(t)$ is the total population which is further divided into five disjoint compartments, namely, vulnerable $S(t)$, exposed $E(t)$, infected $I(t)$, quarantined $Q(t)$, and recovered individuals $R(t)$. To the best of our knowledge, in case of COVID-19, researchers have given less attention to the stochastic modeling with such type of incidence rates. Further, the authors discussed the qualitative analysis and controlling strategies that were completely ignored. Therefore, in this research work, we aim to prove the considerable effect of a stochastic factor on the dynamics and control of the COVID-19 epidemic model. Our model assumes that the total population is constant and that the parameters governing disease transmission and recovery are homogeneous across the entire population. This simplification may not accurately reflect real-world variations, such as changes in population behavior or differing rates of disease transmission across regions. The study relies on specific datasets for parameter estimation and model validation. These datasets are limited to certain regions and time frames, which may not capture all variations in disease dynamics. Consequently, the results may not be fully generalizable to other contexts or populations. While our model effectively captures the impact of stochastic noise on disease dynamics, it is sensitive to the choice of noise intensity parameters. Small changes in these parameters can lead to significantly different outcomes, which may pose challenges in real-world applications where precise parameter estimation is difficult. The stochastic nature of the model requires extensive computational resources for simulations, especially when exploring a wide range of noise intensities and control strategies. This may limit the model's applicability in scenarios where real-time analysis is required. The model and its conclusions are based on certain disease characteristics and may not be directly applicable to other infectious diseases with different transmission mechanisms or population structures. Further validation with different diseases and contexts is necessary to extend the findings.

Chen and Kang ^[29] presented a stochastic vaccination model with backward bifurcation and found that under some conditions, the smaller the intensity of the random perturbations is, the smaller the distance between the solution of the stochastic system and the stable equilibrium of the corresponding deterministic one. Cai et al. ^[30] discussed a stochastic susceptible-infected-recover $(SIRS)$ epidemic model with nonlinear incidence rate and found stochastic perturbations can suppress the disease outbreak. Zhang et al. ^[31] presented a stochastic SIRS model with standard incidence rate and partial temporary immunity. Chang et al. ^[32] presented a stochastic SIRS epidemic model with two different saturated incidence rates and got the thresholds which leads to the extinction or persistence of the disease. These researches reveal that random fluctuations in the environment can restrain the spread of the disease. In other words, the transmission capacity of infectious diseases described by the deterministic model that ignores random perturbations is greater than their true transmission capacity. In this paper, our concern is when random fluctuations are introduced to the proposed system and whether they can also restrain the spread of the disease.

The paper is organized as follows. We developed the stochastic COVID-19 epidemic model in Section 2. The existence of global nonnegative solutions as well as an analysis of their uniqueness are examined in Section 3. Section 4 provides the necessary criteria for a stationary distribution to exist. In Section 5, we deduce conditions of extinction. Moreover, we discuss optimal control strategies in Section 6. To bolster our work in Section 7, numerical work in the form of simulations has been done. Finally, we wrapped up our work in Section 8.

2.   Mathematical model

To formulate the model, we will extend the model of Khan et al. ^[28] by incorporating the incidence rate of the form $\frac{\beta SI(\delta I+1)}{N(t)}$ and we will ignore the delay effect due to the fast spreading nature of the disease. We introduce noise as a stochastic term in our model, where it plays a critical role in capturing the random fluctuations in the system. In our normalized population model, the total population is scaled to unity, meaning that each subpopulation (e.g., susceptible, infected) represents a fraction of the total population. We define excessive noise as a level of stochastic fluctuation that becomes significant compared to the size of the subpopulations. Since the total population is normalized, noise that exceeds the corresponding subpopulation levels can cause unmanageable fluctuations, leading to dynamics that may no longer be reflective of the real-world scenarios the model aims to capture. For example, if the susceptible population constitutes 30% of the total population, noise levels approaching or exceeding this proportion would be deemed excessive. At this point, the noise introduces random fluctuations that overwhelm the deterministic trends in the model, potentially leading to unrealistic results. To capture the dynamics of COVID-19, we will use five ordinary differential equations and initially we formulate a deterministic mathematical model. The overall population is divided into five compartments i.e., $S(t)+E(t)+I(t)+Q(t)+R(t) = N(t)$ with $t$ as the independent variable. Following are the equations that illustrate the model:

The flow for this model is depicted in Figure 1. We assumed the incidence rate of the form $\frac{\beta SI(\delta I+1)}{N}$ because the infection COVID-19 spread too quickly as the infected population tends to increase. Here, the notion $\delta$ is a positive constant and biologically, it is a balance factor for the infected population. This incidence function could reduce the mass action law $\beta SI$ if we set $\delta = 0$, which means the bilinear incidence rate is a special case of this function. To make the problem dimensionally homogenous, we must assume the unit per individual for the parameter $\delta$. The description for the rest of the parameters are presented in Table 1.

The stochastic form of the deterministic system (2.1) is represented as follows:

Here, $W_{i}(t), i = 1, ..., 5$ are independent standard Brownian motions, and $\alpha_{i}, i = 1, ..., 5$ are the intensities of the standard Gaussain white noises, respectively. The interaction between the environment and individuals are presented by $\alpha_{1} S dW_1(t), \alpha_{2} E dW_2(t),$ and $\alpha_{3} I dW_3(t),$ $\alpha_{4} Q dW_4(t),$ and $\alpha_{5} RdW_5(t)$.

3.   Existence and uniqueness

In Section 3, we now emphasize that Theorem 2 establishes the existence and uniqueness of the solution to the stochastic differential equations that describe the model. Specifically, the theorem ensures that under the given assumptions, the stochastic system admits a unique stationary distribution, which implies long-term persistence or eradication of the disease, depending on the model parameters.

Additionally, we have clarified that the existence of the model refers to the mathematical formulation of the system, while existence and uniqueness of the solution refer to the behavior of the system over time under stochastic influences.

Theorem 1. There is a unique solution $(S(t), E(t), I(t), Q(t), R(t))$ of system (2.2) on $t\geq 0$ for any initial value $\zeta(0) = (S(0), E(0), I(0), Q(0), R(0))\in \mathbb{R}^{5}_{+}$., and the solution will remain in $\mathbb{R}^{5}_{+}$ with probability one, namely, $(S(t), E(t), I(t), Q(t), R(t))\in \mathbb{R}^{5}_{+}$ for all $t\geq 0$ almost surely.

Proof. For any $\zeta(0)\in R^{5}_{+}$, the coefficients of our model accurately fulfill the Lipschitz local criterion. Consequently, $\zeta(t)$ represents a unique local solution for $t\in[0, \tau_e)$, where $\tau_e$ denotes the explosion time (^[27]). We shall now demonstrate that the solution is still globally valid and establish $\tau_e = \infty$ a.s. When $k_0$ is large enough to ensure that $k_0\geq 0$, the starting approximation $S(0)$, $E(0)$, $I(0)$, $Q(0)$ and $R(0)$ are in $[\frac{1}{k_0}, k_0]$. We calculate the stopping time for $k\geq k_0$, as

In this study, we employ the concept of $inf\phi = \infty$, where $\phi$ represents null set, in accordance with the work of ^[33]. Since $\tau_k$ increases whenever $k$ approaches to $\infty$, thus, $\lim_{k\rightarrow \infty}\tau_k = \tau_\infty$ along the application of $\tau_\infty\leq\tau_e$ a.s showing that $\tau_\infty = \infty$ a.s. Accordingly, this will ensure the that solution of model (2.2) lies in $R^{5}_{+}$ a.s., $\forall$ $0\leq t$. To show that $\tau_e = \infty$ a.s, if this assertion is false, then there exist a pair of constants $T > 0$ and $\epsilon\in (0, 1)$ such that

As a result, there is an integer $k_1 \geq k_0$ such that

Consider a $C^{2}$-mapping $H:\mathbb{R}^5_{+}\rightarrow \bar{\mathbb{R}}_{+},$ where $\bar{\mathbb{R}}_{+} = \{x\in \mathbb{R}:x\geq 0\},$ by

Using the $It\hat{o}'s$ formula in Eq (3.3) gives us

The differential operator $L$ associated to $H$ in the above relation is given by

where

Thus in Eq (3.4), $LH:\mathbb{R}^5_{+}\rightarrow \mathbb{R}_{+}$, and we can write

Since $K$ is a positive constant therein independent of $S, E, I, Q, R,$ and $t$, we can get

Integrating both sides (3.6) from 0 to $T \wedge \tau_k$ and taking expectations, we can obtain

$\forall$ $k \geq k_1$, assume $T = \Omega_k$ & ${T\geq\tau_k}$, then $P(\Omega_k)\geq\epsilon$. Utilizing $\omega\in \Omega_k,$ for each, we find $S(\omega, \tau_k)$, $E(\omega, \tau_k)$, $I(\omega, \tau_k)$, $Q(\omega, \tau_k)$, $R(\omega, \tau_k)$ values of which are based on $k$ or $\frac{1}{k}$. Therefore, $H(S(\tau_k)$, $E(\tau_k, \omega))$, $I(\tau_k)$, $Q(\tau_k)$, $R(\tau_k)$ is not less than

Consequently,

From Eqs (3.7) and (3.8), we have

Here, $1_{\Omega(\omega)}$ is used to represent the function indicator for $\Omega_k$. The contradiction develops when k approaches infinity, that is, $\infty > H\big(S(0), E(0), I(0), Q(0), R(0))\big)+KT = \infty,$ showing that $\tau_\infty = \infty$ a.s.

4.   The stationary distribution and extinction

We investigate the system's proposed problem of stationary distribution $(1)$. It really is easy to make a deterministic system by putting $\alpha_{1} = \alpha_{2} = \alpha_{3} = \alpha_{4} = \alpha_{5} = 0$, the present version; on the other hand, is separate from its deterministic counterpart. Furthermore, we are aware that there is no endemic equilibrium in the stochastic system. Because of this, a linear stability analysis cannot be used to investigate the disease's permanence; instead, we look at the system's stationary distribution. Finally, this indicates that the illness will continue to exist. We shall use a well-known Khasminskii outcome for this ^[34]. Suppose

and for $M = \int_0^t u(s)dB(s)$, we have

Lemma 1. By (Strong Law for large numbers)^[35], if $F = \{F\}_{t\geq 0}$ is a real-valued continuous local martingale which goes away at $t = 0$, then

4.1.   Stationary distribution

Different from the deterministic system (1.1), the stochastic system does not have the endemic equilibrium. Hence, we cannot study the persistence of the disease by studying their stability of the endemic equilibrium, and we turn to research the existence and uniqueness of the stationary distribution for the system (1.2) which implies the persistence of the disease in some sense. To this end, we cite a well-known result from Hasminskii ^[34]. Let $X(t)$ be a regular time-homogeneous Markov process in Rn + described by

This is the display of the diffusion matrix:

Lemma 2. By ^[34], the unique stationary distribution $m(\cdot)$ of the Markov process $X(t)$ is dependent upon the existence of a bounded domain $U \in \mathbb{R}^d$ with a regular boundary and its closure $\overline{U} \in \mathbb{R}^d$, satisfying the following properties:

$1)$ Within the open domain U and its surrounding areas, the diffusion matrix A(t)'s smallest eigenvalue is bounded away from zero.

$2)$ If $x \in R^d U,$ then $Sup_{x\in k} E^x\tau < \infty$ for each compact subset $K\subset R^n$, and the time average $"\tau"$ at which a region originating from x goes to the finite set U. Furthermore, if $f (.)$ is a function that can be integrated with respect to $\pi$, then

Disease-free equilibrium in the deterministic model & the concerned stochastic system has a threshold value is given by

Following ^[36] and keeping in view the expression from $R^d_0$, we can define a threshold parameter for the stochastic system of the form

Theorem 2. If $R^S_0 > 1$, then the root $(S(t), E(t), I(t), Q(t), R(t))$ of model (2.2) is ergodic and there is only one stationary distribution $\pi(.)$.

Proof. We demonstrate that, under condition (2) of Lemma 2, we need to construct a nonnegative $C_2-$function $V:R^5_+ \rightarrow R_+$. To this end, we first define

where it is necessary to identify the values of the three constants referred to in the problem, namely, $c_{1}$, $c_{2}$, and $c_{3}$. We can calculate system (2.2) using Itô's formula.

Therefore, we have

The above implies that

Let

namely,

Consequently,

Furthermore, we acquire

where $c_{4} > 0$ represents a constant to be determined later. It can be readily obtained that

where $(\frac{1}{p}, p) \times (\frac{1}{p}, p) \times(\frac{1}{p}, p) = U_p.$

Following that, we will show that $V_2(S, E, I, Q, R)$ has unique minimum value $V_2(S_0, E_0, I_0, Q_0, R_0).$ With respect to $Q, R, I, S, E$, the partial derivative of $V_2(S, E, I, Q, R)$ is as follows:

It is easy to demonstrate that $V_2$ has a only one stagnation point.

Furthermore, the Hessian matrix of $V_2(S, E, I, Q, R)$ at $(S(0), E(0), I(0), Q(0), R(0))$ is

This shows that $\Lambda$ is obviously a nonnegative definite matrix. Therefore, $V_2(S, Q, E, I, R)$ has a minimum value $V_2(S_0, E_0, I_0, Q_0, R_0).$

According to Eq (4.9), and the continuity of $V_2(S, E, I, Q, R)$, we can conclude that in the interior of $V_2(S, E, I, Q, R)$ has a only one minimum value $V_2(S(0), E(0), I(0), Q(0), R(0))$ belonging to $\mathbb{R}^5_+$.

Moreover, we take into account a nonnegative $C^2-$ mapping $V : \mathbb{R}^5_+ \rightarrow \mathbb{R}_{+}$ as follows:

Ito's formula can be used to derive the proposed model.

the conclusion of which is as follows

where

The set $D_i$ represents a sequence of domains used to control the behavior of the solution and ensure it remains within a feasible region. The set $D_i$ is chosen based on the dynamics of the system to prevent the solution from approaching critical boundaries, such as zero or infinity. The following step is to define the set.

where $\varepsilon_{i} > 0(i = 1, 2, 3, ......10)$ indicate small constants that will be evaluated hereafter.

To investigate the convenience, $\mathbb{R}_{+}^{5} \backslash D$ can be divided into the ten domains listed below,

Next, we will demonstrate that $LV(S, E, I, Q, R) < 0$ on $\mathbb{R}_{+}^{5} \backslash D$ using the same method from the ten domains mentioned above.

Case 1. According to Eq (4.11), if $(S, E, I, Q, R) \in D_1$, we acquire

Choosing $\epsilon_1 > 0$ yields

Case 2. By Eq (4.11), if $(S, E, I, Q, R)\in D_2$, we acquire

Let $\epsilon_1^2 = \epsilon_2$. We can chose sufficiently enough large $c_4 > 0$ and a small enough $\epsilon_1 > 0$ so that we can obtain

Case 3. By Eq (4.11), if $(S, E, I, Q, R)\in D_3$, we acquire

We can set a small enough $\epsilon_1 > 0$ so that we can obtain

Case 4. By Eq (4.11), if $(S, E, I, Q, R)\in D_4,$ we acquire

We can set a small enough $\epsilon_2 > 0$ so that we can obtain

Case 5. According to Eq (4.11), if $(S, E, I, Q, R) \in D_5,$ we acquire

We choose a small $\epsilon_2 > 0$ so that we can obtain

Case 6. By Eq (4.11), if $(S, E, I, Q, R)\in D_6,$ we acquire

We may select a small enough $\epsilon_2 > 0$ so that we can obtain

Case 7. By Eq (4.11), if $(S, E, I, Q, R)\in D_7,$ we acquire

We choose a small $\epsilon_2 > 0$ so that we can obtain

Case 8. If $(S, E, I, Q, R)\in D_8,$ from Eq (4.11), we have

We can choose a small $\epsilon_2 > 0$ so that we can obtain

Case 9. By Eq (4.11), if $(S, E, I, Q, R)\in D_9,$ we acquire

If $\epsilon_1 > 0$, then we can find

Case 10. If $(S, E, I, Q, R)\in D_10,$ from Eq (4.11), we have

We can choose sufficiently small $\epsilon_2 > 0$ so we can obtain

In conclusion, from 4.13 to 4.21, we can see that there exists a positive constant $W > 0$ such that

Hence,

Let $(S(0), E(0), I(0), Q(0), R(0)) = (x_1, x_2, x_3, x_4, x_5) = x \in R^5_+ \backslash D,$ the time $\tau^x$ when a region originating from x approaches to D, $min\{\tau^x, t, \tau_k\} = \tau^{(k)}(t),$ and $\tau_k = \inf\{t\; \; such that\; \; |X| = k\}$.

Using Dynkins formula, we take expectation and integrate both sides of the Eq (4.22) from zero to $\tau^{(k)}(t)$,

Since $V(x)$ is positive, hence

The proof of the outcome from Theorem 2 implies that $P(\tau_e = \infty) = 1$. In other words, the system (2.2) is regular.

Therefore, we have that $\tau{(k)}(t)$ tends to $\tau^x$almost surely as $k$ tends to $\infty$ and $t \rightarrow \infty.$ As a result, with the assistance of Fatou's lemma, we arrive at

Obviously, $\sup_{x\in K} \mathbb{E}\tau^x < \infty$, where $K$ denotes a compact subset of $\mathbb{R}_{+}^{5}$. Consequently, condition (2) of Lemma 2 holds. Moreover, the model (2.2) diffusion matrix is provided by

Consider $M = \min_{(S, E, I, Q, R)\in \overline{D}\in R^5_+}\{\sigma^2_1Q^2, \sigma^2_2 E^2, \sigma^2_3I^2, \alpha^2_4S^2, \alpha^2_5R^2\},$ and we have

where $\xi_{i}, i = 1, 2, 3, 4, 5 = \xi \in \mathbb{R}_{+}^{5}.$

By applying Lemma 2 and verifying that condition (1) of Lemma 2 holds, we have established that the diffusion matrix of the system is nondegenerate and bounded away from zero in the domain D. Furthermore, we have shown that the stopping time is $\tau_k \to \infty$ as $k \to \infty$, ensuring that the system's solution remains in the positive region $\mathbb{R}^5_+$ with probability one. Thus, conditions (1) and (2) of Lemma 2 are satisfied, implying the existence of a unique stationary distribution for the stochastic system. This stationary distribution reflects the long-term persistence of the disease within the population. Therefore, the statement of Theorem 2 is fully proved.

5.   Extinction

Concerning the disease's extinction, the following information is available to us.

Lemma 3. Assume a solution $(s(t), E(t), I(t), Q(t), R(t))$ of system (2.2) with starting condition $(S(0), E(0), I(0), Q(0), R(0)) \in \mathbb{R}_+^5$, then

Further, the root of the model has the following properties:

and

Furthermore, when $\mu > \frac{1}{2}(\alpha_1^2\vee\alpha_2^2\vee\alpha_3^2\vee\alpha_4^2\vee\alpha_5^2)$ holds, then

Proof. The proof follows the same approach as the proof of Lemma 4.1 in ^[37]. Therefore, we omit the details here.

Theorem 3. Let $\mu > (\frac{1}{2})(\alpha_1^2\vee\alpha_2^2\vee\alpha_3^2\vee\alpha_4^2\vee\alpha_5^2)$, then $(R(0), S(0), Q(0), I(0), E(0)) \in \mathbb{R}_+^5,$ and if

holds, then

Moreover,

Proof. Define a function $W_0,$ which is differentiable as follows:

According to Ito's formula and model (2), we acquire

Divide by $t$ on both sides of (5.8) and integrate from $0$ to $t$. Also by Lemma 3, we have

which illustrates that

Furthermore, the first equation of system $(2)$ can be divided by $t$ on both sides and integrated from zero to $t$ to get a solution.

which follows if we examine (5.10) and Lemma 3

The proof for Theorem 3 has been completed.

6.   Optimal control

Stochastic control of system (2) will take the same form as before, taking the same control variables into account.

with the initial data

We establish a vector of the form for the purpose of simplicity.

and

where the time functions $z_{i}$ and $u_{i}$ are used. It is also possible to represent the initial data in the form of

The following vector components are represented by the functions g and f.

$g_1 = \alpha_1 S, \; \; g_2 = \alpha_2 E, \; \; g_3 = \alpha_3 I, \; \; g_4 = \alpha_4 Q, \; \; g_5 = \alpha_5 R.$ In addition, we take into account the quadratic cost functional

where $A_1, A_2, A_3, k_i for i = 1, 2, ...5$ are nonnegative constants. In order to achieve our main objective, we must find a control $u^*(t) = (u_{2}^*(t), u_{1}^*(t))$ that allows us to perform this.

In this case, the set U denotes permissible controls.

Here, $u_j^{\max},$ where $j = 1, 2$ is real and nonnegative. ˜Hamiltoniana for the given system must first be defined before the stochastic maximum principle can be applied to the system.

where the inner product space of Euclidean is represented by $\langle., .\rangle$ while the adjoint vectors systems $[m_1, m_2, m_3, m_4, m_5]' = m$ and $[n_1, n_2, n_3, n_4, n_5]' = n$ are different. The following is a direct expansion of the maximum principle

where $z^*(t)$ represent the optimal path for $z(t).$ The conditions for the start and end of Eqs (6.4) and (6.5) are

respectively. There is a mapping between "m(t), n(t), " and "z*(t)" in the optimal control, as shown in the Eq (6.6).

As a result, the Hamiltonian in the given case is

It is suggested by the stochastic maximum principle that

We obtain

The auxiliary starting and ending criteria are also listed.

and

where $m_1(t_f) = -k_1S, \; \; m_2(t_f) = -k_2E, \; \; m_3(t_f) = -k_3I, \; \; m_4(t_f) = -k_4Q, \; \; m_5(t_f) = -k_5R, .$ Now, by differentiating the Hamiltonian equation with respect to u1 and u2, we get the following optimal controls $u_1^*$ and $u_2^*$:

7.   Numerical simulation

In the context of our stochastic modeling, we recognize the importance of providing a more comprehensive view of the inherent variability in our system. As such, we have incorporated confidence intervals in the figures to quantify the uncertainty around our estimates. Confidence intervals represent the range within which the system's state variables are expected to fluctuate due to the stochastic effects present in the model.

We will include a detailed explanation of how the level of noise is quantified relative to the sub-populations. Specifically:

● Low noise: Minor fluctuations that don't significantly affect the deterministic dynamics.

● Moderate noise: Fluctuations that introduce variability but still allow the population to follow predictable patterns over time.

● Excessive noise: Fluctuations on the same order of magnitude as the subpopulations (or higher), where random events can push the infected population toward extinction or cause erratic behavior.

For each simulated trajectory, we computed confidence intervals based on multiple realizations of the stochastic process. These intervals offer a visual representation of the range of potential outcomes, allowing for a better understanding of the variability in the disease dynamics. Specifically, 95% confidence intervals have been included in all figures to indicate the range in which the true values of the state variables lie, with 95% certainty.

Additionally, we have added stationary distributions to represent the long-term behavior of the system under stochastic influences. The stationary distribution reflects the equilibrium state that the system tends to settle into over time, despite random fluctuations. This distribution is crucial for understanding the persistent nature of the epidemic under stochastic dynamics.

The inclusion of both confidence intervals and stationary distributions enriches the analysis by highlighting the impact of stochasticity on the system's behavior, providing a clearer picture of both the short-term dynamics and the long-term tendencies of the epidemic model. In the following simulations, we consider three distinct cases to explore the different behaviors of the epidemic model under varying parameter settings. Each case corresponds to different combinations of the infection rate, recovery rate, and other model parameters:

● Case 1: Low infection rate and high recovery rate, simulating conditions leading to disease eradication.

● Case 2: Moderate infection rate and recovery rate, representing an endemic scenario where the disease persists.

● Case 3: High infection rate and low recovery rate, modeling a severe outbreak scenario.

These cases allow us to study the impact of different epidemiological factors on the dynamics of disease spread and control. Based on the deterministic $SEIQR$ system, this stochastic system is a coupled system (2.1). Stochastic first order Runge Kutta approach is used to simulate numerical solution for system (2.2). Initially, in this section, we will show how to derive the stochastic Runge Kutta scheme for the system (2.2).

Here, $\Delta W_{i, t_{n}} = W_{i, t_{n+1}}-W_{i, t_{n}}$, for $i = 1, ..., 5$ represents the independent increments of the Gaussian Brownian motion and $\Delta t_n = t_{n+1}-t_n$ represents the time increment. In our case, we restrict ourself to a constant time step $\Delta t_n = \Delta t$. 1000 equally spaced time increments are used to divide the interval. We numerically solve the $SEIQR$ system $(2)$ using the theoretical results shown above, with a number of randomly generated initial conditions. We run our code for generating three simulations. The first one is based on the choice of the parameter used in ^[28]. Similar choices of the parameter are used for the second and the third tests up to the parameter $\beta$. A list of the parameters that were used is provided in Table 2. For the mean solution, we generate 1000 realizations. The following correlation coefficients $\alpha_{i}$ for $i = 1, 2, ..., 5$ with values in the range $(0, 1)$ are set for the three simulations by applying randomly chosen values to the uniform random generator. The following three scenarios are examined (see Table 2).

All of the individuals' initial values are chosen at random in the range of zero to one in the closed region $[0, 1]$. For each step time, we establish the normatization

It should be noted that model (2) operates based on five distinct independent white noises denoted as $\Delta W_{i}(t)$ for $i = 1, 2, ..., 5$. The fourier series is being used to approximate the multiple stochastic integrals in order to assure the first order of our numerical approach ^[38].

Based on these values of parameters, we estimated both the parameters $R_0^d$ and $R_0^s$. It was obsreved that both the parameters are greater than one, therefore, the endemic equilibrium must exist in case of the deterministic model and the stochastic model should have a unique stationary distribution. To verify these facts numerically, we plotted Figures 2–4, which shows sample realizations of the stochastic $SEIQR$ model (2.2). The corresponding mean solution is generated out of 1000 realizations. In the end, we concluded that outcomes in system (2.1) were satisfied by the simulations of Tests 1–3, namely, $(S(t), E(t), I(t), Q(t), R(t)) \in \mathbb{R}_+^5$ for some $0 \leq t$. In addition, all $SEIQR$ numerical stability has been verified in a series of tests (2.2). One can notice from the figures that the stochastic curves fluctuate around the endemic equilibrium, which in turns, shows the global stability of the endemic equilibrium of the associated deterministic system. Simulated trajectories of the infected population with 95% confidence intervals (shaded area) and corresponding stationary distribution, and simulated trajectories of the recovered population with 95% confidence intervals and stationary distribution are shown.

8.   Conclusions

In this study, we have developed a novel stochastic epidemic model and demonstrated the existence of an ergodic stationary distribution, highlighting the role of noise in affecting disease dynamics. Unlike traditional deterministic models, our stochastic approach provides a more realistic representation of real-world uncertainties by capturing the randomness in disease transmission and recovery processes. This novel framework shows that noise can significantly influence the persistence of an epidemic and, in extreme cases, lead to extinction or recurrence of outbreaks. These findings contribute to the growing body of research on stochastic epidemic models, offering new perspectives on how randomness shapes epidemic outcomes.

Our study builds upon the extensive work done in epidemic modeling, particularly on the stochastic $SEIR$ and $SEIQR$ models. Previous works, such as those by Rihan and Alsakaji ^[23], Zhang et al. ^[31], and Cai et al. ^[30], introduced stochastic perturbations and nonlinear incidence rates to capture the uncertainties in disease transmission. However, these models often did not incorporate quarantine dynamics, which are crucial for COVID-19.

In contrast, our model extends this body of work by incorporating a quarantine class alongside stochastic perturbations driven by Brownian motion. The addition of quarantine in our $SEIQR$ model provides a more realistic framework for understanding the feedback loop between noise intensity, infection rates, and quarantine measures. Similar to the findings by Zhang et al. ^[31], where stochastic noise could suppress disease spread, our results indicate that high noise levels can lead to disease extinction. However, our study goes further by illustrating the pivotal role of quarantine in this dynamic, which has not been explored in-depth in previous models.

By using the Lyapunov function, we also provide mathematical guarantees for the existence of ergodic stationary distributions, adding a layer of robustness to our results. This advances the field by offering more concrete tools for understanding the long-term behavior of epidemic systems under stochastic influences. Our findings suggest that noise and quarantine strategies, when properly managed, can significantly alter the trajectory of an epidemic, offering valuable insights for public health policymakers.

Future work: Moving forward, several extensions to this study can be explored. First, future work could investigate the impact of different types of noise, such as discrete-event noise or environmental fluctuations, which may offer a more accurate representation of real-world uncertainties. Additionally, the model could be expanded to consider heterogeneous populations or spatially structured networks, where the contact patterns between individuals vary based on factors like age, location, or social behavior. This would provide a more nuanced understanding of how epidemics spread in complex social systems.

Moreover, future research could focus on the transient dynamics of the model, particularly during the early stages of an outbreak, to better inform intervention strategies. Finally, validating the model with more case studies from different epidemic scenarios and diseases would enhance its applicability and robustness. By addressing these avenues in future work, the current model could be further refined to provide more comprehensive insights into epidemic dynamics under uncertainty.

Author contributions

R. Ikram: Formal analysis, Writing–original draft; G. Hussain: Validation, Software; I. Inayat: Resources, Software; A. Khan: Conceptualization, Supervision, Visualization, Methodology; G. Zaman: Data curation, Formal analysis; A. Raezah: Writing–review and editing, Investigation.

Acknowledgments

The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University for funding this work through Research Group Project under grant number RGP 1/78/45.

Conflict of interest

The authors declare no conflict of interest.