Comparing stochastic Lotka–Volterra predator-prey models

Applied Mathematics and Computation 360 (2019) 181–189

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Comparing stochastic Lotka–Volterra predator-prey models Fernando Vadillo Department of Applied Mathematics and Statistics and Operations Research, University of the Basque Country (UPV/EHU), Spain

a r t i c l e

i n f o

MSC: 92B09 60H16 65N19 Keywords: Extinction-time Population models Stochastic Differential Equations Finite Elements Method

a b s t r a c t The phenomenology of population extinction is one of the central themes in population biology which it is an inherently stochastic event. In the present investigation, we study this problem for three different stochastic models built from a single Lotka–Volterra deterministic model. More concretely, we study their mean-extinction time which satisfies the backward Kolmogorov differential equation, a linear second-order partial differential equation with variable coefficients; hence, we can only compute numerical approximations. We suggest a finite element method using FreeFem++. Our analysis and numerical results allow us to conclude that there are important differences between the three models. These differences enable us to choose the most “natural way” to turn a the deterministic model into a stochastic model. © 2019 Elsevier Inc. All rights reserved.

Building models is very different from proclaiming truths. It’s a never-ending process of discovery and refinement, not a war to win or destination to reach. Uncertainty is intrinsic to the process of finding out what you don’t know, not a weakness to avoid. Neil Gershenfeld: Truth Is a Model. In [5].

1. Introduction In 1926 Volterra [42] suggested a simple model for the predation of one species by another to explain the oscillatory of fish catches in the Adriatic. If N(t) represents the prey population and P(t) is the predator population at time t, Volterra’s model is

⎧ dN ⎪ ⎪ = N (a − bP ), ⎨ dt

⎪ ⎪ ⎩ dP = P (cN − d ),

(1)

dt

where a, b, c and d are positive constants; a, d represent the intrinsic growth rates and b, c measure the interaction between two species. This ordinary differential system has already been studied in the classical bibliography, see for example [29, p. 41], [31, p. 79], [38, p. 285], [4, p. 494], [43, p. 76], [1, p. 149], [3, p. 333], [19, p. 366] and [21, p. 2].

E-mail address: [email protected] https://doi.org/10.1016/j.amc.2019.05.002 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

182

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189

In general, a classical two-species deterministic Lotka-Volterra model can be expressed as follows

⎧ dx1 ⎪ ⎪ = x1 [ r1 − a11 x1 − a12 x2 ], ⎨ dt

(2)

⎪ ⎪ ⎩ dx2 = x2 [ r2 − a21 x1 − a22 x2 ], dt

where x1 (t) and x2 (t) represents the sizes of populations. Note that the difference between this system and the Volterra (1) is the new intraspeficic coefficients aii . These systems cover three standard classification Cooperation: r1 > 0, r2 > 0, a12 < 0, a21 < 0; Competition: r1 > 0, r2 > 0, a12 > 0, a21 > 0; Predator-prey: r1 > 0, r2 < 0, a12 > 0, a21 < 0. There exists an extensive literature concerning with the dynamics, see for example [31, Chap. 3], [3, Chap. 7] or [32, Part V] for an general introduction. However, in these deterministic models the parameters are constant and do not incorporate many unpredictable factors that be related to environmental factors and demographic variability (food variability, diseases...) In general, environmental variability is treated by modifying the parameters of the model. This approach consists in defining the parameters with a linear function of Gaussian white noise (WN) i.e. a linear function of the derivative of a Wiener process. For example, in [7,18,23,35,44] or [24] it is assumed that the parameters r1 and r2 are replaced by

rj → rj + σj

dW j (t ) , dt

(3)

where σ j > 0 is positive constant and Wj (t) is a Brownian motion. Therefore we obtain a stochastic model described by the following Itô system

⎧ dX1 dW1 ⎪ ⎪ ⎨ dt = X1 [ r1 − a11 X1 − a12 X2 ] + σ1 X1 dt ,

⎪ ⎪ ⎩ dX2 = X2 [ r2 − a21 X1 − a22 X2 ] + σ2 X2 dW2 , dt

(4)

dt

