A Comparison of Deterministic and Stochastic Susceptible-Infected-Susceptible (SIS) and Susceptible-Infected-Recovered (SIR) Models

Abdelmalik Moujahid, Fernando Vadillo^*

Department of Mathematics, University of the Basque Country, Biscay, Spain.

DOI: 10.4236/ojmsi.2021.93016                

Abstract

In this paper we build and analyze two stochastic epidemic models with death. The model assumes that only susceptible individuals (S) can get infected (I) and may die from this disease or a recovered individual becomes susceptible again (SIS model) or completely immune (SIR Model) for the remainder of the study period. Moreover, it is assumed there are no births, deaths, immigration or emigration during the study period; the community is said to be closed. In these infection disease models, there are two central questions: first it is the disease extinction or not and the second studies the time elapsed for such extinction, this paper will deal with this second question because the first answer corresponds to the basic reproduction number defined in the bibliography. More concretely, we study the mean-extinction of the diseases and the technique used here first builds the backward Kolmogorov differential equation and then solves it numerically using finite element method with FreeFem++. Our contribution and novelty are the following: however the reproduction number effectively concludes the extinction or not of the disease, it does not help to know its extinction times because example with the same reproduction numbers has very different time. Moreover, the SIS model is slower, a result that is not surprising, but this difference seems to increase in the stochastic models with respect to the deterministic ones, it is reasonable to assume some uncertainly.

Keywords

Persistence Time, Epidemic Dynamics, Stochastic Epidemic Models, Finite Element Method

A population is like a billiard ball: you get a lot of variability, but the variability is random, in all directions.

Stephen Jay Gould: The Pattern of Life’s History [1].

1. Introduction

The ordinary differential equations in epidemic models has been a well-known topic for some time, there exist classic book like ( [2], Chap. 10) or [3] and more recent monographs: [4] and [5]. All these deterministic models serve as a framework for formulating analogous stochastic models and as a source of comparison with the stochastic models (see for example [6] for a complete reference guide to all important contributions in stochastic epidemic models). In this way, there have been two ways to pass from the ordinary system to the stochastic system, one of them is simply adding new stochastic terms, for example in [7] [8] [9] [10]. The other one is explained in [11] [12] and it is used in this paper. This technique begins by assuming different probabilities of the changes and calculating means and covariance matrix to obtain a stochastic differential system. But later we’ll see that the stochastic parts are different and this will cause great differences in the asymptotic behavior of the solutions (see for example this results in [13] for a simple example about a population model).

The Stochastic Differential Equation (SDE) system for the dynamics of n variables has the form

$dX=\mu \left(t,X\right)+B\left(t,X\right)dW,$ (1)

where $X={\left(X_1\mathrm{,}\cdots \mathrm{,}{X}_n\right)}^{\text{T}}$ and $W={\left(W_1\mathrm{,}\cdots \mathrm{,}{W}_n\right)}^{\text{T}}$ are n independent Wiener processes. The vectorial function $\mu \left(t\mathrm{,}X\right)$ is called the drift, and $B\left(t,X\right)$ the diffusion matrix. Moreover, the components ${\mu }_j=\mathbb{E}\left(X_j\right),j=1,\cdots ,n$ and the diffusion term is the square root of the covariance matrix $D\left(t\mathrm{,}X\right)$ i.e. $B\cdot B=D$.

Obviously, a key question in the epidemic models is to understand the constrains that lead to extinction or not of the disease and when. In order to study this question, let define the random variable T that indicates the persistence time i.e. the time it takes for the size of either variables to reach zero.

$T\equiv \inf \left\{t\ge 0:X_j=0,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\text{for}\text{some}j=1,\cdots ,n\right\},$

obviously T depends on the initial value $X\left(0\right)$ although it is not explicitly indicated.

As discussed in ( [2], p. 150), the mean persistence-time $\tau \equiv \mathbb{E}\left(T\right)$ for (1) satisfies the stationary backward Kolmogorov equation

$\mathfrak{L}\left(\tau \right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\equiv \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\mu }_j\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\frac{\partial \tau }{x_j}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\frac{1}{2}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{D}_{jk}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\partial }^2\tau }{\partial x_j\partial x_k}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-1.$ (2)

