Dynamics of a stochastic COVID-19 epidemic model considering asymptomatic and isolated infected individuals

1.   Introduction

Since December $2019$, many cases of unexplained pneumonia with a history of exposure to the Southern Seafood Market of China have been reported by some hospitals in Wuhan, Hubei province of China. With the rapid growth of confirmed cases, all provinces in China have been affected by the coronavirus pneumonia. On 11 February $2020$, the World Health Organization (WHO) announced that coronavirus disease would be named as "COVID-19". The virus of COVID-19 spreads quite fast and severe cases can lead to dyspnea and even death. Due to the timely and effective measures taken by the government, the transmission of COVID-19 had been basically controlled in China until the end of April $2020$. More seriously, over $200$ countries around the world are affected by the COVID-19 virus ^[1]. The pandemic has a severe influence on the public health and economics in many countries, and it is straining health-care systems to their limits for developing countries such as India ^[2].

The COVID-19 pandemic has attracted worldwide attention, and the governments of most countries have taken strict epidemic prevention and control measures. Public Health and Social Measures (PHSM) suggested by WHO include that: (ⅰ) facial coverings and/or mask wearing; (ⅱ) adaptation or closure of schools and businesses; (ⅲ) limits and restrictions on public and private gatherings; (ⅳ) restrictions on domestic movement and international travel. Moreover, a large proportion of susceptible population around the world has been vaccinated against COVID-19. According to the data of WHO, a total of $10,227,670,521$ vaccine doses had been administered as of 13 February 2022. Nevertheless, due to the high variation and infectivity of coronavirus, the pandemic is still serious. By the middle of April $2021$, the global cumulative cases had ballooned to $130$ million with over $2.96$ million deaths. But ten months later, by the middle of February $2022$, there had been more than $412$ million confirmed cases of COVID-19, including over $5.82$ million deaths reported to WHO ^[3].

Epidemic model is a critical tool to explore the spread and control strategy of infectious diseases ^{[4,5,6,7,8,9]}. Since the pandemic of coronavirus disease, a large number of epidemic models based on COVID-19 have been proposed and investigated ^{[10,11,12,13,14,15,16]}. For instance, Kucharski et al. ^[17] built a model based on all confirmed cases until 5 March $2020$ in Wuhan. Wu et al. ^[18] investigated a SEIR model with the basic reproduction number 2.68, and simulated the epidemic situation in Wuhan from 31 December $2019$ to 28 January $2020$. Din et al. ^[19] studied the spread and control of the novel corona virus (COVID-19) in China through a SIQ model and a control problem using two control measures. Meiksin ^[20] used the SEIR model to constrain the role of contaminated fomites in spreading an epidemic, and applied the approach to the spread of COVID-19 in the UK during the lockdown in the first half of $2020$. Mello et al. ^[21] provided a comprehensive review focused on Covid-19 using the Ising-model and percolation theory, and their main results included that: (ⅰ) the temporal evolution of the accumulative number of infections and fatalities follow a logistic function; (ⅱ) the percolation probability can be associated with the probability of a person being infected with Covid-19. For more details the reader can refer to the literatures ^[21,22] and the references therein.

Recently, Sarkar et al. ^[23] proposed and analyzed a high-dimensional COVID-19 epidemic model, in which the total population is divided into six categories including asymptomatic and isolated infected individuals. In their work, the authors calculated the basic reproduction number, and simulated the model to predict the dynamics of COVID-19 in $17$ provinces of India and the overall India. The model in reference ^[23] was a deterministic epidemic model. However, in the real world, biological populations exist inevitably in a noisy world of random variation. Environmental randomness can have important effects on the growth rate, environmental capacity as well as other parameters in biological models ^[24]. Furthermore, the probability for an individual to be infected by others depends on their own characteristic, such as age, nutritional status, sex and so on ^[25]. Therefore the transmission of disease is inevitably affected by environmental disturbance, and stochastic models are more realistic considering different body constitutions and other factors ^{[26,27,28,29]}.