where the stochastic term is lineal in the stochastic variables X1 and X2 . Under the simple hypothesis σ 1 , σ 2 > 0 the Theorem 2.1 from [27] shows that this solution is positive and global. Another possibility considered in [8,27,28], or [18] where the white noises (WN) appear in the parameters aij such that

a i j → a i j + σi j

dW (t ) , dt

(5)

with the conditions σ ii > 0 and σ ij ≥ 0 if i = j. In this case, the diffusion term in the stochastic differential system would be quadratic. In [25] is studied the difference between these two models for the stochastic non-autonomous logistic equation. The Stochastic Differential Equation (SDE) system for the dynamics of 2 interacting populations has the form

dX = μ(t, X ) + B(t, X ) dW,

(6)

where X = (X1 , X2 )T and W = (W1 , W2 )T are two independent Wiener processes. The vectorial function μ(t, X ) is called the drift, and B(t, X ) the diffusion matrix. A key question in population biology is to understand the constrains that lead to extinction or coexistence of the species. In order to study this question, let define the random variable T that indicates the persistence time i.e.

T ≡ inf{t ≥ 0 : X1 (t ) = 0 or X2 (t ) = 0}, obviously T depends on the initial value X1 (0), X2 (0) although it is not explicitly indicated. As discussed in [1, p. 150], the mean persistence-time τ (r, f ) ≡ E(Tr, f ) for (6) satisfies the stationary backward Kolmogorov equation

L ( τ ) ≡ μ1

2 2 ∂τ ∂τ 1   ∂ 2τ + μ2 + di j = −1, ∂ x1 ∂ x2 2 i=1 j=1 ∂ xi x j

where the diffusion term is the matrix D = B · BT = (di j ) and with boundary conditions

(7)

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189

⎧ τ ( 0, x2 ) = τ ( x1 , 0 ) = 0, ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ∂τ ( M , x ) = 0, ∂r 1 2 ⎪ ⎪ ⎪ ⎪ ∂τ ⎪ ⎩ ( x , M ) = 0, ∂f 1 2

183

(8)

provided that the number of x1 and x2 cannot exceed some values M1 and M2 respectively. Moreover, the operator L is elliptic because the matrix D is positive definite in the domain, in fact

x T · D · x = ( x T B ) · ( B T x ) = ( BT x ) T · ( B T x ) =  BT x

22 ≥ 0,

therefore its eigenvalues are positive (see for example [11]). To solve numerically the boundary problems (7) and (8), we will use the Finite Element Method (FEM). The FEM has been invented by engineers around 1950 for solving the partial differential equations arising in solid mechanics; the main idea was to use the principe of virtual work for designing approximations of the boundary value problems. The more popular reference for FEM applied to solid mechanics is the books [45] and [46]. Generalization to others fields of physics or engineering have been done by applied mathematicians through the concept of variational formulation of partial differential equations (see for example [6]). There is an ample bibliography about finite element method and their implementation (see for instance [12] for a good introduction). As for this study, we computed FEM using FreeFem++, see [13], which is freely available. This paper is organized as follows. In Section 2, we describe the three stochastic Lotka-Volterra predator-prey models and obtain the backward Kolmogorov equation. In Section 3 we describe the computational implementation of the equation. Finally in Section 4 we analyze the numerical results and draw the main conclusions. Our numerical methods were implemented in Matlab and FreeFem++, the codes are available on request. The experiments were carried out in an Intel(R) Core(TM)2 Duo CPU U9300 @ 1.18GHz, 1.91 GB of RAM. 2. Three stochastic models In this paper, we consider a very simple Lotka-Volterra predator-prey system (1) with a = 2, d = 1 and b = c = α . Let suppose the preys are rabbits that enjoy an infinite supply of food and the predator ares foxes. This problem can be modeled by a nonlinear ordinary differential equations

⎧ dr ⎪ ⎪ ⎨ dt = r (2 − α f ),

⎪ ⎪ ⎩ df = f (−1 + α r ),

(9)

dt