and with boundary conditions

$\{\begin{array}{l}\tau \left(\cdots ,x_{j-1},0,x_{j+1},\cdots \right)=0,\\ \frac{\partial \tau }{\partial x_j}\left(\cdots ,x_{j-1},M_j,x_{j+1}\cdots \right)=0,\end{array}$ (3)

assuming that $x_j\le M_j,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}j=1,\cdots ,n$.

The Equation (2) is an Elliptic Partial Differential Equations of Second Order [14], really it is a advection diffusion equation, and as the name suggests, the mean persistence-time will depend on the operator $D$. These comments would explain the results in [13]. Moreover in [7] [8] [9] [10] the matrices $B$ are diagonals which implies that the variables are not correlated, an unreasonable hypothesis. As already commented in the abstract, we will solved numerically the backward Kolmogorov equation using finite element method with FreeFem++. These authors et. al. have used this technique en several paper: [13] [15] [16] [17] or [18] with very hopeful results to spread more complex problems.

This paper is organized as follows. In Section 2, we describe the SIS and SIR models with death and obtain the backward Kolmogorov equations. In Section 3 we explain the numerical methods with their results, with special attention to variational formulation of the problem (10) - (11). Finally in Section 4 we analyze the numerical results and draw the main conclusions.

Our numerical methods were implemented in Matlab© and FreeFem++ which is freely available and particularly efficient, see [19]. The experiments were carried out in an Intel(R) Core(TM)i7-8665U CPU @ 1.90 GHz, 16.0 GB of RAM. The codes for the numerical tests are available on request.

2. Two Models with Death

2.1. Stochastic SIS Model

In the SIS epidemic models the variable $X={\left(S\mathrm{,}I\right)}^{\text{T}}$ where $S=S\left(t\right)$ is the susceptible population size and $I=I\left(t\right)$ is the infected population size. In this model an individual starts off susceptible, at some stage catches the disease and after a short infectious period become susceptible again. Such a model is appropriate for a bacterial disease such as pneumococcus.

The SIS model without demography is considerably simpler than the one considered here (see for example [17] and its references). The SIS model with demography was proposed in [20] or [21], in particular we consider only death of infected individuals with rate $\beta$. The changes and their probabilities to the first order in $\Delta t$ are given in Table 1 with $x={\left(S\mathrm{,}I\right)}^{\text{T}}$.

Following ( [11], p. 148) we obtain the following deterministic model:

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\mu }_1\left(S,I\right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=-\frac{\alpha SI}{S+I}+\gamma I,\\ \frac{\text{d}I}{\text{d}t}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\mu }_2\left(S,I\right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\frac{\alpha SI}{S+I}-\left(\gamma +\beta \right)I,\end{array}$ (4)

This model has a fixed point at $I=0$ because ${\mu }_1\left(S,0\right)={\mu }_2\left(S,0\right)=0$ and the solution approaches the disease-free equilibrium ${\lim }_{t\to \infty }I\left(t\right)=0$ when a

Table 1. Possible change in $x={\left(S\mathrm{,}I\right)}^{\text{T}}$ and their probabilities.

new basic reproduction number satisfies $R=\alpha /\left(\gamma +\beta \right)\le 1$ (see [17] or [22] ).

The covariance matrix

$D\left(S\mathrm{,}I\right)=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(\begin{array}{cc}d\left(S\mathrm{,}I\right)& -d\left(S\mathrm{,}I\right)\\ -d\left(S\mathrm{,}I\right)& d\left(S\mathrm{,}I\right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\beta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}I\end{array}\right)\mathrm{,}$ (5)

with $d\left(S\mathrm{,}I\right)=\frac{\alpha SI}{S+I}+\gamma I$.

Then if $S\left(0\right)=x,I\left(0\right)=y$ its backward Kolmogoronv equation is

${\mu }_1\left(x,y\right)\frac{\partial \tau }{\partial x}+{\mu }_2\left(x,y\right)\frac{\partial \tau }{\partial y}+\frac{1}{2}d\left(x,y\right)\left(\frac{{\partial }^2\tau }{\partial x^2}-2\frac{{\partial }^2\tau }{\partial x\partial y}+\frac{{\partial }^2\tau }{\partial y^2}\right)+\frac{\beta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}y}{2}\frac{{\partial }^2\tau }{\partial y^2}=-1,$ (6)

