On Deterministic and Stochastic Multiple Pathogen Epidemic Models

Abstract

In this paper, we consider a stochastic epidemic model with two pathogens. In order to analyze the coexistence of two pathogens, we compute numerically the expectation time until extinction (the mean persistence time), which satisfies a stationary partial differential equation with degenerate variable coefficients, related to backward Kolmogorov equation. I use the finite element method in order to solve this equation, and we implement it in FreeFem++. The main conclusion of this paper is that the deterministic and stochastic epidemic models differ considerably in predicting coexistence of the two diseases and in the extinction outcome of one of them. Now, the main challenge would be to find an explanation for this result.

1. Introduction

There have been numerous models for the spread of infectious diseases in populations. In the deterministic models, you may consult some excellent books on this subject, from the classic Chapter 10 in [1,2,3], and other excellent and complete references such as Chapters 9 and 10 in [4] and the monographs in [5,6]. All these deterministic models serve as a framework for formulating analogous stochastic models; these models are characterized by randomness and the variables are the solutions of stochastic differential equation (SDE). In this way, in my opinion there have been two ways to pass from the ODE to the SDE, one of them is simply adding new stochastic terms, for example, in [7,8,9,10]. The second 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 get then a stochastic system. It has recently been shown in [13,14] that this difference in the stochastic part cause great differences in the asymptotic behavior of the solutions.

In this paper, we will consider a stochastic epidemic model with two pathogens. There are several studies on the evolution and the dynamics of deterministic epidemic models with multiple pathogens (see for example [15] and references therein). Later, the authors of [16,17] proposed a stochastic model, and their main conclusions are that both models differ considerably in predicting coexistence of the two pathogens, and the coexistence in the stochastic model has a very low probability. In this paper, we analyze the mean persistence time for this stochastic model solving numerically the backward Kolmogorov equation. In order to solve this equation numerically, we will use the Finite Element Method (FEM). This authors have studied this kind of problem in several papers [13,18,19] with very hopeful results to spread more complex problems.

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

$$d\mathbf{X}(t)=\mu (t,\mathbf{X}(t))dt+\mathbf{B}(t,\mathbf{X}(t))d\mathbf{W}(t),$$

(1)

where $\mathbf{X}={(X_1,\cdots ,X_n)}^T$ and $\mathbf{W}={(W_1,\cdots ,W_m)}^T$ are m independent Wiener processes. The vectorial function $\mu (t,\mathbf{X}(t))$ is called the drift, and $\mathbf{B}(t,\mathbf{X}(t))$ is the diffusion matrix, a matrix $n\times m$.

Obviously, a key question in epidemic models is to understand the constraints that lead to extinction or persistence of the disease. In order to study this question, let us 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\{t\ge 0:X_j=0,\mathrm{for}\mathrm{some}j=1,\cdots ,n\},$$

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

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

$$\mathfrak{L}(\tau )\equiv \sum _{j=1}^n{\mu }_j{\displaystyle \frac{\partial \tau }{x_j}}+{\displaystyle \frac{1}{2}}\sum _{j=1}^n\sum _{k=1}^n{d}_{jk}{\displaystyle \frac{{\partial }^2\tau }{\partial x_j\partial x_k}}=-1.$$

(2)

where $\mathbf{D}=\mathbf{B}{\mathbf{B}}^T$ and with boundary conditions

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

(3)

assuming that $x_j\le M_j,j=1,\cdots ,n$.

The Equation (2) is an elliptic partial differential equations of second order [20] or Chapter 8 in [21,22], really it is an advection–diffusion equation, and as the name suggests, the mean persistence time will depend on the operator $\mathbf{D}$. These comments would explain the results in [13]. Moreover, in [7,8,9,10], the matrices $\mathbf{D}$ are diagonals which implies that the variables are not correlated, maybe an unrealistic hypothesis.

