Extinction and stationary distribution of a stochastic COVID-19 epidemic model with time-delay

1. Introduction

The outspread of infectious diseases like COVID-19 have been reported employing mathematical models such as stochastic and deterministic. Almost all models are the offshoots of a classical SIR model by Kermack Mckendrick [1]. SIR model is sub-divided in three groups such as susceptible S, infected I and recovered R population. The primary framework of the disease in a population is associated with the rate of incidence. It is therefore related with the mean of secondary cases evolved by an infected individual in the susceptible population. Many descendent models have been employed after Kermack McKendrick model [2,3]. The variations in our social and environmental differences in our daily life are the justifications of establishing stochastic integration in such models. Stochastic noise of a model reshapes the solution behavior of correlated deterministic system and also changes the threshold level of a system for an epidemic to occur. The noise induction affects the dynamics of the population [4]. An epidemiological infection with a source of noise with memory was employed with a dynamical model [5]. The epidemic dynamics model was analyzed by employing pulse noise model. Elsewhere threshold variation is described to examine the stochastic SIR model [6]. An increasing attention has been noticed for the analysis and control of COVID-19, also for vaccination and treatment policies. The association of vaccination or other treatment strategies and their relation with the transmission of a disease has been a hot topic for theoretical and applied analysis [[7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21]]. The disease transmission modeling in a population where vaccination is under effect, the main issue is the inefficiency of the vaccine in a given population. There is a possibility of low efficacy such as partial induction of immunization. Considering SIR-type disease such as COVID-19 during the vaccination program is in effect, the total population is divided into four classes i.e susceptible, infected, vaccinated and removed represented as S, I, V, and R respectively. (see Table 1 )

The basic reproduction number of this model is $R_0=\frac{\beta }{\mu +\lambda }$. In this model, it is assumed that during a given time, a fraction of the susceptible class is vaccinated. The vaccination may or may not immunize the individual, so the model includes ρ factor which is 0 < ρ < 1, where ρ = 0 refers to the effectiveness of the vaccine and ρ = 1 refers to the non-effectiveness of the vaccine. We also assume that the effect of vaccine is lost at some rate of proportion θ. Therefore, since the immunity is long lasting hence a fraction λ of infective goes to the removed class. It is also supposed that birth rate of the population takes place at a constant rate m of death and the neonates move into susceptible class. Thus, the overall population is constant and variables are normalized. In this study, we observed as a case where the vaccine is effective (θ = 0). This model can be further modified according to the environment of the system. The stochastic version of the above model is presented as:

where W ₁(t), W ₂(t), W ₃(t) and W ₄(t) stand for the independent Brownian motions. ${\sigma }_1^2$, ${\sigma }_2^2,{\sigma }_3^2$ and ${\sigma }_4^2$ are white noises, with ICs:

This paper is presented as follows: The existence and uniqueness have been carried out in Section 2. In Section 3, Extinction analysis of the underlying model is investigated. In Section 4, the existence of ergodic stationary distribution. Numerical simulations by using first order stochastic Runge-Kutta scheme is demonstrated in Section 5. In the last section i.e 6, concluding remarks are presented.

2. Existence and uniqueness

As the solution of SDE (2) has biological significance, it should be nonnegative [22]. Moreover, in order for a stochastic differential equation to have a unique global (i.e., no explosion in a finite time) solution for any given initial value, the coefficients of the equation are generally required to satisfy a linear growth condition and a local Lipschitz condition [23]. However, the coefficients of SDE (2) do not satisfy a linear growth condition, though they are locally Lipschitz continuous. In this section, we will use a method similar to the proof of [24,25], to prove that the solution of SDE (2) is nonnegative and global.

Theorem 1

System (2) has a unique positive solution (S(t), I(t), V(t), R(t)) on t ≥ − τ, and the solution will remain in ${\mathbb{R}}_{+}^4$ for the given initial condition (3) with probability one.

Proof 1

We define a C ² − function $V:{\mathbb{R}}_{+}^4\to {\mathbb{R}}_{+}$ as follows:

where k > 0 will be determined later on. By Ito's formula, we can obtain

where

Let $k=\frac{\mu -\rho \beta }{\beta }$, then we have

where $\mathcal{M}>0$. Hence,

Integrating (6) from 0 to ${\tau }_n\wedge \tilde{T}=\min \{{\tau }_n,\tilde{T}\}$ leads us

implies

According to (7), we get

limiting case leads us

contradiction arises hence τ _∞ = ∞, a.s.

3. Extinction

In this section, we will show that if the noise is sufficiently large, the solution to the associated SDE (2) will become extinct with probability 1 [[26], [27], [28]].

Lemma 1

Let$M={\{M_t\}}_{t\ge 0}$be a real-valued continuous local martingale vanishing at t = 0, and${⟨M,M⟩}_t$be the quadratic variation of M. Then

$\underset{t\to \infty }{\lim }{⟨M,M⟩}_t=\infty ,a.s.\Rightarrow \underset{t\to \infty }{\lim }\frac{M_t}{{⟨M,M⟩}_t}=0a.s$.,

and also.

$\underset{t\to \infty }{\lim }\sup \frac{{⟨M,M⟩}_t}{t}<\infty a.s.\Rightarrow \underset{t\to \infty }{\lim }\frac{M_t}{t}=0,a.s.$,

Lemma 2

Let (S(t), I(t), V(t), R(t)) be the solution of (2) with any $(S(0),I(0),V(0),R(0))\in {\mathbb{R}}_{+}^4$, then