with a singularity at the boundary $y=0$ because for this value the coefficient of the second derivative is cancelled, which is why boundary payer will arise. The boundary conditions are:

$\{\begin{array}{l}\tau \left(0,x\right)=\tau \left(y,0\right)=0,\\ \frac{\partial \tau }{\partial x}\left(M_x,y\right)=0,\\ \frac{\partial \tau }{\partial y}\left(x,M_y\right)=0,\end{array}$ (7)

provided that the number of x and y cannot exceed some values $M_x$ and $M_y$ respectively.

2.2. Stochastic SIR Model

In the SIR epidemic models the variable $X={\left(S\mathrm{,}I\mathrm{,}R\right)}^{\text{T}}$ where $S=S\left(t\right)$ is the susceptible population size, $I=I\left(t\right)$ is the infected population size and $R\left(t\right)$ the recovered individuals because are assumed that individuals are immune after recovering from the disease. For example, [7] [8] [9] [10] proposes a model that in view of the results of [13] does not have any practical value, in my opinion.

Now, the changes and their probabilities to the first order in $\Delta t$ are given in Table 2 with $x={\left(S\mathrm{,}I\mathrm{,}R\right)}^{\text{T}}$.

Quite similar the model in [23] plus a new parameter $\beta$. The deterministic model is

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-\frac{\alpha SI}{S+I},\\ \frac{\text{d}I}{\text{d}t}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\frac{\alpha SI}{S+I}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(\gamma +\beta \right)I,\\ \frac{\text{d}R}{\text{d}t}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\gamma I,\end{array}$ (8)

Table 2. Possible change in $x={\left(S\mathrm{,}I\mathrm{,}R\right)}^{\text{T}}$ and their probabilities.

In this differential system, the surface $I=0$ is fix and its jacobian has three eigenvalues: $\lambda =0$ doble and $\lambda =\alpha -\left(\gamma +\beta \right)$ and then when

$\gamma \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\beta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; 1\; <mtext>\; \hspace{0.17em}\; </mtext>}<\alpha <\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; 1\; <mtext>\; \hspace{0.17em}\; </mtext>}+\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\gamma \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\beta ,$

and the matrix covariance

$D\left(S\mathrm{,}I\mathrm{,}R\right)=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(\begin{array}{ccc}\frac{\alpha SI}{S+I}& -\frac{\alpha SI}{S+I}& 0\\ -\frac{\alpha SI}{S+I}& \frac{\alpha SI}{S+I}+\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(\gamma +\beta \right)I& -\gamma I\\ 0& -\gamma I& \gamma I\end{array}\right)\mathrm{.}$ (9)

Now assuming $S\left(0\right)=x,I\left(0\right)=y,R\left(0\right)=z$ the backward Kolmogoronv equation is

$\begin{array}{l}-d\left(x,y\right)\frac{\partial \tau }{\partial x}+\left(d\left(x,y\right)-\left(\gamma +\beta \right)y\right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\frac{\partial \tau }{\partial y}+\gamma y\frac{\partial \tau }{\partial z}\\ +\frac{1}{2}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}d\left(x,y\right)\left(\frac{{\partial }^2\tau }{\partial x^2}-2\frac{{\partial }^2\tau }{\partial x\partial y}+\frac{{\partial }^2\tau }{\partial y^2}\right)\\ +\frac{1}{2}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(\gamma +\beta \right)y\frac{{\partial }^2\tau }{\partial y^2}+\frac{1}{2}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\gamma y\frac{{\partial }^2\tau }{\partial z^2}-\gamma y\frac{{\partial }^2\tau }{\partial y\partial z}=-1\end{array}$ (10)

where $d\left(x\mathrm{,}y\right)=\frac{\alpha xy}{x+y}$, and the boundary conditions