where r(t) and f(t) represent the rabbit and foxes population, and the parameter α ≥ 0. If α = 0, the two populations do not interact and the foxes die from starvation, and when α > 0, foxes meet rabbits with a probability that is proportional to the product of their numbers. Although there is no analytical solution to this system, we can draw a phase plane using matlab software (see for example [15,37] or [30]). In Fig. 1 from [17] we can observe a center in (1/α , 2/α ) and that the solutions are located in closed and periodic orbits. Let now introduce demographic and environmental variability as discussed in the introduction. Model 1: Demographic variability, sometimes referred as Lagevin approach, was first introduced by Langevin in [22]. The idea is to approach the variation of the process. For further information, please refer to [3, p. 359] or [1, p. 149]. Suppose t a small time interval, the corresponding

R = R(t + t ) − R(t ), F = F (t + t ) − F (t ), there are five possibilities for these changes

Prob{R = 1, F = 0} = 2R(t )t ,

Prob{R = 0, F = 1} = α R(t )F (t )t,

Prob{R = −1, F = 0} = α R(t )F (t )t, Prob{R = 0, F = −1} = F (t )t , Prob{R = 0, F = 0} = 1 − P (t )t , where P (t ) = 2R(t ) + F (t ) + 2α R(t )F (t ). To obtain the stochastic differential equations we compute the expectations and covariance of the stochastic variables R and F (see for example [41]). In this case, the stochastic differential

184

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189 200 180 160 140

f

120 100 80 60 40 20 0 0

20

40

60

80

100

120

140

160

r Fig. 1. Three numerical solutions beginning in R(0 ) = 150, F (0 ) = 50: the red line is the mean of 10 0 0 trials for Model 1 in 0 ≤ t ≤ 2.5. The blue line is the Model 2 with σ = 0.001 and 0 ≤ t ≤ 12 and finally the green line is the of Model 3 with σ1 = σ2 = 0.001 in 0 ≤ t ≤ 12. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

system is





dR = R (2 − α F ) dt + R(2 + α F ) dW1 (t ), dF = F (α R − 1 ) dt + F (α R + 1 ) dW2 (t ),

(10)

where dW1 (t) and dW2 (t) are two independent Wiener processes. This model has the single parameter α . Model 2: As the environmental variability is introduced by modifying the parameters of the model, let consider the following modification of the parameter α

α → α+σ

dW (t ) , dt

where W(t) is a Brownian motion and σ > 0 is a constant. Therefore, if R(t) and F(t) represent the size of the population of rabbits and foxes respectively, they satisfying the stochastic differential system



dR = R (2 − α F ) dt − σ R F dW (t ), dF = F (α R − 1 ) dt + σ R F dW (t ),

(11)

here the term σ R F represents the weighted contact term with noise intensity rate σ . This type of system has been studied in the literature: [28], [8] or [25], however, it has always been assumed that the coefficients ajj > 0 in (2), an assumption that is not reasonable in our model. Model 3: Finally, we consider a third model where we introduce distinct linear perturbations to our variations. We get the following stochastic system



dR = R (2 − α F ) dt + σ1 R dW1 (t ), dF = F (α R − 1 ) dt + σ2 F dW2 (t ),

(12)

where dW1 (t) and dW2 (t) are two independent Wiener process, the terms σ 1 R and σ 2 F represent the noise intensity rates for rabbits and foxes respectively. This model is a particular case of 4 and therefore its solutions are global and positive when σ 1 , σ 2 > 0. This model is also studied in [9,35,36,44, Chap. 8] or [24]. Our numerical simulations for the stochastic Model 1,2 and 3 have been done using the Euler-Maruyama, this method is particulary simple and straightforward to implement (see for example [20] or [14]). Fig. 1 shows three possible solutions with the initial conditions R(0 ) = 150, F (0 ) = 50. For model 1, the red line represents the mean of 10 0 0 trials for 0 ≤ t ≤ 2.5. For t > 2.5 we get negative which imply complex square roots. For model 2, the mean on 10 0 0 trials is represented with a blue line. Recall the solution to model 2 consists in numerical approximations. We set σ = 0.001 and 0 ≤ t ≤ 12, finally the green line is the numerical solution of Model 3 with σ1 = σ2 = 0.001 in 0 ≤ t ≤ 12. This Figure suggests for short periods of time there are small differences between the three models. 2.1. Backward Kolmogorov equations Once the three stochastic models have been calculated, we study their backward Kolmogorov Eq. (7).

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189

