Dynamics and simulations of stochastic COVID-19 epidemic model using Legendre spectral collocation method

1.   Introduction

Everywhere in the world there have been disasters like; tsunamis, earthquakes, floods, droughts and pandemics since the beginning of the time. These disasters are likely to happen at different times and in different areas. For instance, at various points in history, humans have been afflicted with a variety of infectious diseases, including Ebola virus, HIV, dengue virus, bird flu, malaria, tuberculosis, influenza, diarrhea, hepatitis C and B, rubella disease and many others. These highly contagious diseases afflicted both humans and animals. In addition to this, these contagious diseases also spread among plants such as; Moko diseases, Pine wilt, orange rust, sugar cane, Dutchelm and Karnal bunt etc. diseases ^[1]. These widespread diseases not only left the bad impact on world economic situation, but they also results in the death of a huge number of people.

Mathematical models play a key role in the study of infectious diseases. Due to COVID-19 pandemic recently, a lot of work has been done and some good results have been achieved so for. The research on contagious diseases is vital for the better understanding of spread of these diseases and its behavior on biological entities. Human population is rapidly expanding, which has a negative impact on the nature. The spreading of these diseases are far more different than that of the population; as some of the diseases spread more rapidly and some of them spread with a slow rate. The first human case on COVID-19 was reported on 31st December, 2019 in Wuhan, China. First symptoms of COVID-19 were thought to be pneumonia. The infected patient did not recover after being vaccinated and treated for pneumonia, and the medication was ineffective ^[2]. Second, it was discovered that this virus spread quickly from person to person. Moreover, the infected patients were not restricted to Wuhan but this infectious disease further rapidly spread geographically within China ^[3,4]. The prevalence of COVID-19 virus was not restricted to Wuhan, it spread to other parts of the globe as well. Thousands of people have died as this deadly virus has spread rapidly throughout China and the rest of the countries, for instance North America, Italy, France and Germany etc. This infectious virus also hits Asia, particularly India and Pakistan. It was identified in Pakistan on February 27th, 2020 ^[5]. Studies on other contagious diseases like HIV, chickenpox, hepatitis E, dengue fever, rubella disease, measles and tuberculosis are studied in ^[6,7,8,9]. For comprehensive description of the vaccination game for mathematical epidemiology by analytical and multi-agent simulation approaches, we refer the reader to ^[10].

Recently re-emerging disease Monkeypox, the Ebola virus disease (EVD) and the influence of quadratic Lévy noise on the dynamic of an SIC contagious was studied in detail in ^[11,12,13]. For the understanding of the dynamics of the these epidemic diseases, they used some nonstandard analytical methods and find that the noise intensities are correlated with some parameters. They use real data to confirmed their theoretical justifications. Stochastic delay differential equations modeling these pandemic diseases with a time-periodic delay is analyzed using the spectral element approach in ^[14]. A comprehensive study consist of numerical solutions to stochastic differential equations can be found in ^{[15,16,17,18,19,20,21,22,23,24,25]}. Some recent work on the analysis of these epidemic models can also be found in ^{[26,27,28,29]}. As these epidemic models play an important role in helping to understand the spread of diseases and to plan control policies for diseases. Some researcher used statistical analysis of diseases data using epidemic models, which may not be very helpful for the above purpose as data may not be correct. In addition to this, the mathematical descriptions of these epidemic models are coupled and nonlinear in nature and it is very hard or in some cases it is not possible to find the exact solution. For this reason, one must use some numerical techniques for the approximate solution of these models. To the best of our knowledge, most of the authors used low order methods to find the numerical solution, which may not be helpful to study the long time behavior of these models. To this end, we use a higher order numerical scheme based on spectral methods to study the effects of different parameters on the COVID-19 epidemic SEIR model. Spectral methods are well known for their exponential rate of convergence and one can get a very accurate solution for relatively very small collocations points.

The remaining structure of the paper is: Mathematical model is presented in section 2, followed by spectral method in section 3. The stability analysis along with some main results are provided in section 4, followed by the brief numerical simulations in section 5. Conclusions are drawn in section 6.

2.   Mathematical model

