Threshold of Stochastic SIRS Epidemic Model from Infectious to Susceptible Class with Saturated Incidence Rate Using Spectral Method

Abstract

Stochastic SIRS models play a key role in formulating and analyzing the transmission of infectious diseases. These models reflect the environmental changes of the diseases and their biological mechanisms. Therefore, it is very important to study the uniqueness and existence of the global positive solution to investigate the asymptotic properties of the model. In this article, we investigate the dynamics of the stochastic SIRS epidemic model with a saturated incidence rate. The effects of both deterministic and stochastic distribution from infectious to susceptible are analyzed. Our findings show that the occurrence of symmetry breaking as a function of the stochastic noise has a significant advantage over the deterministic one to prevent the spread of the infectious diseases. The larger stochastic noise will guarantee the control of epidemic diseases with symmetric Brownian motion. Periodic outbreaks and re-infection may occur due to the existence of feedback memory. It is shown that the endemic equilibrium is stable under some suitable initial conditions, taking advantage of the symmetry of the large amount of contact structure. A numerical method based on Legendre polynomials that converts the given stochastic SIRS model into a nonlinear algebraic system is used for the approximate solution. Finally, some numerical experiments are performed to verify the theoretical results and clearly show the sharpness of the obtained conditions and thresholds.

1. Introduction

Mathematical models became a vital tool after the fundamental work of McKendrick and Kermack [1]. This study seeks to learn about the causes and the control of infectious diseases. In ecological systems, random changes in the external environment constantly influence the human power. For example, models of the rate of disease transmission during an epidemic have been modified for a number of meteorological factors, such as the survival and spread of many viruses and bacteria in wet conditions and the amount of ultraviolet (UV) light [2,3,4]. Growth of the population and environmental impacts tend to combine with changes in the range of food supplies, which in turn is largely affected by unpredictable precipitation fluctuations [5]. In the literature, a very popular method is to randomly change modes in such a way that the Markov chain continuous-time values are separated with values in space in a finite state, which leads a change in basic parameters with the Markov chain state switching [6,7,8]. The dynamical system becomes regime switching via a mathematical equation in the form of Markov process, which is piecewise deterministic. Based on a biological system in terms of stochastic environmental circumstances, the epidemic models along with their deterministic primary parameters are not likely to be realistic. It is essential to examine the disease effect of irregular mutations in natural psychological states on the spread of the infection in the proposed population. Despite the significance of environmental noise, it has gained little consideration in the literature on infectious diseases. Initially, Gray, Greenhalgh and Mao et al. proposed an epidemic SIS model for Markovian switching [9]. For time discretization, some authors obtained reliable results taking the irregular switching of environmental regimes [10,11]. For the approximate solution of the stochastic differential models, some authors used LSCM [12,13]. Moreover, Legendre polynomials are also used for the numerical approximation of integrals, as well as differential equations [14,15,16,17,18]. The effect of Markov switching on a deterministic SIRS model is investigated in [19]. Due to the COVID-19 epidemic, SIRS stochastic models became more essential and many scholars have adopted them to study the transmission of COVID-19. For example, in [20], the authors used the classical deterministic SIR epidemic model to construct a more compatible model for COVID-19 transmission. Differential equations of fractional order for the mathematical modeling of COVID-19 were used by Khan et al., where they provided a detailed analysis by considering Saudi Arabia as a case study [21]. An approach based on piecewise differential and integral operators was developed by Atangana, Seda and Prolonger [22,23]. The dynamics of COVID-19 with the fuzzy Caputo approach are discussed in [24]. In all these works, the authors used deterministic approaches to formulate mathematical models for the transmission of COVID-19. However, these models have some limitations in describing the reality compared with the stochastic models, because the system parameters are not constant and fluctuate around an average value [25].

The connection between the Lie symmetry of a SIRS stochastic differential equation and its integrability is established in [26]. The derived equation is used to find the admitted Lie point symmetries of the stochastic SIRS model. We provide the message passing system of equations for this model and prove, for the first time, that it has a unique feasible solution. A generalized upper bound is obtained for the expected epidemic size at a fixed time $t>0$, where the graph (network) is symmetric with a homogeneous special case in [27]. In real life, external fluctuation is a key factor in the transmission of epidemic diseases. Therefore, it is important to include the randomness in the model, and hence stochastic SIRS models are considered to be more practical and realistic to study and simulate the evaluation of any epidemic disease. A great deal of work has been done on the application of stochastic SIRS models in epidemiology [15,19,25,28,29,30,31,32,33,34,35,36,37,38,39].

The aim of this work is to improve numerous studies existing in the field of stochastic SIRS models and provide an alternative technique to treat and analyze more complex biological and epidemiological models. To this end, we use a stable, robust and efficient numerical technique based on the Legendre collocation method with theoretical justifications.