185

In Model 1, the diffusion term is defined by the matrix



D1 =



1 r (2 + α f ) 0 2

0 , f (α r + 1 )

(13)

and the mean extinction-time τ (r, f ) ≡ E(Tr, f ) satisfies the stationary backward Kolmogorov equation

L1 ( τ ) = r ( 2 − α f )

∂τ ∂τ 1 ∂ 2τ 1 ∂ 2τ + f (α r − 1 ) + r (2 + α f ) 2 + f (α r + 1 ) 2 = −1. ∂r ∂f 2 2 ∂r ∂f

(14)

with some boundary conditions (8). Moreover, as r (2 + α f ) > 0 and f (α r + 1 ) > 0 when r, f > 0 the operator L1 is locally uniformly elliptic and the solution of the problem (14)–(8) is nonnegative (Theorem 1 in [17]). For Model 2 the diffusion coefficient matrix is



D2 =

1 2 2 2 1 σ r f −1 2



−1 , 1

(15)

and the new stationary backward Kolmogorov equation is



∂τ ∂τ 1 2 2 2 ∂ 2 τ ∂ 2τ ∂ 2τ L2 ( τ ) = r ( 2 − α f ) + f (α r − 1 ) + σ r f −2 + = −1. ∂r ∂f 2 ∂ r∂ f ∂ f 2 ∂ r2

(16)

Now, the eigenvalues in D2 are {0, 1} and the L2 operator is not uniformly elliptic. Consequently we can not prove that solutions of (16)–(8) are always positive. Moreover, it cannot be applied the Feyman-Kac formula in [0, M1 ] × [0, M2 ] (see for example [26, pp. 78]). However, we can apply this formula in [ε , M1 ] × [ε , M2 ] although it is necessary a bit change in the definition (see for example [24]). It is said that a population ξ t dies out if: (a) ξt (ω ) = 0 for t ≥ t0 (ω). (b) limt→∞ P rob(ξt ≥ ε ) = 0 for each ε > 0. It clear (a)⇒(b). Then, we apply this weaker definition in [ε , M1 ] × [ε , M2 ], in fact, in our numerical experiments ε = 0.1 and it seems reasonable to predict that the times will be much longer. Finally for Model 3 the diffusion term matrix is

D3 =

1 2



σ12 r2 0

0



σ22 f 2

,

(17)

and the new stationary backward Kolmogorov equation

L3 ( τ ) = r ( 2 − α f )

∂τ ∂τ 1 2 2 ∂ 2 τ 1 2 2 ∂ 2 τ + f (α r − 1 ) + σ r + σ f = −1. ∂r ∂ f 2 1 ∂ r2 2 2 ∂ f 2

(18)

In this model there exists global and positive solutions for σ 1 , σ 2 > 0 and we compute in [ε , M1 ] × [ε , M2 ]. 3. Numerical results In order to apply the FEM technique, we would need to write a well-suited variational form for (14)–(8), (16)–(8) and (18)–(8). This is a standard technique hence we will skip the details. In order to apply this method, we multiply the equation by a test function and set the boundary constraints. We introduce this new equation in FreeFem++ and run the software. On the left side of Fig. 2, we have represented the numerical solution to Model 1 with α = 0.05, This result obtained with FreeFem++ is similar to the solution to equation in [17] obtained with MATLAB. On the right side of Fig. 2, we plotted the numerical solution to Model 2 with α = 0.05, σ = 0.01. Fig. 3 also represents the numerical solution to Model 2 with different parameters α = 0.05, σ = 0.005 and σ = 0.001 respectively. In the second model, when σ decreases, the diffusion coefficient in Eq. (16) also decreases and then layer boundary problem appear. This leads to the layer boundary problem that requires special methods (see for example [40] or [10] and its references). In particular, the method suggested by [34] may be an excellent tool to overcome these difficulties. In fact, we see that the peak grow quickly and so we can conclude that the second model is extinguished much later. In order to confirm these numerical results, we compare the mean extinction-time for the models (11) and (10) using the Euler-Maruyama method. For each simulation we move forward in time until we verify the conditions Rn ≤ eps or Fn ≤ eps ≈ 2.2204 · 10−16 , the step time t = 10−5 and in each case 10 0 0 trials have been performed. In Table 1 shows the mean and the standard deviation of 10 0 0 simulations. We can appreciate the significant differences between two models as we already see in Figs. 2 and 3. Finally, the long-time behaviour of a stochastic prey-predator (12) (Model 3) can be analyzed using [35]. Let define the coefficients