$\underset{t\to \infty }{\lim }\frac{S(t)}{t}=0,\underset{t\to \infty }{\lim }\frac{I(t)}{t}=0,\underset{t\to \infty }{\lim }\frac{V(t)}{t}=0,\underset{\to \infty }{\lim }\frac{R(t)}{t}=0,a.s$.,

Furthermore, if$\mu >\frac{{\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2\vee {\sigma }_4^2}{2}$, then.

$\underset{t\to \infty }{\lim }\frac{{\int }_0^{t}S(s)d{W}_1(s)}{t}=0,\underset{t\to \infty }{\lim }\frac{{\int }_0^{t}I(s)d{W}_2(s)}{t}=0,\underset{t\to \infty }{\lim }\frac{{\int }_0^{t}V(s)d{W}_3(s)}{t}=0$,

$\underset{t\to \infty }{\lim }\frac{{\int }_0^{t}R(s)d{W}_4(s)}{t}=0$, a.s.

Theorem 2

If$R_0^s<1$and$\mu >\frac{{\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2\vee {\sigma }_4^2}{2},$then (2) obeys:

Proof 2

Let $\mathcal{U}(t)=\lambda (I+V)+(\lambda +\mu )R$, and applying Ito's formula leads us,

From model (2), we have

Integration gives us

where

Using Lemmas 1 and 2,

limit of (13) gives us,

Integrating leads us

where

Incorporating Lemmas 1 and 2.

$\underset{t\to \infty }{\lim }{\psi }_2(t)=0,$a.s.

Since $R_0^s<1$, limit of (15) leads us

implying $\underset{t\to \infty }{\lim }I(t)=0,\underset{t\to \infty }{\lim }V(t)=0,\underset{t\to \infty }{\lim }R(t)=0$ a.s.

From (14), we have $\underset{t\to \infty }{\lim }⟨S⟩=1,a.s$.

4. Stationary distribution

Herein, we construct a suitable stochastic Lyapunov function to study the existence of a unique ergodic stationary distribution [29,30] of the positive solutions to the system (2).

Consider

where $\widehat{\lambda }=\frac{\phi \rho \beta \mu }{(\lambda +\mu +\frac{{\sigma }_2^2}{2})},\widehat{\theta }=\frac{\phi \rho \beta \mu }{(\mu +\theta +\frac{{\sigma }_3^2}{2})},\widehat{\phi }=\mu +\phi +\frac{{\sigma }_1^2}{2}$ and $\widehat{\mu }=\mu +\frac{{\sigma }_4^2}{2}$.

Theorem 3

Assume that$R_0^s>1$and$\mu -\frac{{\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2\vee {\sigma }_4^2}{2}>0$, then for value$(S(0),E(0),I(0),Q(0))\in {\mathbb{R}}_{+}^4$, then (2) possess SD π(.).

Proof 3

To prove the theorem we take the help of two conditions in Lemma 1 of [24]. For this, we consider the diffusion matrix of model (2) as:

$\mathrm{\Lambda }=\left(\begin{array}{cccc}{v}_1^2{S}^2& 0& 0& 0\\ 0& v_2^2{I}^2& 0& 0\\ 0& 0& v_3^2{V}^2& 0\\ 0& 0& 0& v_4^2{R}^2\end{array}\right)$.

It is easy to show that Λ is positive definite, hence the first condition of Lemma 1 in Ref. [24] is satisfied.

Furthermore, consider ${\mathcal{C}}^2$-function $\mathcal{V}:{\mathbb{R}}_{+}^4\to \mathbb{R}$:

where $c_1=\frac{\phi \rho \beta \mu }{\lambda +\mu +\frac{{\sigma }_1^2}{2}},c_2=\frac{\phi \rho \beta \mu }{\mu +\theta +\frac{{\sigma }_3^2}{2}}and{c}_3=\frac{\phi \rho \beta \mu }{\mu +\frac{{\sigma }_4^2}{2}}$. Let

where $(S,I,V,R)\in (\frac{1}{n},n)\times (\frac{1}{n},n)\times (\frac{1}{n},n)\times (\frac{1}{n},n)$ and n > 1.

${\mathcal{V}}_1=-\ln S-c_1\ln I-c_2\ln V-c_3\ln R+\beta {\int }_t^{t+\tau }I(s-\tau )ds$,

${\mathcal{V}}_2=-\ln S+\beta {\int }_t^{t+\tau }I(s-\tau )ds$,

${\mathcal{V}}_3=-\ln V,{\mathcal{V}}_4=-\ln R$,

${\mathcal{V}}_5=\frac{1}{\xi +1}{(S+I+V+R)}^{\xi +1}$,

ξ > 1 is a constant satisfying

and $\mathcal{N}>0$, obeying

Where $\eta =\frac{\phi \rho \beta \mu }{\widehat{\lambda }\widehat{\theta }}-(\mu +\phi +\frac{{\sigma }_1^2}{2})-c_3(\mu +\frac{{\sigma }_4^2}{2})>0$,

and

Applying Itô’s formula to ${\mathcal{V}}_1$, we have

Similarly, we can get

where $\varpi =({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2\vee {\sigma }_4^2)$, Using (22), (23), (24), (25), (26), leads us

Let

a bounded closed set and ε > 0. In the set ${\mathbb{R}}_{+}^4\backslash \mathcal{D}$, consider