In this paper, we extend the COVID-19 epidemic model in reference ^[23] from a deterministic framework to a stochastic epidemic model. The paper is organized as follows. Section 2 describes the formulation of the stochastic epidemic model. In Section 3, we provide the mathematical analysis of the stochastic epidemic model. Firstly, the existence and uniqueness of positive solution is verified. Then the parameter condition for the extinction of the coronavirus disease (COVID-19) is obtained. Finally, we prove that there is a unique ergodic stationary distribution for the stochastic epidemic model under certain conditions. Epidemiologically, the existence of stationary distribution implies that the disease will persist almost surely in the time mean sense. In Section 4, we conduct numerical simulation and also compare the simulation with the real data in India. A brief conclusion is given in the last section.

2.   Model formulation

Sarkar et al. ^[23] proposed a COVID-19 model, where the total population is divided into six categories: susceptible individuals $(S)$, quarantined susceptible individuals $(S_q)$, asymptomatic infectious individuals $(A)$, infected or infectious with symptoms $(I)$, isolated infected individuals $(I_q)$, recovered compartment $(R)$ (no more infectious). And the model is expressed as follows

Here the total population $N(t) = S(t)+S_q(t)+A(t)+I(t)+I_q(t)+R(t)$. In the model: (ⅰ) $\Lambda_s$ represents the net inflow rate of susceptible individuals, and $\delta$ denotes the natural mortality rate; (ⅱ) $\beta_s$ is the probability of disease transmission, and $\rho_s$ represents the quarantined rate of susceptible individuals; (ⅲ) $m_s$ denotes the rate at which quarantined susceptible individuals are released into uninfected class, and $\varepsilon_s$ is the daily contact rate per unit of time; (ⅳ) $\gamma_a$ is the probability rate at which asymptomatic individuals develops clinically symptoms, and $\gamma_i$ is the probability rate at which infected individuals become isolated; (ⅴ) $\xi_a$, $\xi_i$ and $\xi_q$ denote recovery rate of asymptomatic infected individuals, infected individuals, and isolated infected individuals, respectively. All of the parameters are nonnegative. One can refer to reference ^[23] for more detailed biological interpretations of model (2.1).

The spread of diseases in the real world is inevitably subject to random environmental perturbation. There are different possible approaches to introduce random effects in the epidemic models. Due to continuous changes in nature, environmental disturbance always fluctuates around some average values, and it is usually assumed to be types of white noise ^[30]. In this work, motivated by the approaches in references ^[1,24,31], we introduce to system (2.1) Gaussian white noise which is directly proportional to $S(t)$, $S_q(t)$, $A(t)$, $I(t)$, $I_q(t)$ and $R(t)$, and obtain

where $B_i(t)$ $(i = 1, \, 2, \, 3, \, 4, \, 5, \, 6)$ are mutually independent standard Brownian motions, and $\sigma_i$ $(i = 1, \, 2, \, 3, \, 4, \, 5, \, 6)$ are noise intensities. Throughout the paper, let $(\Omega, \:\mathcal{F}, \:\{\mathcal{F}\}_{t\geq0}, \:\mathbb{P})$ be a complete probability space with a filtration ${\{\mathcal{F}_{t}\}_{t\geq0}}$ satisfying the usual conditions (i.e. it is increasing and right continuous while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets), and $B_{i}(t)$ is defined on this complete probability space. We also recall the following notations: $\mathbb{R}_{+}^{d} = \{(x_1, ..., x_d):x_i > 0, 1\leq i\leq d\}$, $a\wedge b = \min\{a, b\}$, $a\vee b = \max\{a, b\}$, and $\langle f\rangle = \frac{1}{t}\int_0^t f(r)dr$.

3.   Mathematical analysis of the model (2.2)

In this section, we aim to theoretically analyze the stochastic epidemic model (2.2). To begin with, we introduce some basic definitions ^[30]. Let $X(t)$ be a regular time-homogeneous Markov process in $\mathbb{R}^{d}$ described by the stochastic differential equation

The diffusion matrix is defined as $A(x) = (a_{ij}(x)), a_{ij}(x) = g^i(x)g^j(x)$, and the diffusion operator $L$ is defined by

If $L$ acts on a function $V\in C^2(\mathbb{R}^{d}, \mathbb{R}_+)$, then

where

With the diffusion operator $L$, Itô's formula is written as

3.1.   The existence and uniqueness of positive solution

