Stochastic Lotka-Volterra Predator-Prey Model

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 7: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 7:42

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] 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] 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]] 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!

• Or maybe NSDSolve makes more sense ;) – Chris K Jan 21 '18 at 16:40
• 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 1: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

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]] 