c1 = 2 −

1 2 σ , 2 1

c2 = 1 +

1 2 σ > 0, 2 2

the Theorem 1 from [35] predicts two different behaviors:

186

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189

Fig. 2. On the left hand side the numerical solution to (14)–(8) with α = 0.05. The right side is the numerical solution to (16)–(8) with α = 0.05, σ = 0.01.

Fig. 3. Numerical solution to (16)–(8) with α = 0.05, σ = 0.005 on the left and σ = 0.001 on the right. Table 1 Results of 10 0 0 simulations with α = 0.05. Model

R(0)

F(0)

σ

Mean

Std

1 2 2 2

120 120 120 120

50 50 50 50

0.01 0.005 0.001

2.0111 40.0851 136.7820 2.6023 × 103

3.1780 16.9583 60.1283 2.1075 × 103

(i) If c1 < 0 ⇔ σ 1 > 2. Both prey and predator population die. In fact, the prey population dies even if there are no predators. (ii) If c1 > 0 ⇔ σ 1 < 2. Both species survive. In other words, this result means that only when the stochastic noise of the prey R is sufficiently large (σ 1 > 2) both species die. In addition, we can note that the stochastic noise in the predator F it is hardly relevant. In order to test this observation, we have plotted the numerical solution to (18)–(8) with α = 0.05, σ1 = 2.1, σ2 = 0 and σ2 = 2, Fig. 4 confirms this observation as in a short time period one of the two populations extinguishes, similarly Model 1, one of the two populations is zero; further considering these plotters, we believe that it is the prey whose hit the zero value. Moreover, Fig. 5 shows the numerical solutions to Model 3 where the σ 1 and σ 2 are very small, in this case, theoretically both species remain. These numerical results confirm with Theorem 1 in [44] and Theorem 3 in [24], nevertheless, it should be noted that the theorems assumed not null intraspeficic coefficients i.e. a11 = 0, a22 = 0 in system (2) which is not true in our example. Therefore, our results suggest a possible extension of these theorems to models with null intraspeficic coefficients.

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189

187

Fig. 4. Numerical solution to (18)–(8) with α = 0.05, σ1 = 2.1 and σ2 = 0 and σ2 = 2 respectively.

Fig. 5. Numerical solution to (18)–(8) with α = 0.05 and several σ 1 , σ 2 .

Fig. 6. Numerical solution of backward Kolmogorov equation in the Model 1 applied to (1) with a = 10, d = 5 and b = c = 0.05, 0.025 on left and right respectively.

188

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189 Table 2 Comparing estimates of τ . Method

N(0)

P(0)

a

b=c

d

τ

[39] Euler-Maruyama Backward Kolmogorov [39] Euler-Maruyama Backward Kolmogorov

100 100 100 200 200 200

200 200 200 400 400 400

10 10 10 10 10 10

0.05 0.05 0.05 0.025 0.025 0.025

5 5 5 5 5 5

12 12.60 12.62 24 25.09 25.58

