Extinction-time for stochastic population models

Journal of Computational and Applied Mathematics 295 (2016) 159–169 Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam Extinction-time for stochastic population models Anna Doubovaa, Fernando Vadillob,∗ aDepartment of Differential Equations and Numerical Analysis, University of Sevilla, Spain bDepartment of Applied Mathematics and Statistics and Operations Research, University of the Basque Country (UPV/EHU), Spain article info Article history: Received 23 September 2014 Received in revised form 19 December 2014 MSC: 92B05 60H10 65N30 Keywords: Extinction-time Population models Stochastic differential equations Finite elements method abstract The analysis of interacting population models is the subject of much interest in mathematical ecology. Moreover, the persistence and extinction of these models is one of the most interesting and important topics, because it provides insight into their behavior. The mean extinction-time for stochastic population models considered in this paper depends on the initial population size and satisfies a stationary partial differential equation, related to the backward Kolmogorov differential equation, a linear second-order partial differential equation with variable coefficients. In this communication we review several papers where we have proposed some numerical techniques in order to estimate the mean extinction-time for stochastic population models. Besides, we will compare the theoretical predictions and numerical simulations for stochastic differential equations (SDEs). This work can be viewed as a unified review of the contributions de la Hoz and Vadillo (2012), de la Hoz et al. (2014) and Doubova and Vadillo (2014). ©2015 Elsevier B.V. All rights reserved. 1. Introduction Stochastic population models are systems of stochastic differential equations (SDE) with a deterministic part plus a stochastic one. In classic books such as for example [1] or [2], we can find results on the evolution, persistence and extinction for deterministic population models, but the stochastic models are more complicated and have quite different persistence and extinction behaviors. The SDE system for the dynamics of ninteracting populations has the following form: dX(t)=µ(t,X)dt +B(t,X)dW(t), X(0)=X0,(1) where X(t)=(X1(t), . . . , Xn(t))tis an n-dimensional random variable, and W(t)=(W1(t), . . . , Wm(t))tare mindependent Wiener processes. The vectorial function µ(t,X)is called the drift, and B(t,X)is the diffusion matrix. Notice that if the matrix B≡0, then this system reduces to a standard deterministic model for the population dynamics. One of the major goals of stochastic population dynamics (1) is to predict finite persistence time, i.e. the time when one or more components of the population go to zero. In this paper, Twill be the random variable that indicates the persistence time: T:= inf{t≥0:Xi(t)=0 for any i=1,...,n}. Notice that there are other definitions of persistence, as for example weak persistence, weak persistence in the mean, strong persistence in the mean, etc. (see [3, p. 444]). The stochastic variable Tobviously depends on the initial population size X0. ∗Corresponding author. E-mail addresses: [email protected] (A. Doubova), [email protected] (F. Vadillo). http://dx.doi.org/10.1016/j.cam.2014.12.037 0377-0427/©2015 Elsevier B.V. All rights reserved. 160 A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 On the other hand, let us also introduce the expectation time until extinction τ:= E(T). It is clear that the variable τalso depends on the initial population size, i.e. τ=τ(X0). Let us notice that this mean persistence time can be computed by solving stationary partial differential equations, which are related to the backward Kolmogorov equation (see [4, p. 443]; we refer to [5, p. 150] or Chapter 8 of [6] for a clear review of these equations). The backward Kolmogorov equations are linear second-order partial differential equations (PDEs) with variable coefficients. It is well known that it is possible to compute an analytical solution only in a few simple cases, see for example [4], for a population model with linear birth and death rates, or [7], where the exact expressions are given for the extinction-time of a class of birth–death processes. In general, we can only compute numerical approximations, see for example [8], where central difference approximations are used. In this paper we use a Finite Element Method (FEM). There is an ample bibliography about FEM, see for instance [9] or [10] for a good introduction. The more popular reference on FEM in solid mechanics are the books [11,12]. Generalization to other fields of physics or engineering have been done by applied mathematicians through the concept of variational formulation of partial differential equations (see for example [13]). This paper is organized as follows. In Section 2, we describe a Lotka–Volterra predator–prey model. We will estimate its mean extinction-time using the numerical approximations of the associated Kolmogorov equation and we will make a direct comparison between predictions and numerical simulations. In Sections 3and 4we will analyze the mean extinction-time for more complicated models. Finally in Section 5we draw the main conclusions. Our numerical methods were implemented in Matlab©and FreeFem++, the experiments were carried out in an Intel(R) Core(TM)2 Duo CPU E6850 @ 3.00 GHz. The codes for the numerical tests are available on request. 2. Extinction-time for a Lotka–Volterra predator–prey model 2.1. The Lotka–Volterra model A classical model in mathematical ecology is the standard Lotka–Volterra predator–prey model. Let us consider a simple ecosystem consisting of preys (for example rabbits) that have an infinite supply of food and predators (for example foxes) that need prey for their food. Following [6, p. 359] or [5, p. 149], this problem is modeled by a nonlinear SDE system:    dR =R(2−αF)dt +R(2+αF)dW1,t>0, dF =F(αR−1)dt +F(αR+1)dW2,t>0, R(0)=r,F(0)=f, (2) where R=R(t)and F=F(t)are random variables corresponding to the number of preys and predators, respectively, rand fare the initial population sizes, αis a positive constant and Wj(t), for j=1,2, are two independent Wiener processes. If α=0, the two populations do not interact and the foxes die off from starvation. If α > 0, the foxes encounter the rabbits with a probability that is proportional to the product of their numbers. The deterministic part is very simple: its phase plane has a center in (1/α, 2/α), and the solutions are always periodic with a period that depends on the initial conditions (see left hand-side of Fig. 1). These solutions cannot be expressed in terms of other known functions and the equations must be solved numerically. Moreover, as was already commented in [14], even tiny amounts of environmental noise cause huge differences between the stochastic model (2) and the associated deterministic model. On the right hand-side of Fig. 1, we have plotted three types of solutions corresponding to the initial values R(0)=100,F(0)=150: the deterministic solution (blue), a stochastic path for 0 ≤t≤20 (black), and the mean of hundred stochastic paths for 0 ≤t≤12 (green). The mean extinction-time τ=τ(r,f)for (2) satisfies the following stationary backward Kolmogorov equation: r(2−αf)∂τ ∂r+1 2r(2+αf)∂2τ ∂r2+f(αr−1)∂τ ∂f+1 2f(αr+1)∂2τ ∂f2= −1.(3) Taking into account that the number of preys and predators cannot exceed some positive values Mrand Mf, respectively, we have the following boundary conditions:            τ(0,f)=τ(r,0)=0, (r,f)∈(0,Mr)×(0,Mf), ∂τ ∂r(Mr,f)=0,f∈(0,Mf), ∂τ ∂f(r,Mf)=0,r∈(0,Mr). (4) We can extend this procedure to higher-order moments τp=τp(r,f). More precisely, it can be easily shown that τpsatisfy the following partial differential equation r(2−αf)∂τp ∂r+1 2r(2+αf)∂2τp ∂r2+f(αr−1)∂τp ∂f+1 2f(αr+1)∂2τp ∂f2= −pτp−1(5) with the above boundary conditions (4) written for τp. A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 161 Fig. 1. Phase plane for the deterministic part of (2) with α=0.01 (left). Three types of solutions for the initial data R(0)=100,F(0)=150 (right): the deterministic solution (blue), a stochastic path for 0 ≤t≤20 (black), and the mean of hundred stochastic paths for 0 ≤t≤12 (green). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) We have that the linear elliptic differential operator corresponding to (3) given by L=r(2−αf)∂ ∂r+1 2r(2+αf)∂2 ∂r2+f(αr−1)∂ ∂f+1 2f(αr+1)∂2 ∂f2,(6) is locally uniformly elliptic, because 1 2r(2+αf) > 0 and 1 2f(αr+1) > 0, when r,f>0 and α > 0. Then, we have the following theorem: Theorem 2.1. The solution τof the problem (3)–(4) is nonnegative. The proof of this theorem can be found in [15]. 2.2. Numerical simulation for the stochastic Lotka–Volterra model (2) Our numerical simulations for the stochastic model (2) have been done using the Euler–Maruyama method, which is a simple generalization of the Euler method for ordinary differential equation to stochastic differential equation. This method is simple and straightforward to implement. It has strong order of convergence equal to 1/2 and weak order equal 1. The traditional convergence theory for numerical methods applied to SDEs requires a global Lipschitz assumption on the drift and diffusion coefficients; however, in practice, many important SDE models satisfy only a local Lipschitz property and, since Brownian paths can make arbitrarily large excursions, the global Lipschitz based theory is not directly relevant. More details, along with further references, can be found in [16]. We refer to [17,18] for a review on the numerical solution of SDEs. We implement the Algorithm 3.1 from [15]. 2.3. Numerical methods for (3)–(4) As we have already mentioned, in order to perform the numerical experiments we will use the finite elements method. We write a well-suited variational form associated to (3)–(4) and we discretize it by means of rectangular elements with cubic polynomials as basic functions. We use a direct solver for the resulting sparse linear system, which can be formulated very conveniently in Matlab by means of sparse matrices. The complete Matlab code is offered in the Appendix of [15]. The implementation of our code for some simple examples showed us that it is convenient to choose the values of Mr and Mflarge enough depending on α. Then, bearing in mind the boundary conditions (4), we anticipate that the biggest errors should appear along r=0 and f=0, which is confirmed by the numerical simulations. Therefore, the easiest way to (partly) offset these errors is to concentrate the nodes toward the left and lower edges of the rectangle. More precisely, after testing different choices, we have found that a nodal distribution of the form ri=Mr(i/Nr)3,0≤i≤Nr, fj=Mf(j/Nf)3,0≤j≤Nf, yields particularly good results. In our experiments, for α=0.1 and α=0.05 we have taken Mr=Mf=200 and Nr=Nf=80. On the other hand, for α=0.01, the parameters are Mr=Mf=2000 and Nr=Nf=160. We have observed that larger values of Mr,Mf,Nrand Nfdo not improve our results. In Fig. 2, we have plotted the numerical approximation of mean extinction-time τand its standard deviation for α=0.1. The peak of τis located approximately at r=9.3,f=19.2, and its height is about 6.5984. As αdecreases, the peak moves toward the point (1/α, 2/α), and increases its height; for example, when α=0.01, the maximum is about 64.0532. 162 A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 8 6 4 2 0 200 150 100 50 0 050 100 150 200 f r 5 4 3 2 1 0 200 150 100 50 0050 100 150 200 f r std τ Fig. 2. Numerical solution of τand τ2−(τ)2, for α=0.1. Table 1 Comparison of Finite Element Method applied to (3)–(4) with Euler–Maruyama method (Algorithm 3.1 from [15]) for the stochastic system (2). In columns four and five we write solutions to (3)–(4), i.e. mean extinction-time and its standard deviations. In columns six and seven, we show the mean and standard deviation of solutions to (2) of 500 000 simulations using Algorithm 3.1 from [15] with △t=10−5and nrun =5×105. αR(0)=r F(0)=fττ2−(τ)2Mean Std 0.1 80 11 1.0121 0.7023 1.0139 0.7252 0.1 5 4 3.6906 2.8870 3.6917 2.8816 0.1 5 5 4.0910 3.1151 4.0924 3.1187 0.1 5 6 4.4119 3.2998 4.4109 3.2971 0.05 120 50 2.0121 3.1776 2.0111 3.1780 0.05 8 134 2.0038 3.9895 1.9968 3.9827 2.4. Numerical test We have computed the mean extinction-time, τ, and its standard deviation, τ2−(τ)2, by solving numerically (3)–(4) by the finite element method described in Section 3.1 from [15]. On the other hand, we have also simulated numerically these two quantities by means of the Algorithm 3.1 from [15], in order to make a comparison between both approaches. In Table 1, we show the results corresponding to several tests. In all of them, we have chosen the time-step △t=10−5, and the maximal number of simulations nrun =5·105.Table 1 contains the values of r=R(0)and f=F(0)for different values of αand also the approximations to the mean extinction-time τand the standard deviation τ2−(τ)2(columns four and five, respectively) obtained from the application of the finite element method, and the mean and the standard deviation (columns six and seven, respectively) obtained from Algorithm 3.1 from [15]. Although we have used two completely different approaches to obtain the mean and the standard deviation, their values coincide reasonably well. Furthermore, it is very important to underline the large difference in runtime between the two numerical methods. For instance, the elapsed time for FEM with α=0.1 was about 190 s, while the required time for the results of the first row (which was the fastest one to execute) took more than 8 h and a half. In addition, the proposed method that uses FEM has to be executed just once for each α, while the stochastic simulation of Algorithm 3.1 from [15] has to be run individually for each set of initial data. 3. Extinction-time for a SIS epidemic models In [19] we study two SIS epidemic models, the first one with a constant population size and the second one with deaths. Following [5, p. 148], the stochastic SIS model has the following form: dS =µ(S,I)dt +ρ(S,I)(dW1−dW2), t>0, dI = −µ(S,I)dt +ρ(S,I)(−dW1+dW2), t>0, S(0)=s,I(0)=y, (7) where S=S(t)is the susceptible population size, I=I(t)is the infected population size, sand yare the initial population sizes, W1(t)and W2(t)are two independent Wiener processes, and the drift and diffusion functions are given by µ(S,I)= − αSI S+I+γI, ρ(S,I)=1 2αSI S+I+γI.(8) The parameter α > 0 is the contact rate for passing the infection, i.e. the average number of individuals with whom an infected individual makes sufficient contact to pass an infection, and γ > 0 is the probability that an infected individual is removed from the infection process during a unit time interval (relative removal rate). A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 163 Fig. 3. The numerical solutions of (10) for N=500 in black, N=1000 in blue and N=2000 in red. Left: α=0.25 and γ=0.75 (R<1). Right: α=0.50 and γ=0.49 (R>1). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) We will also assume that the total population is constant, i.e. S(t)+I(t)=Nis constant for t≥0. Therefore, the system (7) can be reduced to a single stochastic differential equation for I:      dI =α−γ−αI NIdt +I(t)α+γ−αI NdW,t>0, I(0)=y. (9) 3.1. Extinction-time for the one dimensional problem The deterministic part of Eq. (9) is two fixed points: I=0 and I=N(R−1)/R, where R=α/γ is called the basic reproductive number. This number is defined as the average number of secondary infections caused by one infective individual during his/her infective period in an entirely susceptible population, and determines the asymptotic behavior with two different behaviors: (i) If R<1,limt→∞ I(t)=0, i.e. the solutions approach the disease-free equilibrium. (ii) If R>1,limt→∞ I(t)=N(R−1)/R, i.e. the solutions approach a unique positive endemic equilibrium. The time to disease extinction T=inf{t≥0:I(t)=0}is the extinction-time for the SDE (9) and its expected value τ= τ(y):= E(T)is the solution to the following boundary value problem for Ordinary Differential Equation (ODE):    1 2yα+γ−αy Nτ′′ +α−γ−αy Nyτ′= −1,y∈(0,N), τ(0)=τ′(N)=0. (10) We have solved numerically this problem in Matlab using Chebyshev differential matrix as it is explained in the Charter 7 of [20]. In Fig. 3 we have plotted the numerical solutions of (10) for two different couple of values of αand γfor N=500, in black, for N=1000 in blue and for N=2000 in red. The left-hand side of Fig. 3 corresponds to α=0.25 and γ=0.75. In this case, R=α/γ =1/3<1. On the right hand side of Fig. 3,α=0.50, γ=0.49, and then R>1. The difference in the type of solutions is evident. While in the first case (left), when R<1, the values of τbasically depend on y=I(0), in the second case (right), when R>1, the solution changes depend on the values of N, conclusion that has already been noted in [21,22]. Moreover, additional tests show that this situation is exacerbated when the basic reproductive number Rincreases. Our computational results show that in this model, there exist a basic reproduction number Rsuch that: If R≤1 The deterministic part and stochastic model approach the disease-free equilibrium, and the time to infection disappearance depends essentially on initial number of infected. If R>1 While the deterministic model approaches a unique positive endemic equilibrium, the stochastic model always tends to a disease-free state and the time to infection disappears depends essentially on the size of the population. 3.2. Extinction-time for the two dimensional problem In order to compute the extinction-time for (7), we consider the stationary backward Kolmogorov equation µ(s,y)∂τ ∂s−∂τ ∂y+ρ2(s,y)∂2τ ∂s2−2∂2τ ∂s∂y+∂2τ ∂y2= −1,(11) 164 A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 Fig. 4. Case α=0.25 and γ=0.75 (R<1). Left: the numerical solution to (11)–(12). Right: the cut at s=899 (blue) for the numerical solution to (11)–(12) is compared with the red points, which represent the means of 10 000 trials of the Euler–Maruyama method applied to (7). Additionally, we have plotted the cuts at s=121 (green) and s=15 (magenta). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) in Ω= [0,MS]×[0,MI], with boundary conditions            τ(0,y)=τ(s,0)=0,s∈ [0,MS],y∈ [0,MI], ∂τ ∂s(MS,y)=0,y∈(0,MI), ∂τ ∂y(s,MI)=0,s∈(0,MS), (12) with MSand MIpositive constants. Let us multiply (11) by a regular test function w=w(s,y)satisfying Dirichlet boundary conditions w(0,y)=w(s,0)= 0 for s∈ [0,MS],y∈ [0,MI], and integrate over the domain Ω. We have Ω µ∂τ ∂s−∂τ ∂swdsdy +Ω 1τ(ρ2w) dsdy −Ω ρ2∂2τ ∂s∂ywdsdy = −Ω wdsdy. Integrating by parts the second order terms and taking into account the Neumann boundary conditions from (12), we obtain the following weak formulation: Ω µ∂τ ∂s−∂τ ∂ywdsdy −Ω∇τ·∇(ρ2w) dsdy +Ω ∂τ ∂y ∂ ∂s(ρ2w) dsdy −MI 0 ∂τ ∂I(MS,y)(ρ2w)(MS,y)dy = −Ω wdsdy,(13) where the last boundary term cannot be removed or replaced by the boundary conditions in (12). Using similar ideas to those in [15], i.e. discretizing by means of rectangular elements with quadratic polynomials as basis functions, and a nodal distribution of the form si=MS(i/NS)3,0≤i≤NSand yj=MI(j/NI)3, 0 ≤j≤NI, particularly good results are obtained. In Fig. 4, we have taken MS=MI=1000 and NS=NI=50, α=0.25 and γ=0.75. From the right-hand side of Fig. 4, it is straightforward to see that the solution depends essentially on y=I(0), because the contour lines are roughly parallel to the horizontal axis. On the left-hand side of Fig. 4 we have represented three cuts of the numerical solution to (11)–(12): at s=889 (blue), s=121 (green) and s=15 (magenta). Red points represent the means of 10 000 trials of Euler–Maruyama method applied to (7) (Algorithm 2.1 in [19] with Tol1=Tol2=1 for S(0)=s=899), and give, in our opinion, an acceptable approximation. The case R>1 is different and more complicated. In Fig. 5 we have taken MS=MI=2000 and NS=NI=100, for α=0.50 and γ=0.49. Now R=1.0204 >1. The elapsed time was about 138.64 s. Even though we have doubled the length of the intervals with respect to Fig. 4, we can hardly see any improvement. This suggests that more nodes are necessary, which requires a higher computational cost. Taking, for instance NS=NI=200, the elapsed time was about 1140.37 s, but hardly any changes can be appreciated in Fig. 5. Observe that, in this case, the situation is very different, because the contour lines on the right-hand side are roughly S(0)+I(0)=N, i.e. as we discussed in the previous subsection: the extinction-time now depends basically on the size of the population N, while in the case R<1, it depends on the initial number of infected individuals, I(0)=y. On the other hand, we can see on the right-hand side of Fig. 5, that the stationary backward Kolmogorov equation (11) works well for small y=I(0), whereas for large y=I(0), the enforcement of homogeneous boundary conditions make the numerical solution diverge from the results of Algorithm 2.1 of [19]. In general, this type of problems will arise when we truncate an unbounded domain to a bounded one. In order to avoid this situation, a promising idea may be to use transparent boundary conditions (TBCs) as it can be seen in [19]. A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 165 Fig. 5. Case α=0.50 and γ=0.49 (R>1). Left: the numerical solution to (11)–(12). Right: the cut at s=1005 (blue) for the numerical solution to (11)–(12) and the blue points represent the means of 10 000 trials of Algorithm 2.1 of [19]. Additionally, we have plotted the cuts at s=500 (red) and s=101 (green). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 2 Zeros of (15) located in each column of the table. Zeros123456 S∗∗∗0∗0 I1∗0 0 ∗∗∗ I20∗∗∗∗∗ To conclude this section, let us notice that we can consider a generalization of the model (7) by introducing a death factors with rate β. This has been proposed in [23]. Obviously, when β=0, we will get (7). In [19] we made numerical approximations of the model with death factors. For the numerical simulations we have chosen FreeFem++, see [24], which is freely available and greatly facilitates the implementation of FEM, and is particularly efficient. 4. Two pathogens epidemic model An epidemic model with two pathogen strains and random demographics can be formulated as follows (see [25,26]):      dS =µ(S,I1,I2)dt +B1,1dW1+B1,2dW2+B1,3dW3, dI1=µ1(S,I1,I2)dt +B2,1dW1+B2,2dW2+B2,3dW3, dI2=µ2(S,I1,I2)dt +B3,1dW1+B3,2dW2+B3,3dW3, S(0)=s,I1(0)=y1,I2(0)=y2, (14) where S=S(t)and I1=I1(t),I2=I2(t)denote continuous random variables susceptible to infection and infected individuals respectively, Wj(t)are three independent Wiener processes and s>0, y1>0 and y2>0 are the initial conditions. The deterministic members in (14) are                µ(S,I1,I2)=Sb−d(N)−β1I1+β2I2 N+b1I1+b2I2, µ1(S,I1,I2)=I1b−b1−d(N)−α1+β1S N, µ2(S,I1,I2)=I2b−b2−d(N)−α2+β2S N, (15) where b>0 is the per capita birth rate, d(N)=1+0.05Nis the per capita death rate, depending of the total density N(t)=S(t)+I1(t)+I2(t). The per capita birth rate to susceptible class Sfrom class Ijis bj≥0, j=1,2. The parameter βj>0, j=1,2 is the transmission rate and αj>0 is the disease-related, per capita death rate for individuals infected with strain Ij,j=1,2. The Symbolic Math Toolbox of Matlab©computes conveniently six zeros of the system (15) which are in Table 2, where ∗denotes nonzero values; we do not write their values due to lack of interest. Each column in this table is a zero of (15). Let us denote the coexistence equilibrium for the deterministic model by E=(S,I1,I2), where S>0,I1>0,I2>0 which corresponds to the fifth column in Table 2. Obviously its asymptotic behavior depends on eigenvalues of the Jacobian. The stochastic part of the model (14) depends on the matrix B=√C, where Cis the following 3 ×3 positive definite matrix: C=σ1,1σ1,2σ1,3 σ2,1σ2,20 σ3,10σ3,3, 166 A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 Fig. 6. No vertical transmission: b=b1=b2=6, β1=30, β2=15, α1=8.5 and α2=3. Deterministic solutions: S(0)=I1(0)=I2(0)=40 (red), S(0)=40, I1(0)=10, I2(0)=30 (magenta) and S(0)=40, I1(0)=1, I2(0)=40 (green) and the mean of 1000 stochastic simulations using the Euler–Maruyama method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) where σ1,1=Sb+d(N)+β1I1+β2I2 N+b1I1+b2I2, σj+1,j+1=Ijb−bj+d(N)+αj+βjS N,j=1,2, σ1,j+1=σj+1,1= −βj SIj N,j=1,2. There exist numerous theoretical studies on the evolution, persistence or extinction of multiple pathogen strains in deterministic epidemic models, see for example [26] and references therein. It has been proved that the dynamics of a deterministic model depends on the basic reproduction numbers defined by Rj=βj+b−bj b+αj ,j=1,2 and it is known that if these numbers are less than one, then the disease extinction occurs. There have been very few investigations (see [4,8]) involving persistence or extinction of multiple pathogen strains for stochastic epidemic models. In [16–18] the simplest simulations for the stochastic model are performed using the Euler– Maruyama method of strong order of convergence 0.5. We have written an algorithm similar to algorithm from [15] or [19] which we do not repeat here. Let us only note that we have three stopping test for the three variables S,I1and I2and for given time step △t, number of simulation nrun and an initial position (S(0), I1(0), I2(0)). The algorithm is straightforward until one of the variable is less than one, because they are all integer. After all the trials, we computed the mean and standard deviation of the stopping times. The Matlab implementation is very similar to the program offered in the appendix of [15]. It is important to remark that in this implementation at each step it is necessary to compute the matrix B=√Cwith Matlab command sqrtm which slows the execution. Let us assume that in (14) no strain is transmitted vertically, i.e. b=b1=b2. We take b=b1=b2=6, β1=30, β2=15, α1=8.5 and α2=3. The basic reproduction numbers are R1=2.0690,R2=5/3. The deterministic solution reaches the equilibrium point E=(11,8.7273,10.2727), that is globally asymptotically stable, because it satisfies Theorem 1 from [26] and the eigenvalues of its Jacobian are negative: −9.4884,−1.5000 and −0.8652. In Fig. 6 we have compared three deterministic trajectories and the mean of 1000 stochastic solutions of (14), which confirms the comments from [26] about the differences in the asymptotic behavior between deterministic and stochastic models. In Table 3 we show the results of 10 000 trials of the Euler–Maruyama method with △t=10−4and for three initial points from Fig. 6. For each initial population size, we have noted the number of stops of each test: S<1 or I1<1 or I2<1 and we have computed the mean stop time and its standard deviation. In view of these results, the population which disappears is I1despite having a higher reproduction number. This strange result suggests that the numerical simulations using Euler–Maruyama method could be spurious. A. Doubova, F. Vadillo / Journal of Computational and Applied Mathematics 295 (2016) 159–169 167 Table 3 No vertical transmission: b=b1=b2=6, β1=30, β2=15, α1=8.5 and α2=3. Results Euler–Maruyama method with △t=10−4. Initial point Number of stops Mean Std S<1 2326 1.2976 0.6818 (40, 40, 40) I1<1 5378 1.1912 0.6660 I2<1 2296 1.3117 0.6565 S<1 1498 1.4583 0.7154 (40, 10, 30) I1<1 7064 0.9038 0.6963 I2<1 1438 1.4943 0.6918 S<1 62 (40, 1, 40) I1<1 9866 0.0080 0.0055 I2<1 72 Fig. 7. Left: case R1,R2>1. Green represents points such that τ≈0.4, blue τ≈0.3 and red τ≈0.2. Right: case R1>1, R2<1. Green represents points such that τ≈0.15, blue τ≈0.1 and red τ≈0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 4.1. The extinction-time Let us denote by Ω= [0,Ms] × [0,M1] × [0,M2], where Ms,M1and M2are positive constants and its boundary by ∂Ω=Γ1∪Γ2∪Γ3∪Γ4, where      Γ1= {(s,y1,y2)∈Ω|s=0,or y1=0,or y2=0}, Γ2= {(s,y1,y2)∈Ω|s=Ms}, Γ3= {(s,y1,y2)∈Ω|y1=M1}, Γ4= {(s,y1,y2)∈Ω|y2=M2}. The mean persistence time τ=τ(s,y1,y2)for (14) satisfies the following partial differential equation: µ(s,y1,y2)∂τ ∂s+ 2  j=1 µj(s,y1,y2)∂τ ∂yj+1 2σ1,1 ∂2τ ∂s2+1 2 2  j=1 σj+1,j+1 ∂2τ ∂y2 j + 2  j=1 σ1,j+1 ∂2τ ∂s∂yj= −1, (s,y1,y2)∈Ω(16) with following boundary conditions:                  τ(0,y1,y2)=u(s,0,y2)=u(s,y1,0)=0 on Γ1, ∂τ ∂s(Ms,y1,y2)=0 on Γ2, ∂τ ∂y1 (s,M1,y2)=0 on Γ3, ∂τ ∂y2 (s,y1,M2)=0 on Γ4. (17) Remember that the variables s,y1and y2represent the initial populations for S,I1and I2respectively. In order to perform numerical experiments using FreeFem++ software, we write (16)–(17) in the variational form, see [27]. The complete FreeFem+++ code is offered in the Appendix of [27]. We consider again the values of parameters as before, i.e. the case where no strain is transmitted vertically: b=b1= b2=6 and β1=30, β2=15, α1=8.5 and α2=3 with the basic reproduction numbers R1=2.0690, R2=5/3. In Fig. 7