In this subsection, we prove the existence and uniqueness of positive solution for system (2.2), which ensures that the numbers of all counterparts are positive almost surely.

Theorem 3.1. For any initial value $(S(0), S_q(0), A(0), I(0), I_q(0), R(0))\in\mathbb{R}^{6}_{+}$, there is a unique positive solution $(S(t), S_q(t), A(t), I(t), I_q(t), R(t))$ of system (2.2) on $t\geq0$ and the solution will remain in $\mathbb{R}^{6}_{+}$ with probability one.

Proof. Since the coefficients of system (2.2) are locally Lipschitz continuous, then for any initial value $(S(0), S_q(0), A(0), I(0), I_q(0), R(0))\in\mathbb{R}^{6}_{+}$, there is a unique positive solution $(S(t), S_q(t), A(t), I(t),$ $I_q(t), R(t))$ on $t\in[0, \tau_e)$, where $\tau_e$ represents the explosion time ^[30]. To show this solution is global, we only need to prove that $\tau_e = \infty$ a.s. To begin with, let $k_0 > 0$ be sufficiently large such that $S(0), S_q(0), A(0), I(0), I_q(0)$ and $R(0)$ all lie in the interval $[\frac{1}{k_0}, k_0]$. For each integer $k\geq k_0$, define the stopping time

Here we set $\inf\emptyset = \infty$ ($\emptyset$ denotes the empty set). It is apparent that $\tau_k$ is increasing as $k\rightarrow \infty$. Let $\tau_\infty = \lim\limits_{k\to \infty}\tau_k$, then $\tau_\infty\leq\tau_e$ a.s. Thus, $\tau_\infty = \infty$ a.s. implies $\tau_e = \infty$ a.s. and $(S(t), S_q(t), A(t), I(t),$ $I_q(t), R(t))\in\mathbb{R}^{6}_{+}$ for all $t\geq0$ almost surely. Now we state that $\tau_\infty = \infty$. If this assertion is false, then there exist constants $T > 0$ and $\epsilon\in(0, 1)$ such that $\mathbb{P}\{\tau_\infty\leq T\} > \epsilon$. Hence there exists an integer $k_1\geq k_0$ such that

Next we define a $C^2$-function $H$: $\mathbb{R}^{6}_{+}\rightarrow \mathbb{R}_{+}$ by

Apply Itô's formula to $H$, and we obtain

where

and here $K$ is a positive constant. Therefore,

Let $k\geq k_1$. Integrating both sides of the above inequality from $0$ to $\tau_k\wedge T$ and taking the expectation yield

Set $\Omega_k = \{\omega\in\Omega:\tau_k = \tau_k(\omega)\leq T\}$ for $k\geq k_1$, then we have $P(\Omega_k)\geq\epsilon$ by $(3.3)$. Notice that for every $\omega\in\Omega_k$, there exists $S(\tau_k, \omega)$, $S_q(\tau_k, \omega)$, $A(\tau_k, \omega)$, $I(\tau_k, \omega)$, $I_q(\tau_k, \omega)$ or $R(\tau_k, \omega)$ equals either $k$ or $\frac{1}{k}$. Thus,

By virtue of (3.6) and (3.7), one has

and here $I_{\Omega_k(\omega)}$ is the indicator function of $\Omega_k$. Let $k \rightarrow \infty$, then

which is a contradiction. So it is obvious that $\tau_\infty = \infty,$ a.s. The proof is completed.

3.2.   Extinction

In this subsection, we will explore the parameter conditions for the disease extinction, which is of great help to disease control. We first present the following lemma.

Lemma 3.2. For any initial value $(S(0), S_q(0), A(0), I(0), I_q(0), R(0))\in \mathbb{R}^6_+$, the solution $(S(t), \\S_q(t), A(t), I(t), I_q(t), R(t))$ of system (2.2) satisfies

Furthermore, if $\xi_a+\delta > \frac{1}{2}(\sigma_1^2\vee\sigma_2^2\vee\sigma_3^2\vee\sigma_4^2\vee\sigma_5^2\vee\sigma_6^2)$ holds, then

The proof is similar to Lemma 2.1 and 2.2 in reference ^[32], so we omit it here.

