Stochastic Lotka-Volterra Predator-Prey Model

Question

I am struggling with writing a stochastic version of Lotka-Volterra predator-prey model. This is as far as I have gotten:

a = 0.1;
b = 0.2;
c = 0.3;
d = 0.4;
proc = 
  RandomFunction[
    ItoProcess[
      {\[DifferentialD]X[t] == (a \[DifferentialD]t X[t] + \[DifferentialD]W[t]) - 
         b \[DifferentialD]t X[t] Y[t],
       \[DifferentialD]Y[t] == -c \[DifferentialD]t Y[t] + 
       d \[DifferentialD]t Y[t] X[t]}, 
      {X[t], Y[t]}, {{X, Y}, {0.2, 0.2}}, t, 
      W \[Distributed] WienerProcess[0, 0.1]], {0, 100, 0.01}, 
      Method -> "StochasticRungeKuttaScalarNoise"];
ListLinePlot[proc, PlotRange -> All]

For some reason, only one variable depends on the stochastic noise and it does not bear influence on the other part of the equation, I have tried to correct it, but no matter what I do, no luck.

First thing would be to get it working correctly, then I would worry about excluding "extinction Values > 0", since we are not interested in those.

Also, in regards to the graph, is there other way I would create a phase plot (X against Y)? My current way:

ListLinePlot[proc["Values"], PlotRange -> All]

plots X and Y, against time. How would I plot only X in time?

As to the 2nd question, you might try ListLinePlot[Transpose@proc["Values"], PlotRange -> All] — m_goldberg, Jan 20 '18 at 07:12
Oh right, that plots amplitude of both functions (x and y) in time, on the same graph. Would there be a way to separate those functions, and create two graphs, first with only function X values in time, and second with only function Y values in time? @m_goldberg — Peter, Jan 20 '18 at 07:42

score 2 · Answer 1 · answered Jan 21 '18 at 16:15

There are a few questions wrapped up here. Let me try to take them one-by-one.

only one variable depends on the stochastic noise and it does not bear influence on the other part of the equation

I don't think that's actually correct. After running your code,

ListLinePlot[proc, PlotRange -> All]

Mathematica graphics

The predator (gold) seems to smooth over the variation in the prey (blue) but take a look around t=85: there is a random bump in prey which results in a little shoulder in the decline of the predators. So I think the simulation results are accurate for this model and parameter set.

then I would worry about excluding "extinction Values > 0", since we are not interested in those.

That prey can become negative is a problem with the model formulation. I don't think RandomFunction[ItoProcess has the equivalent of WhenEvent in NDSolve.

Two ideas: 1) replace your additive noise \[DifferentialD]W[t] with a multiplicative noise X[t] \[DifferentialD]W[t], which seems like it should prevent the prey population from becoming negative. 2) Instead of using a stochastic differential equation as your model, you could use a more basic continuous time stochastic process and simulate it with @IstvánZachar's GillespieSSA function from this answer (specifically, see his example 2).

is there other way I would create a phase plot (X against Y)?

Funny, when I run your code, it works as desired (v11.2):

ListLinePlot[proc["Values"], PlotRange -> All]

Mathematica graphics

How would I plot only X in time?

This recent answer by @kglr addresses this point, using an undocumented property "PathComponent":

ListLinePlot[proc["PathComponent", 1]]

Mathematica graphics

Anyone want to write a wrapper (call it SNDSolve) for RandomFunction[ItoProcess that mimics NDSolve in input and output? Seems like it'd be appreciated by many!

Thank you very much for the explanation.
In regards to the phase graph, I was just wondering whether mine was a correct way of doing it, since I was not able to find "Values" or "PathComponent" anywhere online...

I will try to tweak it for the extinction event case.

Once again, thank you! — Peter, Jan 22 '18 at 01:48
If this answers your question, you can accept it by clicking on the check mark beside the answer to toggle it from greyed out to filled in. — Chris K, Jan 22 '18 at 14:23

score 1 · Answer 2 · answered Jan 20 '18 at 14:11

This solution will first define the ItoProcess and then simulate it; care is taken to initialize all the available parameters at their appropriate values (using With):

With[{a = 0.1, b = 0.2, c = 0.3, d = 0.4, Xo = 0.2, Yo = 0.2, wm = 0.,ws = 0.1},
 proc = ItoProcess[{
    \[DifferentialD]X[t] == (a - b Y[t]) X[t] \[DifferentialD]t + \[DifferentialD]W[t],
    \[DifferentialD]Y[t] == (-c + d X[t]) Y[t] \[DifferentialD]t
   }, {X[t], Y[t]}, {{X, Y}, {Xo, Yo}},
   t, W \[Distributed] WienerProcess[wm, ws]
  ]
 ];

Next, we simulate 3 paths for the process (proc):

td = RandomFunction[proc, {0, 100, 0.01}, 3, Method -> "StochasticRungeKuttaScalarNoise"]

Finally we plot the derived paths:

Grid[
 Partition[
  MapIndexed[

   With[{val = Part[#1, All, -1]},
     ListLinePlot[
      Transpose[val],
      PlotRange -> All,
      PlotLegends -> {X[t], Y[t]},
      PlotLabel -> Row[{Path, Null, #2[[-1]]}]
    ]] &, td["Paths"]], 2, 2, {1, 1}, Null]]

Stochastic Lotka-Volterra Predator-Prey Model

2 Answers2