The paper is organized as follows. In Section 2.1, first we will present a stochastic epidemic model with two pathogens. Then, using the Symbolic Math Toolbox of Matlab©, we compute the equilibrium states for the deterministic part of the model and its numerical simulations by the classical Euler–Maruyama method. In Section 2.4, we will describe this techniques, based on the resolution of an associated backward Kolmogorov type equation for the mean persistence time, $\tau$. In order to do the FreeFem++ implementation, we will first write a well-suited variational formulation and then we will present the numerical results. In Section 3, we compare the dynamics of the two models in tree examples, and finally in Section 4, we draw the main conclusions and some future researches.

Our numerical methods were implemented in Matlab© and FreeFem++, which are freely available and particularly efficient, see in [23]. 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. Materials and Methods

2.1. Derivation of Stochastic Epidemic Model with Two Pathogen Strains

In this section, we will present an epidemic model with two pathogen strains and random demographics. The changes and their probabilities to the first order in $▵t$ are given in

Table 1. Possible change in $\mathbf{X}={(S,I_1,I_2)}^T$ and their probabilities.

Changes	Probabilities
$▵{\mathbf{X}}^{(1)}={[1,0,0]}^T$	$p_1=bS▵t$
$▵{\mathbf{X}}^{(2)}={[0,1,0]}^T$	$p_2=b{I}_1▵t$
$▵{\mathbf{X}}^{(3)}={[0,0,1]}^T$	$p_3=b{I}_2▵t$
$▵{\mathbf{X}}^{(4)}={[-1,0,0]}^T$	$p_4=d(N)S▵t$
$▵{\mathbf{X}}^{(5)}={[0,-1,0]}^T$	$p_5=({\alpha }_1+d(N))I_1▵t$
$▵{\mathbf{X}}^{(6)}={[0,0,-1]}^T$	$p_6=({\alpha }_2+d(N))I_2▵t$
$▵{\mathbf{X}}^{(7)}={[1,-1,0]}^T$	$p_7=b_1{I}_1▵t$
$▵{\mathbf{X}}^{(8)}={[1,0,-1]}^T$	$p_8=b_2{I}_2▵t$
$▵{\mathbf{X}}^{(9)}={[-1,1,0]}^T$	$p_9={\beta }_{1}S{I}_1/N▵t$
$▵{\mathbf{X}}^{(10)}={[-1,0,1]}^T$	$p_{10}={\beta }_{2}S{I}_2/N▵t$
$▵{\mathbf{X}}^{(11)}={[0,0,0]}^T$	$p_{11}=1-{\sum }_{i=1}^{10}{p}_i$

with $\mathbf{x}={(S,I_1,I_2)}^T$. In this model, it is assumed that the infection with on strain immunizes for another disease.

Here, $b>0$ is the per capita birth rate and $d(N)=1+0.05N$ is the per capita death rate, depending of the density $N(t)=S(t)+I_1(t)+I_2(t)$. Parameter ${\beta }_j>0$ is the transmission rate and ${\alpha }_j>0$ is the disease-related per capita death rate for individuals infected with strain $I_j$. Moreover, the per capital birth rate is divided into two parts: $b_j$ and $b-b_j$, if the strain j is passed from mother to offspring (vertical transmission), the per capital birth rate to the infected is $b-b_j$, in other case, if there is no vertical transmission, the newborn enter the susceptible class and $b_j=b$.

Fixing $\mathbf{X}(t)$ at time t, we calculate the expected change for the change $\mathbf{X}={(S,I_1,I_2)}^T$

$$\mathrm{E}(▵\mathbf{X})=\sum _{j=1}^{11}{p}_j▵{\mathbf{X}}^{(j)}=\left(\begin{array}{c}{\mu }_1(S,I_1,I_2)\\ {\mu }_2(S,I_1,I_2)\\ {\mu }_3(S,I_1,I_2)\end{array}\right)▵t,$$

(4)

where