$\{\begin{array}{l}\tau \left(0,y,z\right)=\tau \left(x,0,z\right)=\tau \left(x,y,0\right)=0,\text{on}{\Gamma }_1,\\ \frac{\partial \tau }{\partial x}\left(M_x,y,z\right)=0,\text{on}{\Gamma }_2,\\ \frac{\partial \tau }{\partial y}\left(x,M_y,z\right)=0,\text{on}{\Gamma }_3,\\ \frac{\partial \tau }{\partial z}\left(x,y,M_z\right)=0,\text{on}{\Gamma }_4,\end{array}$ (11)

provided that the number of x, y and z cannot exceed some values $M_x\mathrm{,}{M}_y$ and $M_z$ respectively.

3. Numerical Results

3.1. Deterministic Models

In the deterministic models, the value of the basic reproductive number ${\Re }_0=\frac{\alpha }{\beta +\gamma }$ determines persistence or extinction of the disease. If ${\Re }_0<1$, the disease is eliminated, whereas if ${\Re }_0>1$, the disease persists in the population (see for example [22] when the total population size remains constant). However, this parameter does not inform about the speed towards one situation or another.

Because these kinds of models are appropriate for a bacterial disease such as pneumococcus consider in [21], they study the pneumococcus in the population of Scottish children under two years, first we test this realistic example (Case 1) and later two example with similar value for ${\Re }_0$.

Let’s consider the following example.

Case 1: In [17] we studied the pneumococcus in the population of Scottish children with $\alpha =0.0214,\gamma =0.0211,\beta =0.00137\Rightarrow {\Re }_0=0.9986<1$.

Case 2: For $\alpha =0.1997,\gamma =0.1,\beta =0.1\Rightarrow {\Re }_0=0.9986<1$.

Case 3: For $\alpha =0.3994,\gamma =0.3,\beta =0.1\Rightarrow {\Re }_0=0.9986<1$.

We solved the two deterministic models (4) and (8) using Matlab. In Figure 1 we have plotted two paths for $0\le t\le 25000$ in blue is the SIS Model with $S\left(0\right)=1000,I\left(0\right)=10$ and in red the SIR Model with $S\left(0\right)=1000,I\left(0\right)=10,R\left(0\right)=1$, in both path $I\left(t\right)\to 0$ but in very different directions, a result that we would expect. Moreover, it would be more natural to ask Matlab to return the time $t^{\mathrm{<mo>\; *\; </mo>}}$ such that $0<I\left(t^{*}\right)=I^{*}\ll 1$, which can be done using the event location facility. The highlighted numerical results are summarized in Table 3. In this table we have noted the time $t_0$ where $I\left(t^{\mathrm{<mo>\; *\; </mo>}}\right)\approx 0.99$ in the first line, and $I\left(t^{\mathrm{<mo>\; *\; </mo>}}\right)\approx {10}^{-3}$ in the second.

These three examples with the same reproduction number already have slight differences in their extinction, the third case is by far the fastest as well as the second model SIR.

Figure 1. Case 1 for $0\le t\le 25000$, blue SIS Model for $S\left(0\right)=10000,I\left(0\right)=10$ and red SIR model with $R\left(0\right)=1$.

Table 3. Time $t^{\mathrm{<mo>\; *\; </mo>}}$ for $S\left(0\right)=1000,I\left(0\right)=10$.

3.2. Mean Extinction-Time in Each Model

These authors have resolved the boundary problems (6) - (7) using the Finite Element Method (FEM) and the software FreeFem++ which is freely available [19], this is a standard technique hence we will skip the details [17]. In Figure 2 we can see the numerical solution of (6) - (7) with $M_x=1000,M_y=100$. Note that the three graphs are very similar but the scale is quite different, the value of its peaks can be read in Table 4.

Figure 2. The numerical solution of (6) - (7) with $M_x=1000,M_y=100$.

Table 4. Maximum $\tau$ for $0\le I\left(0\right),R\left(0\right)\le 100$.

On the other hand, the boundary problem (10) - (11) has three independent variables which further complicates its numerical study, although the technique is similar. Let us multiply multiple (10) by a regular function $\varphi \left(x\mathrm{,}y\mathrm{,}z\right)$ satisfying the homogeneous Dirichlet boundary conditions on ${\Gamma }_1$. Integrating over the domain $\Omega =\left[0,M_x\right]\times \left[0,M_y\right]\times \left[0,M_z\right]$, the following terms with will appear:

${\displaystyle {\int }_{\Omega }}d\left(x,y\right)\frac{{\partial }^2\tau }{\partial x^2}\varphi =-{\displaystyle {\int }_{\Omega }}d\frac{\partial \tau }{\partial x}\frac{\partial \varphi }{\partial x}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-{\displaystyle {\int }_{\Omega }}\frac{\partial d}{\partial x}\frac{\partial \tau }{\partial x}\varphi \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+{\displaystyle {\int }_{\partial \Omega }}d\left(x,y\right)\varphi \left(\frac{\partial \tau }{\partial x}{n}_x\right),$

${\displaystyle {\int }_{\Omega }}c\left(x,y\right)\frac{{\partial }^2\tau }{\partial y^2}\varphi =-{\displaystyle {\int }_{\Omega }}c\frac{\partial \tau }{\partial y}\frac{\partial \varphi }{\partial y}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-{\displaystyle {\int }_{\Omega }}\frac{\partial c}{\partial y}\frac{\partial \tau }{\partial y}\varphi +{\displaystyle {\int }_{\partial \Omega }}c\left(x,y\right)\varphi \left(\frac{\partial \tau }{\partial y}{n}_y\right),$

${\displaystyle {\int }_{\Omega }}\gamma y\frac{{\partial }^2\tau }{\partial z^2}\varphi =-{\displaystyle {\int }_{\Omega }}\gamma y\frac{\partial \tau }{\partial z}\frac{\partial \varphi }{\partial z}+{\displaystyle {\int }_{\partial \Omega }}\gamma y\varphi \left(\frac{\partial \tau }{\partial z}{n}_z\right),$

$-{\displaystyle {\int }_{\Omega }}d\left(x,y\right)\frac{{\partial }^2\tau }{\partial x\partial y}\varphi ={\displaystyle {\int }_{\Omega }}d\frac{\partial \tau }{\partial y}\frac{\partial \varphi }{\partial x}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+{\displaystyle {\int }_{\Omega }}\frac{\partial d}{\partial x}\frac{\partial \tau }{\partial y}\varphi \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-{\displaystyle {\int }_{\partial \Omega }}d\varphi \left(\frac{\partial \tau }{\partial y}{n}_x\right),$

$-{\displaystyle {\int }_{\Omega }}\gamma y\frac{{\partial }^2\tau }{\partial y\partial z}\varphi ={\displaystyle {\int }_{\Omega }}\gamma y\frac{\partial \tau }{\partial z}\frac{\partial \varphi }{\partial y}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+{\displaystyle {\int }_{\Omega }}\gamma \frac{\partial \tau }{\partial z}\varphi \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}-{\displaystyle {\int }_{\partial \Omega }}\gamma y\varphi \left(\frac{\partial \tau }{\partial z}{n}_y\right),$