The classical deterministic SIR epidemic model is given by:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{dS\left(t\right)}{dt}=\mathrm{\Lambda }-\frac{\beta I\left(t\right)S\left(t\right)}{1+a{I}^2\left(t\right)}+\nu R\left(t\right)-\kappa S\left(t\right),\hfill \\ \frac{dI\left(t\right)}{dt}=\frac{\beta I\left(t\right)S\left(t\right)}{1+a{I}^2\left(t\right)}-I\left(t\right)(\tau +\psi +\kappa ),\hfill \\ \frac{dR\left(t\right)}{dt}=(\tau +\psi )I\left(t\right)-(\kappa +\nu )R\left(t\right),\hfill \end{array}\right.\end{array}$$

(1)

In the proposed model, we observe that the whole population $N\left(t\right)$ is divided into a further three categories: the first category is susceptible S(t), the second one is infectious I(t) and the last one is recovered R(t). All the parameters are assumed to be positive. Constant recruitment of new susceptible individuals is denoted by $\mathrm{\Lambda }$, where $\beta$ is the transmission rate; $\kappa$ shows the per-capita natural mortality rate; disease causing mortality is denoted by $\tau$; the rate of recovery of the infected individuals per capita is denoted by $\psi$; hosts initially recovered are immune. The loss of immunity is denoted by $\nu >0$, by the symmetry of Brownian motion [40]. Using the invariance principle of LaSalle, the basic reproduction number $R_0$ are influenced by threshold dynamics. If $R_0<1$, then Equation (1) would indicate disease-free equilibrium $E_0=(\frac{\mathrm{\Lambda }}{\kappa },0,0)$ and will be stable asymptotically, and otherwise, endemic equilibrium of the proposed system will be stable when reproduction number $R_0>1$. The stochastic perturbations are assumed to be white noise. The system given in Equation (1) can be extrapolated to the following stochastic SIRS model:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}dS\left(t\right)=\left(\mathrm{\Lambda }-\frac{\beta I\left(t\right)S\left(t\right)}{1+a{I}^2\left(t\right)}+\nu R\left(t\right)-\kappa S\left(t\right)\right)dt+\frac{S\left(t\right)\sigma }{1+a{I}^2\left(t\right)}dW\left(t\right),\hfill \\ dI\left(t\right)=\left(\frac{\beta S\left(t\right)I\left(t\right)}{1+a{I}^2\left(t\right)}-I\left(t\right)(\tau +\kappa +\psi )\right)dt+\frac{\sigma I\left(t\right)}{1+a{I}^2\left(t\right)}dW\left(t\right),\hfill \\ dR\left(t\right)=\left(I\left(t\right)(\tau +\psi )-(\nu +\kappa )R\left(t\right)\right)dt+\frac{\sigma R\left(t\right)}{1+a{I}^2\left(t\right)}dW\left(t\right),\hfill \end{array}\right.\end{array}$$

(2)

where the initial conditions are

$$\begin{array}{c}\hfill S\left(0\right)={\tau }_1,I\left(0\right)={\tau }_2,R\left(0\right)={\tau }_3,\end{array}$$

$W\left(t\right)$ are independent standard Brownian motion, where the intensities of $W\left(t\right)$ are denoted by ${\sigma }^2>0$. This study includes the solution of the stochastic SIRS model in Equation (2) using the Legendre spectral collocation method (LSCM). Many authors have been used Legendre polynomials and their operational matrix for the approximation of integral and differential equations [12,14].

Sections of this research article are organized as follows: Section 2 includes a concise description of Legendre polynomials and the description of the scheme to solve the proposed model. Section 3 consists of stability analysis of the SIRS model. Numerical results and discussion are given in Section 4, while conclusions are drawn in Section 5.

2. Description of the Method

Before we describe the proposed scheme in detail for the stochastic SIRS model given in (2), we first briefly discuss Legendre polynomials.

Legendre Polynomials

Consider the Nth order Legendre polynomial, which is denoted by $P_N\left(t\right)$. We use Legendre interpolating polynomials to approximate any function $y\left(t\right)$ in $[-1,1]$ interval, which are defined as

$$\begin{array}{c}\hfill y\left(t\right)=\sum _{n=0}^N{P}_n\left(t\right)x_n\end{array}$$

(3)

where $t_n,n=0,1,\dots ,N$ are interpolating points and satisfy

$$\begin{array}{c}\hfill -1=t_0<\dots <t_N=1,\end{array}$$

and $y_n=y\left(t_n\right)$ denotes unknown Legendre coefficients, and Nth Legendre polynomials are denoted by $P_N\left(t\right)$ and defined as

$$\begin{array}{c}\hfill P_n\left(t\right)=\frac{1}{2^n}\sum _{m=0}^{\left[\frac{n}{2}\right]}{(-1)}^m{(}_n^{2n-2m}){(}_m^n)t^{n-2m},-1\le t\le 1,(n=0,1,\dots ,N)\end{array}$$

(4)