$$\left\{\begin{array}{c}{\mu }_1(S,I_1,I_2)=S\left(b-d(N)-{\displaystyle \frac{{\beta }_1{I}_1+{\beta }_2{I}_2}{N}}\right)+b_1{I}_1+b_2{I}_2,\hfill \\ {\mu }_2(S,I_1,I_2)=I_1\left(b-b_1-d(N)-{\alpha }_1+{\displaystyle \frac{{\beta }_{1}S}{N}}\right),\hfill \\ {\mu }_3(S,I_1,I_2)=I_2\left(b-b_2-d(N)-{\alpha }_2+{\displaystyle \frac{{\beta }_{2}S}{N}}\right).\hfill \end{array}\right.$$

(5)

and the covariance matrix

$$\mathrm{E}(▵\mathbf{X}{(▵\mathbf{X})}^T)=\sum _{j=1}^{11}{p}_j(▵{\mathbf{X}}^{(j)}){(▵{\mathbf{X}}^{(j)})}^T=\mathbf{D}(S,I_1,I_2)▵t,$$

(6)

where $\mathbf{D}=(d_{i,j})$ is the diffusion matrix, a matrix $3\times 3$ with $d_{2,3}=d_{3,2}=0$ and

$$\begin{array}{cc}\hfill d_{1,1}& =S\left(b+d(N)+\frac{{\beta }_1{I}_1+{\beta }_2{I}_2}{N}\right)+b_1{I}_1+b_2{I}_2,\hfill \\ \hfill d_{j+1,j+1}& =I_j\left(b+b_j+d(N)+{\alpha }_j+\frac{{\beta }_{j}S}{N}\right),j=1,2,\hfill \\ \hfill d_{1,j+1}& =d_{j+1,1}=-b_j{I}_j-{\beta }_j\frac{S{I}_j}{N},j=1,2.\hfill \end{array}$$

Finally, the stochastic differential system (SDS) is

$$d\mathbf{X}(t)=\mu (t,\mathbf{X}(t))dt+{\mathbf{D}}^{1/2}(t,\mathbf{X}(t))d\mathbf{W}(t).$$

(7)

2.2. The Deterministic Model

The deterministic part is the following ordinary differential system (ODEs)q:

$$\left\{\begin{array}{c}{\displaystyle \frac{dS}{dt}}={\mu }_1(S,I_1,I_2),\hfill \\ {\displaystyle \frac{d{I}_1}{dt}}={\mu }_2(S,I_1,I_2),\hfill \\ {\displaystyle \frac{d{I}_2}{dt}}={\mu }_3(S,I_1,I_2).\hfill \end{array}\right.$$

(8)

There are numerous theoretical studies on the evolution, persistence, or extinction of multiple pathogen strains for the deterministic epidemic models, see, for example, in [16,17] and their references. It has been proved that the dynamics of a deterministic model depends on the basic reproduction numbers defined by

$$R_j=\frac{{\beta }_j+b-b_j}{b+{\alpha }_j},j=1,2,$$

and it is known that if $R_j<1$, then the disease extinction occurs. Using the Symbolic Math Toolbox of Matlab©, we can compute six equilibrium values of the deterministic model which can be found in the

Table 2. Zeros of (8): each column corresponds to a zero.

Zeros	1	2	3	4	5	6
S	∗	∗	∗	0	∗	0
$I_1$	∗	0	0	∗	∗	∗
$I_2$	0	∗	∗	∗	∗	∗

, each column of this table is a zero of (8). We have denoted by ∗ nonzero values, whose exact formulas we did not write for their long expressions and lack of interest. Obviously, the asymptotic behavior of each equilibrium depends on the eigenvalues of the Jacobian.

2.3. Simulation Using the Euler–Maruyama Method

The numerical simulation for the stochastic differential system (7) implements the classic Euler–Maruyama numerical method, although it has strong order 1/2 and weak order 1 (see for example [24] or ([25] and more recently [26]). This method is simple and straightforward to implement [27].

We have written Euler–Maruyama algorithm, similar to an algorithm from in [18,19], with three stopping test, one for each of the variables S, $I_1$, and $I_2$, and with $▵t$ the time step. Essentially, given $▵t$, the number of simulation and an initial position $(S(0),I_1(0),I_2(0))$, the algorithm is straightforward until one of the variable is less than one, i.e., $S_n<1$ or ${(I_1)}_n<1$ or ${(I_2)}_n<1$. After all the trials, we computed the mean and standard deviation of stopping times. Matlab implementation is very similar to the program offered in the appendix of [18]. It is important to remark that in this implementation, at each step we have to compute the matrix ${\mathbf{D}}^{1/2}$ using the Matlab’ command sqrtm.

2.4. The Mean Persistence Time for the Model

In this section, in order to predict the behavior of two pathogens in the stochastic model (7), we present an alternative technique based on the estimation of the mean existence time $\tau$, by numerical resolution of the Kolmogorov type equation using the FEM and the FreeFem++ implementation.

Let us denote by $s=S(0)$, $y_1=I_1(0)$, $y_2=I_2(0)$ and $\mathsf{\Omega }=[0,M_s]\times [0,M_1]\times [0,M_2]$, where $M_s,M_1$ and $M_2$ are positive constants and its boundary by $\partial \mathsf{\Omega }={\mathsf{\Gamma }}_1\cup {\mathsf{\Gamma }}_2\cup {\mathsf{\Gamma }}_3\cup {\mathsf{\Gamma }}_4$, where

$$\left\{\begin{array}{c}{\mathsf{\Gamma }}_1=\{(s,y_1,y_2)\in \mathsf{\Omega }|s=0,\mathrm{or}{y}_1=0,\mathrm{or}{y}_2=0\},\hfill \\ {\mathsf{\Gamma }}_2=\{(s,y_1,y_2)\in \mathsf{\Omega }|s=M_s\},\hfill \\ {\mathsf{\Gamma }}_3=\{(s,y_1,y_2)\in \mathsf{\Omega }|y_1=M_1\},\hfill \\ {\mathsf{\Gamma }}_4=\{(s,y_1,y_2)\in \mathsf{\Omega }|y_2=M_2\}.\hfill \end{array}\right.$$

The mean persistence time $\tau =\tau (s,y_1,y_2)$ for the stochastic model (7) satisfies the following stationary partial differential equation of Kolmogorov type:

$$\begin{array}{cc}\hfill & {\mu }_1(s,y_1,y_2){\displaystyle \frac{\partial \tau }{\partial s}}+\sum _{j=1}^2{\mu }_{j+1}(s,y_1,y_2){\displaystyle \frac{\partial \tau }{\partial y_j}}\hfill \\ \hfill & \phantom{\mu (s,y_1,y_2)}+{\displaystyle \frac{1}{2}}{d}_{1,1}{\displaystyle \frac{{\partial }^2\tau }{\partial s^2}}+{\displaystyle \frac{1}{2}}\sum _{j=1}^2{d}_{j+1,j+1}{\displaystyle \frac{{\partial }^2\tau }{\partial y_j^2}}+\sum _{j=1}^2{d}_{1,j+1}{\displaystyle \frac{{\partial }^2\tau }{\partial s\partial y_j}}=-1\mathrm{in}\mathsf{\Omega }\hfill \end{array}$$

(9)

with the following boundary conditions:

$$\left\{\begin{array}{c}\tau (0,y_1,y_2)=\tau (s,0,y_2)=\tau (s,y_1,0)=0\mathrm{on}{\mathsf{\Gamma }}_1,\hfill \\ {\displaystyle \frac{\partial \tau }{\partial s}}(M_s,y_1,y_2)=0\mathrm{on}{\mathsf{\Gamma }}_2,\hfill \\ {\displaystyle \frac{\partial \tau }{\partial y_1}}(s,M_1,y_2)=0\mathrm{on}{\mathsf{\Gamma }}_3,\hfill \\ {\displaystyle \frac{\partial \tau }{\partial y_2}}(s,y_1,M_2)=0\mathrm{on}{\mathsf{\Gamma }}_4.\hfill \end{array}\right.$$

(10)

As we said at the beginning of the paper, in order to preform the FreeFem++ implementation, we have to write a variational formulation for the boundary value problem (9) and (10). The next subsection will be devoted to do this.

2.5. Variational Formulation

In order to perform numerical experiments using FreeFem++, which is a partial differential equation solver and has its own language, we will write (9) and (10) in variational form. For this, we will make formal computations. More precisely, let us write (9) as follows:

$$\begin{array}{cc}\hfill -{\mu }_1\frac{\partial \tau }{\partial s}-{\mu }_2\frac{\partial \tau }{\partial y_1}& -{\mu }_3\frac{\partial \tau }{\partial y_2}-\frac{1}{2}{d}_{1,1}\frac{{\partial }^2\tau }{\partial s^2}-\frac{1}{2}{d}_{2,2}\frac{{\partial }^2\tau }{\partial y_1^2}\hfill \\ \hfill & -\frac{1}{2}{d}_{3,3}\frac{\partial \tau }{{\partial }^2{y}_2^2}-d_{1,2}\frac{{\partial }^2\tau }{\partial s\partial y_1}-d_{1,3}\frac{{\partial }^2\tau }{\partial s\partial y_2}=1\mathrm{in}\mathsf{\Omega }.\hfill \end{array}$$

(11)

Let us multiply multiple (11) by a regular “test” function $\varphi (s,y_1,y_2)$ satisfying the homogeneous Dirichlet boundary conditions on ${\mathsf{\Gamma }}_1$. Integrating over the domain $\mathsf{\Omega }$, the following terms with will appear:

$$\begin{array}{ccc}\hfill -\frac{1}{2}{\int }_{\mathsf{\Omega }}{d}_{1,1}\frac{{\partial }^2\tau }{\partial s^2}\varphi & =& \frac{1}{2}{\int }_{\mathsf{\Omega }}{d}_{1,1}\frac{\partial \tau }{\partial s}\frac{\partial \varphi }{\partial s}+\frac{1}{2}{\int }_{\mathsf{\Omega }}\frac{\partial d_{1,1}}{\partial s}\frac{\partial \tau }{\partial s}\varphi -\frac{1}{2}{\int }_{\partial \mathsf{\Omega }}{d}_{1,1}\varphi \left(\frac{\partial \tau }{\partial s}{n}_s\right),\hfill \\ \hfill -\frac{1}{2}{\int }_{\mathsf{\Omega }}{d}_{2,2}\frac{{\partial }^2\tau }{\partial y_1^2}\varphi & =& \frac{1}{2}{\int }_{\mathsf{\Omega }}{d}_{2,2}\frac{\partial \tau }{\partial y_1}\frac{\partial \varphi }{\partial y_1}+\frac{1}{2}{\int }_{\mathsf{\Omega }}\frac{\partial d_{2,2}}{\partial y_1}\frac{\partial \tau }{\partial y_1}\varphi -\frac{1}{2}{\int }_{\partial \mathsf{\Omega }}{d}_{2,2}\varphi \left(\frac{\partial \tau }{\partial y_1}{n}_{y_1}\right),\hfill \\ \hfill -\frac{1}{2}{\int }_{\mathsf{\Omega }}{d}_{3,3}\frac{{\partial }^2\tau }{\partial y_2^2}\varphi & =& \frac{1}{2}{\int }_{\mathsf{\Omega }}{d}_{3,3}\frac{\partial \tau }{\partial y_2}\frac{\partial \varphi }{\partial y_2}+\frac{1}{2}{\int }_{\mathsf{\Omega }}\frac{\partial d_{3,3}}{\partial y_2}\frac{\partial \tau }{\partial y_2}\varphi -\frac{1}{2}{\int }_{\partial \mathsf{\Omega }}{d}_{3,3}\varphi \left(\frac{\partial \tau }{\partial y_2}{n}_{y_2}\right),\hfill \\ \hfill -{\int }_{\mathsf{\Omega }}{d}_{1,2}\frac{{\partial }^2\tau }{\partial s\partial y_1}\varphi & =& {\int }_{\mathsf{\Omega }}{d}_{1,2}\frac{\partial \tau }{\partial y_1}\frac{\partial \varphi }{\partial s}+{\int }_{\mathsf{\Omega }}\frac{\partial d_{1,2}}{\partial s}\frac{\partial \tau }{\partial y_1}\varphi -{\int }_{\partial \mathsf{\Omega }}{d}_{1,2}\varphi \left(\frac{\partial \tau }{\partial y_1}{n}_s\right),\hfill \\ \hfill -{\int }_{\mathsf{\Omega }}{d}_{1,3}\frac{{\partial }^2\tau }{\partial s\partial y_2}\varphi & =& {\int }_{\mathsf{\Omega }}{d}_{1,3}\frac{\partial \tau }{\partial y_2}\frac{\partial \varphi }{\partial s}+{\int }_{\mathsf{\Omega }}\frac{\partial d_{1,3}}{\partial s}\frac{\partial \tau }{\partial y_2}\varphi -{\int }_{\partial \mathsf{\Omega }}{d}_{1,3}\varphi \left(\frac{\partial \tau }{\partial y_2}{n}_s\right),\hfill \end{array}$$

where the boundary terms are of the following form:

$$\left\{\begin{array}{ccc}{\displaystyle {\int }_{\partial \mathsf{\Omega }}{d}_{1,1}\varphi \left(\frac{\partial \tau }{\partial s}{n}_s\right)}\hfill & =\hfill & {\displaystyle {\int }_{\partial \mathsf{\Omega }}{d}_{2,2}\varphi \left(\frac{\partial \tau }{\partial y_1}{n}_{y_1}\right)={\int }_{\partial \mathsf{\Omega }}{d}_{3,3}\varphi \left(\frac{\partial \tau }{\partial y_2}{n}_{y_2}\right)=0,}\hfill \\ {\displaystyle {\int }_{\partial \mathsf{\Omega }}{d}_{1,2}\varphi \left(\frac{\partial \tau }{\partial y_1}{n}_s\right)}\hfill & =\hfill & {\displaystyle {\int }_{{\mathsf{\Gamma }}_2}{d}_{1,2}\varphi \frac{\partial \tau }{\partial y_1},}\hfill \\ {\displaystyle {\int }_{\partial \mathsf{\Omega }}{d}_{1,3}\varphi \left(\frac{\partial \tau }{\partial y_2}{n}_s\right)}\hfill & =\hfill & {\displaystyle {\int }_{{\mathsf{\Gamma }}_2}{d}_{1,3}\varphi \frac{\partial \tau }{\partial y_2},}\hfill \end{array}\right.$$

with $n_s$, $n_{y_1}$ and $n_{y_2}$ are the normal vectors exterior to the boundary ${\mathsf{\Gamma }}_2$, ${\mathsf{\Gamma }}_3$ and ${\mathsf{\Gamma }}_4$, respectively.

3. Numerical Results

We present three numerical examples for the deterministic and stochastic epidemic models, it is remarkable that in Examples 2 and 3 the coexistence dynamics differ between the deterministic and stochastic models.

Example 1.

In this first example, the per capital death and birth rates are$d(N)=1+N/100$and$b=2$. In addition,${\beta }_1=7,{\beta }_2=5.25,{\alpha }_1={\alpha }_2=0.5$and$b_1=b_2=3$, so that$R_1=2;R_2=1.5$and according Theorem 2 in [16] the solutions to deterministic model (8) convergent to

$$\underset{t\to \infty }{lim}{I}_1(t)=39.35,\underset{t\to \infty }{lim}{I}_2(t)=0.$$

In Figure 1, we have plotted the solution of (8) usingMatlab’command ode45 for $S(0)=1000,I_1(0)=50$, and $I_2(0)=50$ in $0\le t\le 7$. The final solutions are $S(7)\approx 33.9912$, $I_1(7)\approx 39.5387$, and $I_2(7)\approx 0.0712$.

On the other hand, in Figure 2 we have represented the numerical solution of the problem (9) and (10) with $M_s=1000,M_1=M_2=100$. More precisely, we have represented $\tau (1000,I_1(0),I_2(0)$ for $0\le I_1(0),I_2(0)\le 100$, also we have highlighted that $\tau (1000,50,50)\approx 2.717$.

In

Table 3. Example 1: Euler–Maruyama Method.

Initial Point		Number of Stops	Mean	Std
(1000, 50, 50)	$S<1$	0
	$I_1<1$	120	1.6929	0.8579
	$I_2<1$	9880	2.5880	1.0005

, we show the results of 10,000 trials of the Euler–Maruyama method with $▵t={10}^{-4}$, and the value is represented in Figure 2. For each initial population size, we have written down the number of stops when $S<1$ or $I_1<1$ or $I_2<1$, and we have computed the mean stop time (mean) and its standard deviation (std). This result, which is close to the figure’s estimate, and it does not seem too far from the deterministic model, to get an idea in the deterministic problem $S(3)\approx 38.87,I_1(3)\approx 42.12,I_2(3)\approx 1.987$.

Example 2.

Vertical transmission of both strains.Let us suppose now that both strains are transmitted vertically. This corresponds to take$b_1=b_2=0$ and $b=6$, ${\beta }_1=15,{\beta }_2=1$, ${\alpha }_1=2.5,{\alpha }_2=2$. In this case, the basic reproduction numbers are

$$R_1=2.4706>1,R_2=0.8750<1.$$

The equilibrium point is $E=(2.1684,4.3367,54.2092)$, and the eigenvalues of the Jacobian are $-{10}^{-5}\pm 0.6682i$ and $-3.0357$. This means that the deterministic solution cycles closer and closer to each equilibriums as we can see in Figure 5 from [17], however its size fluctuates greatly.

In Figure 3, we have represented the numerical solution of the problem (9) and (10) with $M_s=1000,M_1=M_2=100$. More precisely, we have represented $\tau (1000,I_1(0),I_2(0))$ for $0\le I_1(0),I_2(0)\le 100$, also we have highlighted that $\tau (1000,49,51)\approx 0.3328$.

On the other hand, in

Table 4. Example 2: Euler–Maruyama Method.

Initial Point		Number of Stops	Mean	Std
(1000, 49, 51)	$S<1$	4109	0.5716	0.0955
	$I_1<1$	14
	$I_2<1$	5877	0.3276	0.1237

we show the results of 10,000 trials of the Euler–Maruyama method with $▵t={10}^{-4}$ and the value is represented in Figure 3. For each initial population size, we have written down the number of stops when $S<1$ or $I_1<1$ or $I_2<1$, and we have computed the mean stop time (mean) and its standard deviation (std). This result, which is close to the figure’s estimate, apparently differs from the deterministic solution, but not so much because, for example, with these initial values, the numerical solution is

$$S(0.8186)\approx 0.0381,I_1(0.8186)\approx 56.3229,I_2(0.8186)\approx 2.8159,$$

and therefore it is not strange that a slight deviation or disturbance of the trajectory ends up at zero.

Example 3.

A study of coexistence.In this example, the per capital death and birth rates are$d(N)=1+5N/100$ and $b=6$. In addition, ${\beta }_1=30,{\beta }_2=15,{\alpha }_1=4,{\alpha }_2=5/3$, and $b_1=12,b_2=6$, so that $R_1=2;R_2=1.5$, and according Theorem 2 in [16] the solutions to deterministic model (8) converge to

$$\underset{t\to \infty }{lim}{I}_1(t)=15.13,\underset{t\to \infty }{lim}{I}_2(t)=22.35.$$

In Figure 4, we have plotted the solution of (8) using Matlab for $S(0)=1000$, $I_1(0)=52$ and $I_2(0)=51$ in $0\le t\le 20$. The final solutions are $S(20)\approx 35.7715$, $I_1(20)\approx 15.1724$, and $I_2(30)\approx 22.1665$.

In Figure 5, we have represented the numerical solution of the problem (9) and (10) with $M_s=1000,M_1=M_2=100$. More precisely, we have represented $\tau (1000,I_1(0),I_2(0))$ for $0\le I_1(0),I_2(0)\le 100$, also we have highlighted that $\tau (1000,52,51)\approx 0.5838$.

In

Table 5. Example 3: Euler–Maruyama Method.

Initial Point		Number of Stops	Mean	Std
(1000, 52, 51)	$S<1$	0
	$I_1<1$	9543	0.9348	0.3894
	$I_2<1$	457	0.5024	0.3119

, we show the results of 10,000 trials of the Euler–Maruyama method with $▵t={10}^{-4}$ and the value representing in Figure 5. For each initial population size, we have written down the number of stops when $S<1$ or $I_1<1$ or $I_2<1$, and we have computed the mean stop time (mean) and its standard deviation (std). This result, which is close to the figure’s estimate, has no relation to the behavior of the deterministic system; in this example, the asymptotic behaviors of the two models are quite different.

4. Conclusions and Discution

In this paper, we have studied the coexistence of two pathogens model proposed in [16,17]. They showed that the deterministic and stochastic models differ considerably in prediction coexistence of the two pathogens because the probability of coexistence in the stochastic model is very small. We propose an alternative strategy to analyze the behavior of the stochastic model. More precisely, we have computed the mean persistence time for the stochastic model that satisfies a partial differential Kolmogorov type equation with degenerate coefficients. In order to do this, the finite element method has been used and the numerical implementation was performed using FreeFem++.

From our example we have showed the following.

• Example 1: In this case, competitive exclusion occurs because in the deterministic model ${lim}_{t\to \infty }{I}_2(t)=0$, while the stochastic model disappears somewhat more quickly.
• Example 2: In this example, with vertical transmission of both strains, the deterministic solution cycles closer and closer while the stochastic solution is extinguished very quickly. The difference in the asymptotic behavior of deterministic and stochastic is very important.
• Example 3: The difference with respect to the previous one is that now the first infection disappears as we can see in the Figure 5 and Table 5.

The main conclusion of this paper is evident: the deterministic and stochastic epidemic models differ considerably in predicting coexistence of the two diseases and in the extinction or not of one of them. Now, the main challenge would be to find an explanation for this result, in ([4], p. 3) we can read: “If the initial population size is small then a stochastic model is more appropriate, since the likelihood that the population becomes extinct due to chance must be considered.” Deterministic models often provide useful ways of gaining sufficient understanding about the dynamics of populations whenever they are large enough. This may be true in some cases but I propose the following analysis: in the elliptical differential Equation (2) (backward Kolmogorov equation) there exists two part, the advection with the vector $\mu$ and the diffusion with the matrix $\mathbf{D}$, the first tends to move the solution while the second wears it down, in another words: which dominates, the advection or the diffusion? Let is define the ratio

$$R(t)={\displaystyle \frac{\left|\right|\mu (S(t),I_1(t),I_2(t)){\left|\right|}_2}{\left|\right|\mathbf{D}(S(t),I_1(t),I_2(t)){\left|\right|}_2}},$$

(12)

comparing the evolution of these two terms of the equation.

In Figure 6, we have plotted the evolution of $R_j$ using Euler–Maruyama in (7) with $▵t={10}^{-4}$ and $S(0)=1000,I_1(0)=50$ and $I_2(0)=50$. Clearly the red trajectory (Example 1) decreases more slowly than the other two (Examples 2 and 3), this could explain its delay in extinction time. Obviously this does not prove anything, it only indicates a possible line of future research.

Finally, in my opinion we would need to make a previous estimate of the mean persistence time to fully understand the dynamics of a complete stochastic model. In my opinion, the stochastic model seems more realistic because, although it starts out swinging, it disappears reasonably after some time. The environment is constantly evolving, and as the philosopher Heraclitus wrote over 25 centuries ago: Everything changes and nothing remains still … and … you cannot step twice into the same stream.