and $c\left(x\mathrm{,}y\right)=d\left(x\mathrm{,}y\right)+\left(\gamma +\beta \right)y$ and where the boundary terms are of the following form:

$\begin{array}{l}{\displaystyle {\int }_{\partial \Omega }}d\varphi \left(\frac{\partial \tau }{\partial x}{n}_x\right)={\displaystyle {\int }_{\partial \Omega }}c\varphi \left(\frac{\partial \tau }{\partial y}{n}_y\right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}=\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\displaystyle {\int }_{\partial \Omega }}\gamma y\varphi \left(\frac{\partial \tau }{\partial z}{n}_z\right)=0,\\ {\displaystyle {\int }_{\partial \Omega }}d\varphi \left(\frac{\partial \tau }{\partial y}{n}_x\right)={\displaystyle {\int }_{{\Gamma }_2}}d\varphi \frac{\partial \tau }{\partial y},\\ {\displaystyle {\int }_{\partial \Omega }}\gamma y\varphi \left(\frac{\partial \tau }{\partial z}{n}_y\right)={\displaystyle {\int }_{{\Gamma }_3}}\gamma y\varphi \frac{\partial \tau }{\partial z}.\end{array}$.

In Figure 3 we can see the numerical solution (10) - (11) in $S\left(0\right)=x=1000$ and $M_y=100,M_z=100$, again the three graphs are very similar but the scale is quite different and the value of its peaks can be read in Table 4. Note that the value of its peaks are close to half of the SIS model, it now appears that the distance between the two models is greater than the deterministic case.