Let $N(t)$ denotes the overall population, which is further classified into four subclasses. Susceptible, exposed, infected, and recovered individuals $S(t)$, $E(t)$, $I(t)$ and $R(t)$ respectively. According to proposed COVID-19 model, the recruitment rate denoted by $\Lambda$ of susceptible individual, $\eta$ the death rate, $\xi$ denote the interaction rate which is equal to $\xi SI/1+aI^2(t)$, where $\kappa$ denotes the rate of incubation period of the exposed class and enter into the infected class, $\delta$ represents the recovery rate. In mathematical modeling of epidemic diseases, incident rate play an important rule. In many classical models the bilinear rate is most frequently used, especially when the susceptible population is huge in number. For this reason we also choose the incident rate to be bilinear. Another important aspect of the stochastic models is the the persistence time in the mean of the stochastic model, which is the time during which an epidemic population survived. One can think about this in the context of COVID-19 pandemic. Researcher trying to figure out how long a population model persist, which is consider in detail in ^[30]. The deterministic COVID-19 model along with bilinear incidence rate is given by:

where the initial values $S(0) = S^*(t)\geq 0, E(0) = E^*(t)\geq 0, I(0) = I^*(t)\geq 0$ and $R(0) = R^*(t)\geq 0.$

The model given in Eq (2.1) in stochastic form is given below:

where $\epsilon$ is the intensity of $B(t)$, and $B(t)$ is a Wiener process. The aim of the present research work is to apply a reliable and stable numerical technique to the systems given in Eqs (2.1) and (2.2), along with theoretical justifications.

3.   Method description

Let $Q_{n}({\tau_a})$, denotes the Legendre polynomials of $n^{th}$ order. The function $g(\tau_a)$ given on the interval $\big[-1, 1\big]$ is approximated by

$X_i = g(\tau_{ai})$ indicates Legendre coefficients which are unknowns, $\tau_{ai}, i = 0, 1, ..., n$ are interpolating points $-1 = \tau_{a0} < \tau_{a1} < ... < \tau_{an} = 1$. The Legendre polynomial is given by:

The proposed numerical scheme for model Eq (2.2) is based on the usage of Lagendre-Gauss quadrature with a weight function and Legendre-Gauss Lobatto points $\{\tau_k\}_{k = 0}^{n}$. Now, to solve the proposed model Eq (2.2), we begin by integration from $[0, t]$ of Eq (2.2), to get:

where $S(0), E(0), I(0), R(0)$ are initial conditions for each class respectively. Since the proposed scheme is orthogonal on a standard interval $[-1, 1]$. Therefore, by using simple transformations $s = \frac{t}{2}(\vartheta +1)$, Eq (3.3) becomes:

The semi-discretised spectral system with weights functions of Eq (3.4) is given by:

where deterministic weighted function is:

While for stochastic terms the weighted function is given by:

Here the numerical solutions for each class using Eq (3.1), is given as follow:

where $S_{n}$, $E_{n}$, $I_{n}$, $R_{n}$ are unknown Legendre coefficients. Using the above approximation in a system Eq (3.5), one get:

for simplification, we use $\xi_k = \big[\frac{t}{2}(\vartheta_k +1)\big].$ Hence, Eq (3.7) consists of $(4N+ 4)$ unknowns with $4N$ nonlinear equations. Now, the initial conditions are given by:

therefore, Eq (3.7) along with Eq (3.8) forms a system of $(4N+4)$ non-linear equations. We first find these unknown values from Eqs (3.7) and (3.8), then using it in Eq (3.6) gives the solution to the system given in Eq (2.2).

4.   Stability analysis

In this section, we examine the stability of the proposed scheme for the COVID-19 systems given in Eqs (2.1) and (2.2). Keeping in mind that the system's essential reproduction number of the system in Eq (2.1) represented by $R_0$ calculated using next generation method:

Theorem 4.1. If $R_0 < 1,$ then system Eq (2.1) satisfy the disease free equilibrium, $\big(S^*(t), E^*(t), I^*(t), R^*(t)\big)\rightarrow (\frac{\Lambda}{\eta}, 0, 0, 0)$. If $R_0 > 1,$ then proposed system has stable endemic equilibrium $E^* = \big(S^*(t), E^*(t), I^*(t), R^*(t)\big)$, where

Proof. Consider $E^*$ be the endemic equilibrium of Eq (2.1), then $S^*(t), E^*(t), I^*(t), R^*(t)$ the equations in the system are satisfied:

Here we consider two cases, that is $I^*(t) = 0$ and $I^*(t) > 0$. First, if $I^*(t) = 0$, from the last equation of system Eq (4.1), we are simply get the recovery class becomes zero therefore, $R^*(t) = 0$, using this in second equation of proposed system Eq (4.1), the exposed class become zero as, $E^*(t) = 0$. Moreover, using the first equation of system Eq (4.1), we found the susceptible class as $S^*(t) = \Lambda/\eta$. Now we can get the complete disease free equilibrium, $\big(S^*(t), E^*(t), I^*(t), R^*(t)\big)\rightarrow (\frac{\Lambda}{\eta}, 0, 0, 0)$. If $I^*(t) > 0$, for such calculation we simply use Maple-13 and found suitable endemic equilibrium.

Lemma 4.2. Assume the given region $\Pi = \big\{\big(S(0), E(0), I(0), R(0)\big)\in R^4_+: S(t) > 0, E(t) > 0, I(t) > 0, > 0, R(t) > 0;$ $S(t)+E(t)+I(t)+R(t)\leq \frac{\Lambda}{\eta}\big\}$ is a positively invariant set for system Eq (2.1).

Proof. The observed total population are denoted by $N(t)$ by Eq (2.1), we obtain

This implies that

As $S(t)+E(t)+I(t)+R(t)\leq \frac{\Lambda}{\eta}$, thus the region $\Pi$ is positively invariant.

Lemma 4.3. Each given initial values $\big\{S(0), E(0), I(0), R(0)\in R^4_+\big\}$, then solution $\big\{S(t), E(t), I(t), R(t)\in R^4_+ \big\}$ of Eq (2.2) has the properties below:

for each class. For the proof one can see Meng et al. ^[31]

Definition 4.4. Given stochastic system Eq (2.2), the infected class $I(t)$ is called extinctive if $\lim_{t\rightarrow \infty}I(t) = 0$.

The basic reproduction number for the stochastic system Eq (2.2) denoted by $\bar R_0$ is given by:

Theorem 4.5. If $\epsilon^2 > \max \big\{\frac{\xi\eta}{\Lambda}, \frac{\xi^2}{2(\eta+\delta)(\eta+\kappa)}\big\}$ or $\epsilon^2 < \frac{\xi\eta}{\Lambda}$ with $\bar R_0 < 1$, then the infected class of stochastic system Eq (2.2), exponentially approaches to zero.

Proof.} Let us assume $\big\{S(t), E(t), I(t), R(t)\big\}$ be the solution of stochastic system Eq (2.2) and satisfied the initial values $\big\{S(0), E(0), I(0), R(0)\big\}$. Then according to the It$\hat{o}$ formula, we have the following form:

Integral from $0$ to $t$ of system Eq (4.2) we have

Now we discuss the following two different cases, if $\epsilon^2 > \frac{\xi\eta}{\Lambda}$, then

Divide Eq (4.4) by a positive $t$, then

After this taking the $\lim_{t\rightarrow \infty}$ and by Lemma 4.3, then Eq (4.5) converted to:

Therefore when $\epsilon^2 > \frac{\xi^2}{2(\eta+\delta)(\eta+\kappa)}$, which show that $\lim_{t\rightarrow \infty}I(t) = 0.$

The second case if $\epsilon^2 < \frac{\xi\eta}{\Lambda}$, then from Eq (4.3), we have

Dividing Eq (4.6) both sides by $t > 0$, then we obtain

Again by $\lim_{t\rightarrow \infty}$ and use the Lemma 4.3, Eq (4.7) becomes:

For $\bar R_0 < 1$, then

This shows that $\lim_{t\rightarrow \infty}I(t) = 0.$

Theorem 4.6. The infected population of humans and birds, therefore $I(t)$ respectfully are presents in Eq (2.2). For proof, we refer the reader ^[30].

5.   Numerical discussions

In this section, we discuss the stochastic model numerically using the Legendre spectral collocation method (LSCM) for both deterministic and stochastic systems. Figures 1–8 demonstrate the approximate results to the proposed models. Using LSCM, we can precisely describe the behavior of the fundamental reproduction number $R_0$. We established several parameters values for the current systems in order to demonstrate the simple and effectiveness of the suggested strategy for stochastic COVID-19. Our initial conditions are all set to $1$. For the sake of simplification, we fix the noise intensity $\epsilon$ from $0$ to $0.6$ in all our simulations. This is also due to fact that noise intensity dramatically change the dynamics of stochastic models, which are used for mathematical description of epidemic diseases. If the intensity is large, it can shift the bifurcation points and produce noise-induced transitions and making it very hard to find the bifurcation boundaries.