Theorem 3.3. Set $\xi_a+\delta > \frac{1}{2}(\sigma_1^2\vee\sigma_2^2\vee\sigma_3^2\vee\sigma_4^2\vee\sigma_5^2\vee\sigma_6^2)$. For any given initial value$(S(0), S_q(0), A(0), I(0), I_q(0), R(0))\in \mathbb{R}^6_+$, if

then

Furthermore,

Proof. Define a $C^2$-function $U$ by

By Itô's formula, we have

Integrating from $0$ to $t$ and dividing by $t$ on both sides of (3.13), we obtain

Making use of Lemma 3.2, we have

which implies that

Furthermore, one can obtain that

Then integrating from $0$ to $t$ and dividing by $t$ on both sides of the first equation in system (2.2) yield

Hence,

The proof is finished.

The result in Theorem 3.3 shows that the disease will go to extinction almost surely if $R^* < 1$. According to the formula of $R^*$, one can see that the noise intensities $\sigma_3$ and $\sigma_4$ play a key role in the condition $R^* < 1$. In other words, the value of $R^*$ decreases with increasing noise intensities conducted on the population of asymptomatic infectious and infected with symptoms. Epidemiologically, in the practical policies to combat the COVID-19 spread, more attention should be paid to the fluctuations on the number of asymptomatic infectious individuals and infected with symptoms. It also indicates that large environmental noises may help to bring about elimination of the disease.

Remark 1. It has been showed, in a large number of works on stochastic epidemic models, that random fluctuations can suppress disease outbreak. Li et al. ^[33] and Cai et al. ^[34] studied the dynamics of stochastic epidemic models by adding Gaussian white noise to the transmission coefficient $\beta$. They showed that there exists a crucial noise intensity $\sigma_*$, and the diseases tend to extinction almost surely when $\sigma > \sigma_*$. The difference is that there is only one noise intensity in their stochastic models. Recent works on stochastic COVID-19 epidemic models also reported that increasing the value of the noise intensity will decrease the infection ^[35,36]. Such observation coincides well with our result. Specifically, Khan et al. ^[35] investigated a stochastic four-dimensional COVID-19 epidemic model by adding randomness to transmission rates, with two noise intensities in the model. A. El Koufi et al. ^[36] studied a stochastic Covid-19 epidemic model for a population with five compartments. They introduced white noise in the same way as model (2.2) of the present work, with five noise intensities in their model. However, there was no theoretical analysis on the extinction and persistence of the stochastic model, and it was not figured out which noise intensities play a role ^[36].

3.3.   Stationary distribution

Since there is no endemic equilibrium for a stochastic system, we usually focus on the stationary distribution to study the persistence of the disease. First of all, we recall the following well-known lemma.

Lemma 3.4. (^[37]) The Markov process $X(t)$ described by Eq (3.1) has a unique ergodic stationary distribution $\pi(\cdot)$ if a bounded domain $D\subset \mathbb{R}^d$ with regular boundary $\Pi$ exists and

$(H.1)$ there exists a constant $W > 0$ satisfying $\Sigma_{i, j = 1}^{d}a_{ij}(x)\xi_i\xi_j\geq W|\xi|^2$, $x\in D$, $\xi\in\mathbb{R}^{d}$.

$(H.2)$ there is a $C^2$-function $V\geq0$ such that $LV$ is negative for any $\mathbb{R}^{d}\setminus D$. Then

for all $x\in \mathbb{R}^{d}$, where $f(\cdot)$ is a function integrable with respect to the measure $\pi$.

Theorem 3.5. For any initial value $(S(0), S_q(0), A(0), I(0), I_q(0), R(0))\in \mathbb{R}^{6}_{+}$. If

then system (2.2) has a unique stationary distribution and it has ergodic property.

Proof. Firstly, we construct the following functions

where $\theta$ and $a_i\ (i = 1, 2, 3)$ are positive constants satisfying

Moreover, define

where $B > 0$ satisfies the following condition

with

and $\Gamma$ will be determined later. It is clear that