3.3. Numerical Simulation of the Stochastic Models

In view of the results of the two previous subsections: the great difference in extinction time between deterministic and stochastic models, we thought it reasonable to use some numerical simulation of stochastic differential equations in order to have more information. Our numerical simulations for the two stochastic

Figure 3. The numerical solution of (10)-(11) with My=100, Mz=100.

models have been done using the classic Euler-Maruyama numerical method, although it has strong order 1/2 and weak order 1 (we refer to [24] or [17] for a review on the numerical solutions of SDEs). This method applied to stochastic SIS model is straightforward and may be described as follows:

Algorithm 3.1. Given $\alpha \mathrm{,}\beta$ and $\gamma$ , let $\Delta t$ be the time-step, nrun, the number of simulations, and Tol the tolerances. Then, we have:

This Algorithm obtains the mean and standard deviation of the stopping time $T^{\mathrm{<mo>\; *\; </mo>}}=inf\left\{I_n\text{or}{S}_n\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}<Tol\right\}$. The numerical results are in Table 5 with $S\left(0\right)=1000$, $I\left(0\right)=10$, $\Delta t={10}^{-3}$, $nrum=100$. Note first that deviations from mean are very high, to obtain conclusions it would be necessary to decrease the step and increase the number of trials, which would cost a lot of computing time with doubtful results due to possible on numerical stability problems. In other words, we needed more precise and reliable numerical methods. However, note that the times are less than the estimates in the previous subsection and much further away from the deterministic problem. For the SIR model the results are very similar.

Table 5. Euler-Maruyama for $S\left(0\right)=1000,I\left(0\right)=10$.

4. Conclusions

In summary, in this paper first we have found the stochastic SIS and SIR models with death and computed the mean of the extinction-time analyzing and simulating its backward Kolmogorov differential equation. The more important conclusions are the following:

1) Although the reproduction number concludes the extinction or not of the disease, however, this number it does not help to know its extinction times because example with the same reproduction numbers has quite different time.

2) The cases that increase the value of the first parameter $\alpha$ also disappear more quickly, although always with the reproduction number ${\Re }_0<1$.

3) Infection disappears much earlier in stochastic models than in the corresponding deterministic models.

4) Finally, we again encounter a stochastic model, as we saw in [18], which is verydifferent from the deterministic model.

Perhaps the most important problem now is to describe techniques to computer the parameter values using data set, i.e. a calibration of the model, in [17] these authors have researched an academic example.

Acknowledgements