$$\left[\frac{n}{2}\right]=\left\{\begin{array}{cc}\frac{n}{2},\hfill & neven,\hfill \\ \frac{n-1}{2},\hfill & nodd.\hfill \end{array}\right.$$

To solve the proposed model given in Equation (2), we use Legendre Gauss quadrature along with a weight function. We consider, for the LSCM, the Legendre Gauss Lobatto points ${\left\{t_j\right\}}_{j=0}^N$.

Our aim in this study is to develop an efficient approximate solution to Equation (2). Integrating Equation (2) from $[0,t]$, we obtain

$$\begin{array}{cc}\hfill S\left(t\right)& =S\left(0\right)+{\int }_0^t\left(\mathrm{\Lambda }-\frac{\beta S\left(s\right)I\left(s\right)}{1+a{I}^2\left(s\right)}-\kappa S\left(s\right)+\nu R\left(s\right)\right)ds+{\int }_0^t\frac{\sigma S\left(s\right)}{1+a{I}^2\left(s\right)}dW\left(s\right),\hfill \\ \hfill I\left(t\right)& =I\left(0\right)+{\int }_0^t\left(\frac{\beta S\left(s\right)I\left(s\right)}{1+a{I}^2\left(s\right)}-(\kappa +\tau +\psi )I\left(s\right)\right)ds+{\int }_0^t\frac{\sigma I\left(s\right)}{1+a{I}^2\left(s\right)}dW\left(s\right),\hfill \\ \hfill R\left(t\right)& =R\left(0\right)+{\int }_0^t\left((\psi +\tau )I\left(s\right)-(\kappa +\nu )R\left(s\right)\right)ds+{\int }_0^t\frac{\sigma R\left(s\right)}{1+a{I}^2\left(s\right)}dW\left(s\right),\hfill \end{array}$$

(5)

where $S\left(0\right),I\left(0\right)$ and $R\left(0\right)$ are initial conditions for the functions $S\left(t\right),I\left(t\right)$ and $R\left(t\right)$, respectively. Taking linear transformation $s=\frac{t}{2}(1+\theta )$ to analyze the LSCM over standard interval $[-1,1]$ on Equation (5), we obtain

$$\begin{array}{cc}\hfill S\left(t\right)& =S\left(0\right)+\frac{t}{2}{\int }_{-1}^1\left(\mathrm{\Lambda }-\frac{\beta S\left(\frac{t}{2}(\theta +1)\right)I\left(\frac{t}{2}(1+\theta )\right)}{1+a{I}^2\left(\frac{t}{2}(\theta +1)\right)}-\kappa S\left(\frac{t}{2}(1+\theta )\right)+\nu R\left(\frac{t}{2}(\theta +1)\right)\right)d\theta \hfill \\ \hfill & +\frac{t}{2}{\int }_{-1}^1\frac{\sigma S\left(\frac{t}{2}(1+\theta )\right)}{1+a{I}^2\left(\frac{t}{2}(1+\theta )\right)}dW\left(\theta \right),\hfill \\ \hfill I\left(t\right)& =I\left(0\right)+\frac{t}{2}{\int }_{-1}^1\left(\frac{\beta S\left(\frac{t}{2}(\theta +1)\right)I\left(\frac{t}{2}(\theta +1)\right)}{1+a{I}^2\left(\frac{t}{2}(\theta +1)\right)}-(\kappa +\tau +\psi )I\left(\frac{t}{2}(\theta +1)\right)\right)d\theta \hfill \\ \hfill & +\frac{t}{2}{\int }_{-1}^1\frac{\sigma I\left(\frac{t}{2}(\theta +1)\right)}{1+a{I}^2\left(\frac{t}{2}(\theta +1)\right)}dW\left(\theta \right),\hfill \\ \hfill R\left(t\right)=& R\left(0\right)+\frac{t}{2}{\int }_{-1}^1\left((\psi +\tau )I\left(\frac{t}{2}(\theta +1)\right)-(\kappa +\nu )R\left(\frac{t}{2}(\theta +1)\right)\right)d\theta \hfill \\ \hfill & +\frac{t}{2}{\int }_{-1}^1\frac{\sigma R\left(\frac{t}{2}(\theta +1)\right)}{1+a{I}^2\left(\frac{t}{2}(\theta +1)\right)}dW\left(\theta \right),\hfill \end{array}$$

(6)

The spectral equations (semi-discretized) of Equation (6) are given by

$$\begin{array}{cc}\hfill S\left(t\right)& =S\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left(\mathrm{\Lambda }-\frac{\beta S\left(\frac{t({\theta }_k+1)}{2}\right)I\left(\frac{t(1+{\theta }_k)}{2}\right)}{1+a{I}^2\left(\frac{t(1+{\theta }_k)}{2}\right)}-\kappa S\left(\frac{t({\theta }_k+1)}{2}\right)+\nu R\left(\frac{t({\theta }_k+1)}{2}\right)\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N\left(\frac{\sigma S\left(\frac{t({\theta }_k+1)}{2}\right)}{1+a{I}^2\left(\frac{t({\theta }_k+1)}{2}\right)}\right){\omega }_k^{*},\hfill \\ \hfill I\left(t\right)& =I\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left(\frac{\beta S\left(\frac{t(1+{\theta }_k)}{2}\right)I\left(\frac{t(1+{\theta }_k)}{2}\right)}{1+a{I}^2\left(\left(\frac{t(1+{\theta }_k)}{2}\right)\right)}-(\tau +\kappa +\psi )I\left(\frac{t(1+{\theta }_k)}{2}\right)\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N\left(\frac{\sigma I\left(\frac{t(1+{\theta }_k}{2})\right)}{1+a{I}^2\left(\frac{t(1+{\theta }_k)}{2}\right)}\right){\omega }_k^{*},\hfill \\ \hfill R\left(t\right)& =R\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left((\tau +\psi )I\left(\frac{t(1+{\theta }_k)}{2}\right)-(\nu +\kappa )R\left(\frac{t(1+{\theta }_k)}{2}\right)\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N\left(\frac{\sigma R\left(\frac{t(1+{\theta }_k)}{2}\right)}{1+a{I}^2\left(\frac{t(1+{\theta }_k)}{2}\right)}\right){\omega }_k^{*},\hfill \end{array}$$

(7)

where Legendre Gauss quadrature with weights is

$$\begin{array}{c}\hfill {\omega }_k=\frac{2}{{\left[L_{N+1}^{{}^{\prime }}\left(s_k\right)\right]}^2(1-s_k^2)},0\le k\le N.\end{array}$$

Similarly, ${\omega }_k^{*}=\sqrt{{\omega }_k}\times randn(1,N)$ is the stochastic weight function.

To find the approximate solution for each class, we use the Legendre polynomials given in Equation (3):

$$\begin{array}{c}\hfill S\left(t\right)=\sum _{n=0}^N{S}_n{P}_n\left(t\right),I\left(t\right)=\sum _{n=0}^N{I}_n{P}_n\left(t\right),R\left(t\right)=\sum _{n=0}^N{R}_n{P}_n\left(t\right)\end{array}$$

(8)

In the above equation, Legendre coefficients of each class are denoted by $S_n=S\left(t_n\right),I_n=I\left(t_n\right)$ and $R_n=R\left(t_n\right)$, respectively. Incorporating Equation (8) into Equation (7), we have

$$\begin{array}{cc}\hfill \sum _{n=0}^N{S}_n{P}_n\left(t\right)& =\sum _{n=0}^N{S}_n{P}_n\left(0\right)+\frac{t}{2}\sum _{k=0}^N(\mathrm{\Lambda }-\frac{\beta {\sum }_{n=0}^N{S}_n{P}_n\left({\eta }_k\right){\sum }_{n=0}^N{S}_n{P}_n\left({\eta }_k\right)}{1+a{\left({\sum }_{n=0}^N{I}_n{P}_n\left({\eta }_k\right)\right)}^2}\hfill \\ \hfill & -\kappa \sum _{n=0}^N{S}_n{P}_n\left({\eta }_k\right)+\nu \sum _{n=0}^N{S}_n{P}_n\left({\eta }_k\right)){\omega }_k+\frac{t}{2}\sum _{k=0}^N\left(\frac{\sigma {\sum }_{n=0}^N{S}_n{P}_n\left({\eta }_k\right)}{1+a{\left({\sum }_{n=0}^N{I}_n{P}_n\left({\eta }_k\right)\right)}^2}\right){\omega }_k^S,\hfill \\ \hfill \sum _{n=0}^N{I}_n{P}_n\left(t\right)& =\sum _{n=0}^N{I}_n{P}_n\left(0\right)+\frac{t}{2}\sum _{k=0}^N(\frac{\beta {\sum }_{n=0}^N{S}_n{P}_n\left({\eta }_k\right){\sum }_{n=0}^N{S}_n{P}_n\left({\eta }_k\right)}{1+a{\left({\sum }_{n=0}^N{I}_n{P}_n\left({\eta }_k\right)\right)}^2}\hfill \\ \hfill & -(\kappa +\tau +\psi )\sum _{n=0}^N{I}_n{P}_n\left({\eta }_k\right)){\omega }_k+\frac{t}{2}\sum _{k=0}^N\left(\frac{\sigma {\sum }_{n=0}^N{I}_n{P}_n\left({\eta }_k\right)}{1+a{\left({\sum }_{n=0}^N{I}_n{P}_n\left({\eta }_k\right)\right)}^2}\right){\omega }_k^I,\hfill \\ \hfill \sum _{n=0}^N{R}_n{P}_n\left(t\right)& =\sum _{n=0}^N{R}_n{P}_n\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left((\psi +\tau )\sum _{n=0}^N{I}_n{P}_n\left({\eta }_k\right)-(\kappa +\tau )\sum _{n=0}^N{R}_n{P}_n\left({\eta }_k\right)\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N\left(\frac{\sigma {\sum }_{n=0}^N{R}_n{P}_n\left({\eta }_k\right)}{1+a{\left({\sum }_{n=0}^N{I}_n{P}_n\left({\eta }_k\right)\right)}^2}\right){\omega }_k^R,\hfill \end{array}$$

(9)

taking ${\eta }_k=\left(\frac{t}{2}({\theta }_k+1)\right)$ for simplicity. Thus, there are $3N+3$ unknowns in the system given in Equation (9) with $3N$ non-linear algebraic equations. After incorporating the initial conditions, we obtain

$$\begin{array}{c}\hfill \sum _{n=0}^N{S}_n{P}_n\left(0\right)={\tau }_1,\sum _{n=0}^N{I}_n{P}_n\left(0\right)={\tau }_2\sum _{n=0}^N{R}_n{P}_n\left(0\right)={\tau }_3.\end{array}$$

(10)

These two equations—Equations (9) and (10)—consist of $3N+3$ non-linear systems of algebraic equations having $3N+3$ variables (unknowns). We obtain a numerical solution to the system given in Equation (2) by incorporating the values of these unknowns into Equation (8).

3. Stability Analysis

In this section, we take into account both the stochastic SIRS model of Equation (2) and the deterministic SIRS model of Equation (1) considering saturated incident rate $\frac{\beta S\left(t\right)I\left(t\right)}{1+a{I}^2\left(t\right)}$ and a time describing inherent period. The basic reproduction number for Equation (1) is defined by $R_0=\frac{\mathrm{\Lambda }\beta }{\kappa (\kappa +\tau +\psi )}.$ Below, we have given some important results about the basic reproduction number.

Theorem 1.

The solution to the system given in Equation (1) gives a disease-free equilibrium $E_0=(\frac{\mathrm{\Lambda }}{\kappa },0,0)$ if reproduction number $R_0$ is less than 1 and the solution to the system given in Equation (1) gives unique endemic equilibrium $E^{*}=(S^{*},I^{*},R^{*})$ if reproduction number $R_0$ is greater than 1. Here,

$$\begin{array}{cc}\hfill & S^{*}=\frac{(\tau +\kappa +\psi )\left[1+(1+a)a^2{(\tau +\psi +\kappa )}^2\right]}{\beta \left[1+{(\tau +\psi +\kappa )}^2{a}^2\right]},\hfill \\ \hfill & I^{*}=\frac{(\tau +\psi +\kappa )a}{\sqrt{\left(1+{(\kappa +\tau +\psi )}^2{a}^2\right)}},\hfill \\ \hfill & R^{*}=\frac{(\tau +\psi +\kappa )(\psi +\tau )a}{(\kappa +\nu )\sqrt{\left[1+a^2{(\tau +\psi +\kappa )}^2\right]}}.\hfill \end{array}$$

Proof. 

Let the endemic equilibrium of the system given in Equation (1) be $E^{*}$, and the given classes satisfy the stationary system

$$\begin{array}{cc}\hfill & \mathrm{\Lambda }-\frac{\beta S^{*}\left(t\right)I^{*}\left(t\right)}{1+a{I}^{2^{*}}\left(t\right)}-\kappa S^{*}\left(t\right)+\nu R^{*}\left(t\right)=0,\hfill \\ \hfill & \frac{\beta S^{*}\left(t\right)I^{*}\left(t\right)}{1+a{I}^{2^{*}}\left(t\right)}-(\kappa +\tau +\psi )I^{*}\left(t\right)=0,\hfill \\ \hfill & (\tau +\psi )I^{*}\left(t\right)-(\nu +\kappa )R^{*}\left(t\right)=0.\hfill \end{array}$$

(11)

We will consider two cases. Case 1: $I^{*}\left(t\right)=0$; case 2: $I^{*}\left(t\right)>0$.

• Case 1: If $I^{*}\left(t\right)=0$, using the third equation of the system given in Equation (11) gives $R^{*}\left(t\right)=0$, and it follows from the first equation of the system in Equation (11) that $S^{*}\left(t\right)=\frac{\mathrm{\Lambda }}{\kappa }$. Thus, the solution of Equation (11) reduces to $E_0=(\frac{\mathrm{\Lambda }}{\kappa },0,0)$.
• Case 2: Here, $I^{*}\left(t\right)>0$. In this case, the calculation has been done using Maple software version 13. We successfully obtain the following stable endemic equilibrium $S^{*},I^{*},R^{*}$:

  $$\begin{array}{c}\hfill R^{*}=\frac{a(\kappa +\tau +\psi )(\tau +\psi )}{(\kappa +\nu )\sqrt{}(1+a^2{(\kappa +\psi +\tau )}^2)}.\end{array}$$

  $$\begin{array}{c}\hfill S^{*}=\frac{(\tau +\kappa +\psi )\left(1+(1+a)a^2{(\tau +\kappa +\psi )}^2\right)}{\beta \left(1+{(\tau +\kappa +\psi )}^2{a}^2\right)},\end{array}$$

  and

  $$\begin{array}{c}\hfill I^{*}=\frac{(\kappa +\tau +\psi )a}{\sqrt{}(1+a^2{(\kappa +\tau +\psi )}^2)},.\end{array}$$

Thus, if $I^{*}>0$, we must have $\frac{(\tau +\kappa +\psi )a}{\sqrt{}{(1+(\tau +\kappa +\psi )}^2{a}^2)}>0$, which implies that $(\tau +\kappa +\psi )>0$. This means that we must ensure that $R_0>1$. This completes the proof of the theorem. □

Lemma 1.

Initial values $(S\left(0\right),I\left(0\right),R\left(0\right))\in R_{+}^3$, for the system given in Equation (2), have almost surely a unique positive solution $(S\left(t\right),I\left(t\right),R\left(t\right))\in R_{+}^3$ for $t\ge 0$.

Lemma 2.

Total region $\mathrm{\Omega }=\{(S\left(t\right),I\left(t\right),R\left(t\right))\in R_{+}^3:$ are positive; $S\left(t\right)+I\left(t\right)+R\left(t\right)\le \frac{\mathrm{\Lambda }}{\kappa }\}$ have a positively invariant solution for the proposed system given in Equation (2).

Proof. 

As the total population are $S\left(t\right)+I\left(t\right)+R\left(t\right)=N\left(t\right)$, differentiating and using Equation (2), we obtain the following:

$$\begin{array}{cc}\hfill & \frac{dN\left(t\right)}{dt}=\mathrm{\Lambda }+\kappa N\left(t\right)-\tau I\left(t\right)\hfill \\ \hfill & \frac{dN\left(t\right)}{dt}\le \mathrm{\Lambda }+\kappa N\left(t\right)\hfill \end{array}$$

(12)

Integrating and using the initial values gives

$$\begin{array}{cc}\hfill N\left(t\right)& \le \frac{\mathrm{\Lambda }}{\kappa }+N\left(0\right)e^{-\kappa t}\hfill \\ \hfill & \le \frac{\mathrm{\Lambda }}{\kappa }.\hfill \end{array}$$

 □

Lemma 3.

Using the initial conditions, the solution of $\left(S\right(t),I(t),R(t\left)\right)$ of the system in Equation (2) satisfies the following properties:

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^t\frac{\sigma S\left(s\right)}{1+a{I}^2\left(s\right)}dW\left(s\right)=0,\underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^t\frac{\sigma I\left(s\right)}{1+a{I}^2\left(s\right)}dW\left(s\right)=0,\hfill \\ \hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^t\frac{\sigma R\left(s\right)}{1+a{I}^2\left(s\right)}dW\left(s\right)=0,\hfill \end{array}$$

(13)

Definition 1.

The infected class “I(t)” is said to be extinctive if ${lim}_{t\to \infty }I\left(t\right)=0$.

Let us propose

$$\begin{array}{c}\hfill {\overline{R}}_0=R_0-\frac{{\sigma }^2{\mathrm{\Lambda }}^2{e}^{-\kappa \tau }}{2{\kappa }^2(\kappa +{\rho }_1+{\rho }_2+\tau )}\end{array}$$

and we have the theorem as given below.

Theorem 2.

When $max\left\{\frac{\beta \kappa e^{-\kappa \tau }}{\mathrm{\Lambda }},\frac{{\beta }^2{e}^{-2\kappa \tau }}{2(\kappa +{\rho }_1+{\rho }_2+\tau )}\right\}<{\sigma }^2$ or ${\sigma }^2<\frac{\beta \kappa e^{-\kappa \tau }}{\mathrm{\Lambda }}$ along ${\overline{R}}_0<1$, then infected class $I\left(t\right)$ of the above stochastic system given in Equation (2) exponentially approaches zero. Here,

$$\begin{array}{c}\hfill {\overline{R}}_0=R_0-\frac{{\sigma }^2{\mathrm{\Lambda }}^2{e}^{-\kappa \tau }}{2{\kappa }^2(\kappa +{\rho }_1+{\rho }_2+\tau )}\end{array}$$

Proof. 

Consider the solution of the above system given in Equation (2) entertaining the initial conditions $(S\left(0\right),I\left(0\right),R\left(0\right))\in R_{+}^3$. Then, by using the It $\widehat{o}$ formula, we find

$$\begin{array}{cc}\hfill dln\left(I\right(t\left)\right)=& \left(\frac{\beta e^{-\kappa \tau }}{I\left(t\right)\tau +1}S\left(t\right)-({\rho }_1+\kappa +{\rho }_2+\tau )-\frac{{\sigma }^2{e}^{-\kappa \tau }}{2{[I\left(t\right)\tau +1]}^2}{S}^2\left(t\right)\right)dt\hfill \\ \hfill & +\sigma S\left(t\right)dW\left(t\right).\hfill \end{array}$$

(14)

Now, taking the integration of the system in Equation (3) from 0 to t shows

$$\begin{array}{cc}\hfill lnI\left(t\right)=& lnI\left(0\right)+{\int }_0^t\left(\frac{\beta e^{-\kappa \tau }}{1+\tau I\left(t\right)}S\left(t\right)-(\kappa +{\rho }_1+{\rho }_2+\tau )-\frac{{\sigma }^2{e}^{-\kappa \tau }}{2{(1+\tau I\left(t\right))}^2}{S}^2\left(t\right)\right)dt\hfill \\ \hfill & +{\int }_0^t\sigma S\left(t\right)dW\left(t\right).\hfill \end{array}$$

(15)

Here, we discuss two cases:

• Case 1: when ${\sigma }^2>\frac{\beta \kappa e^{-\kappa \tau }}{\mathrm{\Lambda }}$, Equation (15) becomes

  $$\begin{array}{c}\hfill lnI\left(t\right)\le lnI\left(0\right)+\left(\frac{{\beta }^2{e}^{-\kappa \tau }}{{\sigma }^2}-(\kappa +{\rho }_1+{\rho }_2+\tau )\right)t+{\int }_0^t\sigma S\left(t\right)dW\left(t\right).\end{array}$$

  (16)

  Dividing both sides of Equation (16) by $t>0$, we obtain

  $$\begin{array}{c}\hfill \frac{lnI\left(t\right)}{t}\le \frac{1}{t}lnI\left(0\right)-\left(\kappa +{\rho }_1+{\rho }_2+\tau -\frac{{\beta }^2{e}^{-\kappa \tau }}{{\sigma }^2}\right)+\frac{1}{t}{\int }_0^t\sigma S\left(t\right)dW\left(t\right).\end{array}$$

  (17)

  Taking ${lim}_{t\to \infty }$ on Equation (17) and using Lemma 3, we have

  $$\begin{array}{c}\hfill \underset{t\to \infty }{lim}\frac{lnI\left(t\right)}{t}\le -\left(\kappa +{\rho }_1+{\rho }_2+\tau -\frac{{\beta }^2{e}^{-\kappa \tau }}{{\sigma }^2}\right)<0.\end{array}$$

  when ${\sigma }^2>\frac{{\beta }^2{e}^{-2\kappa \tau }}{2(\kappa +{\rho }_1+{\rho }_2+\tau )}$, which shows that ${lim}_{t\to \infty }I\left(t\right)=0.$
• Case 2: When ${\sigma }^2<\frac{\beta \kappa }{\mathrm{\Lambda }}$, then from Equation (15), we have

  $$\begin{array}{c}\hfill lnI\left(t\right)\le lnI\left(0\right)+\left(\frac{\beta \mathrm{\Lambda }{e}^{-\kappa \tau }}{\kappa }-\frac{{\sigma }^2{\mathrm{\Lambda }}^2{e}^{-\kappa \tau }}{2{\kappa }^2}-(\kappa +{\rho }_1+{\rho }_2+\tau )\right)t+{\int }_0^t\sigma S\left(t\right)dW\left(t\right).\end{array}$$

  (18)

  Dividing both sides of Equation (18) by $t>0$, we obtain

  $$\begin{array}{cc}\hfill \frac{1}{t}lnI\left(t\right)& \le \frac{1}{t}lnI\left(0\right)+(\kappa +{\rho }_1+{\rho }_2+\tau )\left(\frac{\beta \mathrm{\Lambda }{e}^{-\kappa \tau }}{\kappa (\kappa +{\rho }_1+{\rho }_2+\tau )}-\frac{{\sigma }^2{\mathrm{\Lambda }}^2{e}^{-\kappa \tau }}{2{\kappa }^2(\kappa +{\rho }_1+{\rho }_2+\tau )}-1\right)\hfill \\ \hfill & +\frac{1}{t}{\int }_0^t\sigma S\left(t\right)dW\left(t\right).\hfill \end{array}$$

  (19)

  Again, taking ${lim}_{t\to \infty }$ on Equation (17) and using Lemma 3, we have

  $$\begin{array}{c}\hfill \underset{t\to \infty }{lim}\frac{lnI\left(t\right)}{t}\le \left({\rho }_1+\kappa +\tau +{\rho }_2\right)({\overline{R}}_0-1).\end{array}$$

  When ${\overline{R}}_0<1$,

  $$\begin{array}{c}\hfill \underset{t\to \infty }{lim}\frac{lnI\left(t\right)}{t}<0.\end{array}$$

  which shows that ${lim}_{t\to \infty }I\left(t\right)=0.$ This completes the proof. □

Theorem 3.

If ${\overline{R}}_0>1$, then $I\left(t\right)$ are present in Equation (2).

Proof is given in [31].

Remark 1.

After the above discussion, we can make the following remarks. Making use of Theorem 1, Equation (1) is disease-free if $R_0<1$ and it has a unique endemic equilibrium if $R_0>1$. By Theorem 2, if ${\sigma }^2>max\left\{\frac{\beta \kappa e^{-\kappa \tau }}{\mathrm{\Lambda }},\frac{{\beta }^2{e}^{-2\kappa \tau }}{2(\kappa +{\rho }_1+{\rho }_2+\tau )}\right\}$, the infectious disease in the system of Equation (2) is almost surely extinct, whereas the system in (2) has a disease-free equilibrium for small intensities and ${\overline{R}}_0<1$, which shows that ${\overline{R}}_0$ is the threshold dynamics along with the extinction of of infectious diseases.

4. Numerical Discussion

Numerical simulations are performed using our proposed numerical scheme for both the deterministic and stochastic SIR models to show the different behaviors of reproduction number $R_0$ using different parameter values involved in both mathematical models shown in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8. In all simulations, we use $S_0=0.8,I_0=0.15,R_0=0.05$. In Figure 1, we use the parameter values $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =1,\tau =0.2,\psi =0.5$, $\sigma =0.8$. For each class, we compare both the models given in Equations (1) and (2), respectively. Using the above parameter values, simple calculation shows that both reproduction numbers $R_0=0.58823<1$ and ${\overline{R}}_0=0.4983<1$. Using Theorems 1 and 2, we clearly observe in Figure 1 that the given systems have a stable disease-free equilibrium $\left(\mathrm{\Lambda }/\kappa ,0,0\right)$, where the susceptible class is $\mathrm{\Lambda }/\kappa =1$. The parameter values $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =2,\tau =0.2,\psi =0.5$, $\sigma =0.8$ were used to obtain Figure 2. For each class, we compare both the models given in Equations (1) and (2), where calculation gives $R_0=1.17647>1$ and ${\overline{R}}_0=1.12148>1$. Applying Theorems 1 and 2, both systems in Equations (1) and (2) have stable endemic equilibrium $E^{*}$. Moreover, Figure 2 shows that all the classes are greater than zero. Using $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =1,\tau =0.2,\psi =0.5,\sigma =0$ parameter values in Equation (1), Figure 3 is obtained with reproduction number ${\overline{R}}_0=0.58823<1$. Using Theorem 1, we notice from Figure 3 that the given systems satisfy the disease-free equilibrium. Moreover, the susceptible class is given by $\mathrm{\Lambda }/\kappa =1,$ and all the remaining classes approach zero. For Figure 4, we use $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =2,\tau =0.2,\psi =0.5,\sigma =0$ parameter values and we draw the numerical solution of the deterministic system given in Equation (1), with reproduction number ${\overline{R}}_0=1.17647>1$. Again, using Theorem 1, Figure 4 satisfies the stable endemic equilibrium. Moreover, all the classes become greater than zero. Figure 5 is obtained using the parameter values $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =1,\tau =0.2,\psi =0.5,\sigma =0.8$ in Equation (1), with reproduction number ${\overline{R}}_0=0.4983<1$. Now, using Theorem 2, we observe from Figure 5 that the given systems satisfy the disease-free equilibrium. Moreover, the susceptible class is given by $\mathrm{\Lambda }/\kappa =1,$ and all the remaining classes approach zero. Using $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =2,\tau =0.2,\psi =0.5,\sigma =0.8$ parameter values, the result is as shown in Figure 6. The numerical solution of the stochastic system given in Equation (1) is obtained with ${\overline{R}}_0=1.12148>1$. Again, by Theorem 2, Figure 6 shows that the systems satisfy the stable endemic equilibrium. Moreover, all the classes become greater than zero. Using the parameter values $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =1,\tau =0.2,\psi =0.5,\sigma =0.8$ and $R_0=0.58823<1$ and ${\overline{R}}_0=0.4983<1$, Figure 7 is obtained. By Theorem 2, Figure 7 satisfies the disease-free equilibrium. Moreover, the susceptible class is given by $\mathrm{\Lambda }/\kappa =1,$ and all the remaining classes approach zero. For the parameter values $\mathrm{\Lambda }=1,\kappa =1,a=0.5,\nu =0.4,\beta =2,\tau =0.2,\psi =0.5,\sigma =0.8$, the numerical solution for both systems given in Equations (1) and (2) with reproduction numbers $R_0=1.17647>1$ and ${\overline{R}}_0=1.12148>1$ is as shown in Figure 8. Now, by using Theorem 2, we evidently notice from Figure 8 that the given systems satisfy the stable endemic equilibrium. Moreover, all the classes become greater than zero.

5. Conclusions

The Legendre spectral method and its orthogonal properties are used to obtain the numerical solution to the proposed SIRS model along with stochastic/random phenomena. To evaluate the integral term more efficiently, Legendre Gauss quadrature along with the weight function are used. Different behavior of the basic reproduction number was observed. It is shown that infected classes approach zero when the reproduction number $R_0<1$ and the disease-free equilibrium $E_0$ is found. It is observed that the system gives unique stable equilibrium $E^{*}$ when the reproduction number $R_0>1$. Very good agreement between the solutions of both the stochastic and deterministic systems is found. The assumption of the strong symmetry of the initial conditions and parameter values plays a key role in the analysis of the SIRS stochastic model. Our simulations show that the extinction conditions of the disease are weaker in the case of stochastic white noise as compared to its deterministic counterpart. Our findings also show that the rate of saturated incidence has a crucial role in controlling the spread of the disease within the population, which confirms that stochastic epidemic models based on virus dynamics are more realistic as compared to their deterministic counterparts. Our theoretical results provide a solid foundation for the study of similar diseases and has a wide range of applications in the field of biomedical sciences. Our future goal is to apply spectral methods to stochastic SIRS models of fractional order.