where $U_n = (\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)$. Since $V(S, S_q, A, I, I_q, R)$ is a continuous function, then $V(S, S_q, A, I, I_q, R)$ has a minimum point $(S^*, S^*_q, A^*, I^*, I^*_q, R^*)$ in the interior of $\mathbb{R}^{6}_{+}$. Next, we define a nonnegative $C^2$-function $F: \mathbb{R}^6_+ \rightarrow \mathbb{R}_+$ as follows

By Itô's formula, we have

Similarly,

where

Moreover,

Combining (3.16)–(3.18), we have

Furthermore, define the following set

where $\varepsilon$ is a sufficiently small positive constant such that

Thus,

and

Now we will verify that $LV(S, S_q, A, I, I_q, R)\leq-1$ on $\mathbb{R}^{6}_{+}\setminus D$.

Case 1. If $(S, S_q, A, I, I_q, R)\in D_1$, according to (3.20)

Case 2. If $(S, S_q, A, I, I_q, R)\in D_2$, then by inequality (3.21)

Case 3. If $(S, S_q, A, I, I_q, R)\in D_3$, inequality (3.22) implies that

Case 4. If $(S, S_q, A, I, I_q, R)\in D_4$, by (3.23)

Case 5. If $(S, S_q, A, I, I_q, R)\in D_5$, by (3.24)

Case 6. If $(S, S_q, A, I, I_q, R)\in D_6$, from (3.25)

Case 7. If $(S, S_q, A, I, I_q, R)\in D_7$, by (3.26)

Case 8. If $(S, S_q, A, I, I_q, R)\in D_8$, it then from (3.27) that

Case 9. If $(S, S_q, A, I, I_q, R)\in D_9$, then from (3.27)

Case 10. If $(S, S_q, A, I, I_q, R)\in D_{10}$, by inequality (3.28)

Case 11. If $(S, S_q, A, I, I_q, R)\in D_{11}$, from (3.29)

Case 12. If $(S, S_q, A, I, I_q, R)\in D_{12}$, according to (3.30)

Therefore, we have proven that $LV\leq-1$ for all $(S, S_q, A, I, I_q, R)\in\mathbb{R}^{6}_{+}\setminus D$. In other word, the condition $(H.2)$ in Lemma 2 is verified. In addition, the diffusion matrix of model $(2.2)$ is

Set $W = \min_{(S, S_q, A, I, I_q, R)\in D}\{\sigma_1^2S^2, \sigma_2^2S_q^2, \sigma_3^2A^2, \sigma_4^2I^2, \sigma_5^2I_q^2, \sigma_6^2R^2\}$, then

for $(S, S_q, A, I, I_q, R)\in D$, $\xi\in\mathbb{R}^{6}$. That is, the condition $(H.1)$ in Lemma 3.4 is satisfied. This completes the proof.

According to Theorem 3.5, there is a unique ergodic stationary distribution for system (2.2) if $R_0^S > 1$. From the perspective of biology, the existence of stationary distribution implies that the coronavirus disease (COVID-19) will persist almost surely in the time mean sense. Epidemiologically, it means the stochastically persistence of the disease in the long term. Moreover, let $\sigma_i = 0$ $(i = 1, \, 2, \, 3, \, 4, \, 5, \, 6)$, then $R_0^S$ is consistent with the basic reproduction number $R_0$ of the deterministic model (2.1), which indicates that we generalize the results of system (2.1). It should be mentioned again that $R_0^S$ decreases with the increase of one of the noise intensities $\sigma_1, \ \sigma_3$ and $\sigma_4$, and these three noise intensities play a crucial role in the condition $R_0^S > 1$. This observation coincides with the result that noise can suppress the disease outbreak ^[34].

4.   Numerical simulations

In this section, to illustrate our theoretical results more intuitively, we conduct numerical simulations to stochastic system (2.2) by making use of Milstein's high order method ^[38]. Firstly we find the data from 3 July to 16 September 2020 and 8 March to 7 April 2021 in India, as is shown in Figures 1(a) and 2(a), respectively. One can see that the number of the infectious is showing an upward trend in general. In Figures 1(b) and 2(b), we fit the model with actual data to illustrate the reliability of the stochastic model, in which the selection of parameter values is referred to reference ^[23]. In addition, the long-term behavior of the infected cases predicted by the stochastic model is shown in Figure 3, where the black solid line represents the actual data of existing cases displayed in Figure 1(a). Moreover, we also find the data of existing cases in India from 11 July to 4 October 2021 and fit the model with actual data, which is shown in Figure 4(a), (b). Figure 5 shows the simulation for the actual data of existing cases in India from October $2021$ to January $2022$. One can see that the number of existing cases tends to be relatively stable.