Figure 1 is obtained using the parameters values $\Lambda = 4; a = 1; \eta = 1; \xi = 1; \kappa = 0.6; \delta = 1; \epsilon = 0$, in deterministic system given in Eq (2.1). For above parameter values the reproduction numbers $R_0 < 1$, thus using Theorem 4.1, we can surly see in Figure 1, that the proposed systems consists disease free equilibrium $E_0$, that is $\big(S^*(t), E^*(t), I^*(t), R^*(t)\big)\rightarrow \big(\frac{\Lambda}{\eta}, 0, 0, 0\big)$, where the susceptible class is $\Lambda/\eta = 4$. In Figure 2, the parameter values are $\Lambda = 4; a = 1; \eta = 1; \xi = 1.7; \kappa = 0.6; \delta = 1; \epsilon = 0$, while he reproduction number for these parameters is $R_0 = 1.275 > 1$. Now again applying Theorem 4.1, the system Eq (2.1), has stable endemic equilibrium $E^*$. We clearly see that all the classes are non-zero.

For Parameter values $\Lambda = 4; a = 1; \eta = 1; \xi = 1; \kappa = 0.6; \delta = 1; \epsilon = 0.6$, we obtain Figure 3. Using this the stochastic reproduction number becomes $\bar{R}_0 < 1$. Now by applying the Theorem 4.6, we can clearly observe from Figure 3, that the systems Eq (2.2) has stable disease free equilibrium, where the susceptible class is $\Lambda/\eta = 4$ and all other classes are zeros. Using same parameter values given in Figure 4 and only choose $\xi = 1.7$, for the solution to the system Eq (2.2). Simulation present that $\bar{R}_0 > 1$, again by using Theorem 4.6, we see that given systems has stable endemic equilibrium and all classes of the system Eq (2.2) are non zero. Taking the parameters values $\Lambda = 4; a = 1; \eta = 1; \xi = 1; \kappa = 0.6; \delta = 1; \epsilon = 0.6$, we get Figure 5. We compare the numerical solution for each class of the systems Eqs (2.1) and (2.2). By applying the Theorem 4.6 again to the parameter values $R_0 < 1$ and $\bar{R}_0 < 1$, we can see from Figure 5 that both systems have a stable disease-free equilibrium $E_0$.

Figure 5 shows a comparison of the numerical solutions of deterministic and stochastic COVID-19 systems using the parameter values listed above. Using the above parameter values the simple calculation appears that and both the reproduction numbers are less then 1 clearly observe in figure 6. A stable endemic equilibrium $E_0$ may be found utilizing the preceding results. For parameter values, $\Lambda = 4; a = 1; \eta = 1; \xi = 1; \kappa = 0.6; \delta = 1; \epsilon = 0.6$ in Figure 7, we plot both humans and bird classes and also we compare the numerical solution for class to the systems Eqs (2.1) and (2.2) of COVID-19 model. With these parameters, both reproduction numbers are greater than 1, the system is in a stable endemic equilibrium, as seen in figure 7. We draw the comparison of the numerical solutions to Eqs (2.1) and (2.2), using same parameter values given in Figure 7. Since both reproduction numbers in these models are greater than 1, they have a stable infectious equilibrium $E^*$ as shown in Figure 8.

6.   Conclusions

A mathematical concept called differentiation and integration is used to model the dynamic behavior of a disease, like COVID-19. These models can be used to forecast how a disease will propagate in a particular population like how many people will get sick by splitting the overall populations into several subgroups for this assignment. Two equilibrium points DFE and EE are derived for the proposed steady state models. The fundamental reproduction number $R_0$ is determined using the next generation approach. We find through stability assessment that the provided model is locally asymptotically stable when $R_0 < 1$. A number of numerical simulations are performed to confirm the theoretical justifications.

Acknowledgments

This work was financially supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [Project No. GRANT784].

Conflict of interest

The authors declare that they have no competing interests.