This work was supported by Spanish Ministry of Sciences Innovation and Universities with the project PGC2018-094522-B-100 and by the Basque Government with the project IT1247-19. The authors would like to thank the referees for the constructive and helpful comments and suggestions on the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Brockman, J. (1995) The Third Culture: Beyond the Scientific Revolution. Simon and Shuster, New York.
[2]	Murray, J.D. (2002) Mathematical Biology I: An Introduction. 3rd Edition, Springer-Verlag, New York. https://doi.org/10.1007/b98868
[3]	Daley, D.J. and Gani, J. (1999) Epidemid Modelling: An Introduction. Cambridge University Press, Cambridge.
[4]	Brauer, F. and Castillo-Chaves, C. (2012) Mathematical Models in Population Biology and Epidemiology. Springer, New York.
[5]	Brauer, F., Castillo-Chaves, C. and Feng, Z. (2019) Mathematical Models in Epidemiology. Springer, New York. https://doi.org/10.1007/978-1-4939-9828-9
[6]	Britton, T. (2010) Stochastic Epidemic Model: A Survey. Mathematical Biosciences, 225, 24-35. https://doi.org/10.1016/j.mbs.2010.01.006
[7]	Jiang, D., Yu, J., Ji, C. and Shi, N. (2011) Asymptotic Behavior of Global Positive Solution to a Stochastic SIR Model. Mathematical and Computer Modelling, 54, 221-232. https://doi.org/10.1016/j.mcm.2011.02.004
[8]	Zhao, Y. and Jiang, D. (2014) The Threshold of a Stochastic SIRS Epidemic Model with Saturated Incidence. Applied Mathematics Letter, 34, 90-93. https://doi.org/10.1016/j.aml.2013.11.002
[9]	Rajasekar, S.P. and Pitchaimani, M. (2020) Ergodic Stationary Distribution and Extinction of a Stochastic SIRS Epidemic Model with Logistic Growth and Nonlinear Incidence. Applied Mathematics Letter and Computation, 337, Article ID: 125143. https://doi.org/10.1016/j.amc.2020.125143
[10]	Liu, Q. and Jiang, D. (2020) Threshold Behavior in a Stochastic SIR Epidemic Model with Logistic Birth. Physica A: Statistical Mechanics and its Applications, 540, Article ID: 123488. https://doi.org/10.1016/j.physa.2019.123488
[11]	Allen, E. (2007) Modeling with Itô Stochastic Differential Equations. Vol. 22, Springer, Dordrecht. https://doi.org/10.1007/978-1-4020-5953-7
[12]	Allen, L.J.S. (2010) An Introduction to Stochastic Processes with Applications to Biology. 2nd Edition, CRC Press, New York. https://doi.org/10.1201/b12537
[13]	Vadillo, F. (2019) Comparing Stochastic Lotka-Volterra Predator-Prey Models. Applied Mathematics and Computation, 360, 181-189. https://doi.org/10.1016/j.amc.2019.05.002
[14]	Gilbarg, D. and Trudinger, N. (1983) Elliptic Partial Differential Equations of Second Order. 2nd Edition, Vol. 224, Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-61798-0
[15]	de la Hoz, F. and Vadillo, F. (2012) A Mean Extinction-Time Estimate for a Stochastic Lotka-Volterra Predator-Prey Model. Applied Mathematics and Computation, 219, 170-179. https://doi.org/10.1016/j.amc.2012.05.060
[16]	de la Hoz, F., Doubova, A. and Vadillo, F. (2015) Persistence-Time Estimation for Some Stochastic SIS Epidemic Models. Discrete and Continuous Dynamical Systems Series B, 20, 2933-2947. https://doi.org/10.3934/dcdsb.2015.20.2933
[17]	Vadillo, F. (2020) On a Stochastic SIS Epidemic Model with Demography and its Calibration. Submitted.
[18]	Vadillo, F. (2020) On Multiple Pathogen Models: Deterministic or Stochastic? Submitted.
[19]	Hecht, F. (2012) New Development in Freefem++. Journal of Numerical Mathematics, 20, 251-265. https://doi.org/10.1515/jnum-2012-0013
[20]	Nasell, I. (2002) Stochastic Models for Some Epidemical Infections. Mathematical Biosciences, 179, 1-19. https://doi.org/10.1016/S0025-5564(02)00098-6
[21]	Greenhalgh, D., Liang, Y. and Mao, X. (2016) SDE SIS Epidemic Model with Demographic Stochasticity and Varying Populationsize. Applied Mathematics and Computation, 276, 218-238. https://doi.org/10.1016/j.amc.2015.11.094
[22]	Allen, L.J.S. and Burgin, A.M. (2000) Comparison of Deterministic and Stochastic SIS and SIR Models in Discrete Time. Mathematical Biosciences, 163, 1-33. https://doi.org/10.1016/S0025-5564(99)00047-4
[23]	Allen, L.J.S. (2017) A Primer on Stochastic Epidemic Models: Formulation, Numerical Simulation, and Analysis. Infectious Disease Modeling, 2, 128-142. https://doi.org/10.1016/j.idm.2017.03.001
[24]	Kloeden, P.E. and Platen, E. (1998) Numerical Solution of Stochastic Differential Equations. Cambridge University Press, Cambridge.