In order to verify the theoretical findings of our stochastic COVID-19 model, two more simulations are presented in the following examples.

Example 4.1. Choose $\Lambda_s = 0.5, \delta = 0.03, \beta_s = 0.07, \rho_s = 0.9, m_s = 0.05, \xi_q = 0.5, \varepsilon_s = 0.85,$ $\gamma_a = 0.9, \xi_a = 0.5, \xi_i = 0.05, \gamma_i = 0.04, \sigma_1 = 0.01, \sigma_2 = 0.3,$ $\sigma_3 = 0.2, \sigma_4 = 0.2, \sigma_5 = 0.2, \sigma_6 = 0.1$, and the initial values $(S(0), S_q(0), A(0), I(0), I_q(0), R(0)) = (0.5, 0.5, 0.5, 0.5, 0.5, 0.3)$. Then it is easy to calculate that $R^* = 0.9454 < 1$, that is, the condition of Theorem 3.2 is satisfied. Numerical simulation is shown in Figure 6(a), from which one can see that the disease eventually goes to extinction.

Example 4.2. Choose $\Lambda_s = 0.5, \delta = 0.11, \beta_s = 0.6, \rho_s = 0.3, m_s = 0.05, \xi_q = 0.33, \varepsilon_s = 9.1,$ $\gamma_a = 0.02, \xi_a = 0.03, \xi_i = 0.05, \gamma_i = 0.02, \sigma_1 = 0.02, \sigma_2 = 0.03,$ $\sigma_3 = 0.02, \sigma_4 = 0.02, \sigma_5 = 0.02, \sigma_6 = 0.01$, and the initial values $(S(0), S_q(0), A(0), I(0), I_q(0), R(0)) = (0.5, 0.5, 0.5, 0.5, 0.5, 0.3)$. Then one can compute that $R_0^S = 2.6431 > 1$. Therefore, the condition of Theorem 3.5 is satisfied and system (2.2) has a unique ergodic stationary distribution. The sample path for the solution is shown in Figure 6(b), from which one can clearly see that all the compartments coexist and the disease will be persistent in the long term. The corresponding density functions are presented in Figure 7, which is a simulation of stationary distribution.

5.   Conclusions

In the paper, we investigate a six-dimensional stochastic COVID-19 epidemic model. Compared with the general epidemic models, it includes the compartment of asymptomatic and isolated infection and is more in line with the characteristics of COVID-19. As for the stochastic model, we prove the existence and uniqueness of the positive solution. When we analyzing biological models, the first mission is to ensure that all counterparts of the population are positive. The existence of global positive solution guarantees that the stochastic model is meaningful from the viewpoint of biology. Moreover, using the method of constructing Lyapunov function, the conditions on the extinction and stationary distribution of the disease are obtained. Epidemiologically, the disease COVID-19 will go to extinction if $R^* < 1$. The stochastic system will admit a unique ergodic stationary distribution if $R_0^S > 1$. The existence of stationary distribution indicates the stochastically persistence of the disease in the time mean sense. In addition, the values of $R^*$ and $R_0^S$ will decrease with increasing the noise intensities conducted on the population of asymptomatic infection and infected with symptoms. Thus, we should pay more attention to the fluctuations on these two compartments in the practical policies to combat the COVID-19 spread. However, it is difficult to get a threshold for the disease control, which will be a topic in future research.

Numerical simulation is provided to illustrate our theoretical result. We also compare the real data in India with the simulation of the stochastic model, and it shows that our model fits well to the real data of a relatively short period. Nevertheless, in view of the complicated trend of the disease spread in the real world, it is difficult for our model to simulate the data of a quite long period. At the moment, the global pandemic is taking a turn for the better with increasing vaccination rates. We hope to formulate a more precise model considering the effect of vaccination in future.

Acknowledgments

The authors are grateful for the anonymous referees for their valuable comments and suggestions which led to the improvement of this paper.

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this manuscript.