4. Conclusions In this paper, we have considered three stochastic Lotka-Volterra predator-prey. Starting with a classical Volterra predator-prey model we obtained the differential (9). The Model 1 following [3, p. 359] or [1, p. 149], the Models 2 introduce white noise through the parameter α and Model 3 adds white noise to the ordinary differential system (9). Although the models look similar at first sight (Fig. 1). We have computed the mean extinction-time solving numerically the backward Kolmogorov equations for each model with completely different results: This can lead to interesting biological interpretation. The first contribution of the present investigation arises from the conclude that there are important differences in their behaviors regarding their extinctions or not. 1. The Model 1, both population extinguish in a reasonable time, obviously depending of the parameter α and the initial values of the populations. 2. In the Model 2, the mean extinction-time depends on new parameter σ , the white noise. However, for α = 0.05 and a reasonable σ ≈ 0.001 this mean large extinction-time. 3. Finally for the Model 3, we observed great differences depending whether σ 1 < 2 or not. A reasonable situation would think that 0 < σ 1 , σ 2 0.1 and this mean very large extinction-time again. At this point, the question that arises is which is the most pertinent model that suggests a reasonable extension of the initial deterministic model (9) or more general (1). To answer this question, we compare our results with some estimations in [33,39]. Please note that their results are obtained studying the evolution of conservation quantities in (1), more specifically, they suggest an algorithm to compute the probability density function of extinction time. In Fig. 6 we plot the numerical solution of backward Kolmogorov equation in the Model 1 applied to (1) with a = 10, d = 5 and b = c = 0.05 on the left and, b = c = 0.025 on right, they are examples (a) and (e) in Table 1 from [39]. These figures suggest that their estimations are very close to our numerical results. Table 2 summarises the results, first the result from [39], second mean of 20 0 0 trials using Euler-Maruyama method and finally the value of the backward Kolmogorov equation. Although we have used three completely different approaches to obtain τ , their values do coincide reasonably well. At this point and based on the results my recommendation woulds be to use model 1. Indeed, models 2 and 3 present limited validity that should be studied in future investigations. On the order hand, recently in [2] has been defended the mean-reverting processes, however, in my opinion with unclear arguments. Finally, in this paper we have demonstrated how we made in [16,17], that the numerical solution for the backward Kolmogorov equation can provide valuable information about extinction forecasts. Acknowledgments This work was supported by the Spanish Ministry of Economy and Competitiveness with the project MTM2014-53145 and by the Basque Government with the project IT-IT-641-13. The author would like to thank the referees for the thorough, constructive and helpful comments and suggestions on the manuscript. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

E. Allen, Modeling with It Stochastic Differential Equations, Springer, 2007. E. Allen, Enviromental variability and mean-reverting processes, Discrete Cont. Dyn. Syst. Series B 21 (7) (2016) 2073–2089. L.J.S. Allen, An Introduction to Stochastic Processes with Applications to Biology, Person Prentice Hall, 2003. W. Boyce, R. Diprima, Elementary Differential Equations and Boundary Value Problems, 4rd, John Wiley and Sons, 1986. J. Brockman, This Will Make You Smarter. New Scientific Concepts to Improve Your Thinking, Harper Perennial, 2016. P.G. Ciarlet, Basic error estimates for elliptic problems, Hand Book of Numerical Analysis, Vol. II, North-Holland, 1991, pp. 17–351. C. Dong, L. Liu, Y. Sun, Partial permanence and extinction on stochastic Lotka-Volterra competitive systems, Adv. Differ. Equ. (2015) 1–17. N.H. Du, V. Sam, Dynamics of a stochastic Lotka-Volterra model perturbed by white noise, J. Math. Anal. Appl. 324 (1) (2006) 82–97. R. Durretr, Stochastic Calculus: A Practical Introduction, CRC Press, 1996. A. Fortin, J. Urquiza, R. Bois, A mesh adaptation method for 1D-boundary layer problems, Int. J. Numer. Anal. Model. Series B 3 (4) (2012) 408–428. D. Gilbarg, N. Trudinger, Elliptic Partial Differential Equations of Second Order, 2nd, Springer, 2001. M. Gockenbach, Understanding and Implementing the Finite Element Method, SIAM, 2006. F. Hecht, New development in freefem++, J. Numer. Math. 20 (3–4) (2012) 251–265. D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (3) (2001) 525–546.

F. Vadillo / Applied Mathematics and Computation 360 (2019) 181–189

189

[15] D.J. Higham, N.J. Higham, MATLAB Guide, SIAM, 20 0 0. [16] F. de la Hoz, A. Doubova, F. Vadillo, Persistence-time estimation for some stochastic SIS epidemic models, Discrete Count. Dyn. Syst. Series B 20 (8) (2015) 2933–2947. [17] F. de la Hoz, F. Vadillo, A mean extinction-time estimate for a stochastic Lotka-Volterra predator-prey model, Appl. Math. Comput. 219 (1) (2012) 170–179. [18] D. Jiang, C. Ji, X. Li, D. O’Regan, Analysis of autonomous Lotka-Volterra competition systems with random perturbation, J. Math. Anal. Appl. 390 (2) (2012) 582–595. [19] F.C. Klebaner, Introduction to Stochastic Calculus with Applications, Imperial College Press, 2005. [20] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Cambridge University Press, 1998. [21] H. Kokko, Modelling for Field Biologists and Other Interesting People, Cambridge University Press, 2007. [22] P. Langevin, Sur la théori du mouvement brownien, Comptes-rendus de l’Académie des Sciences (146) (1908) 530–533. [23] X. Li, X. Mao, Population dynamical behavior of non-autonomous Lotka-Volterra competitive system with random perturbation, Discrete Count. Dyn. Syst. 24 (2009) 523–545. [24] M. Liu, M. Fan, Permanence of stochastic Lotka-Volterra systems, J. Nonlinear Sci. 27 (2017) 425–452. [25] M. Liu, K. Wanga, Persistence and extinction in stochastic non-autonomous logistic systems, J. Math. Anal. Appl. 375 (2) (2011) 443–457. [26] X. Mao, Stochastic Differential Equations and Applications. Second Edition, Woodhead Publishing, 2001. [27] X. Mao, G. Marion, E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stoch. Proceses Appl. 977 (1) (2002) 95–110. [28] X. Mao, S. Sabanis, E. Renshaw, Asymptotic behaviour of the stochastic Lotka-Volterra model, J. Math. Anal. Appl. 287 (1) (2003) 141–156. [29] R. May, Stability and Complexity in Model Ecosystems, Cambridge University Press, 1973. [30] C. Moler, Numerical Computing with MATLAB, SIAM, 2004. [31] J. Murray, Mathematical Biology I: An Introduction, 3rd, Springer-Verlag, 2002. [32] D. Neal, Introduction to Population Biology, Cambridge University Press, 2004. [33] M. Parker, A. Kamenev, Mean extinction time in predator-prey model, J. Stat. Phys. (141) (2010) 201–216. [34] B. Reutera, V. Aizingera, M. Wielanda, F. Frankb, P. Knabnera, FESTUNG: A MATLAB /GNU Octave toolbox for the discontinuous Galerkin method (2016). [35] R. Rudnicki, Long-time behaviour of a stochastic prey-predator model, Stoch. Processes Appl. 108 (2003) 93–107. [36] R. Rudnicki, K. Pichór, Influence of stochastic perturbation on prey-predator systems, Math. Biosci. 206 (2007) 108–119. [37] L. Shampine, I. Gladwell, S. Thompson, Solving ODEs with MATLAB, Cambridge University Press, 2003. [38] G.F. Simmons, Differential Equations with Applications and Historical Notes, McGraw-Hill International Editions, 1972. [39] A. Skvortsov, B. Ristic, A. Kamenev, Predicting population extinction from early observations of the Lotka-Volterra system, Appl. Math. Comput. 320 (2018) 371–379. [40] B.G. Smith, L. Vaughan-Jr., D.L. Chopp, The extended finete element method for boundary layer problems in biofilm growth, Comm. App. Math. and Comp. Sci. 2 (1) (2007) 35–56. [41] F. Vadillo, ¿ Modelos deterministas o estocásticos ? Matematicalia 7 (4) (2011) 1–7. [42] V. Volterra, Variazionie fluttuazioni del numero de individui in specie animali conviventi, Mem. Acad. Lincei. 2 (1926) 31–113. [43] G. de Vries, T. Hillen, M. Lewis, J. Müller, B. Schönfisch, A Course in Mathematical Biology, SIAM, 2006. [44] Y. Zhu, M. Liu, Permanence and extinction in a stochastic service-resource mutualism model, Appl. Math. Lett. 69 (2017). [45] O.C. Zienkiewick, R.L. Taylor, The finite element method, Basis Formulation and Linear Problems, 1, McGraw-Hill, 1989. [46] O.C. Zienkiewick, R.L. Taylor, The finite element method, Solid and Fluid Mechanics: Dynamics and Non-linearity, 2, McGraw-Hill, 1991.