Stochastic Predator-Prey model using Gillespie's SSA

Question

I'm trying to model the predator-prey situation using SSA.

I am very new to Mathematica so my handling of matrices is probably terrible, and that's probably adding to my confusion:

ssa[steps_, aa_, bb_, cc_, dd_, R_, F_] :=
 Module[{tr = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}}, time = {0}, 
   data = {{R, F}}, ksteps = steps, a = aa, b = bb, c = cc, d = dd},
  For[k = 1, k <= ksteps, k++,
    time = 
     Append[time, 
      Last[time] + 
       RandomVariate[
        ExponentialDistribution[
         aa*First[Last[data]] + 
          bb*First[Last[data]]*Last[Last[data]] + 
          cc*Last[Last[data]] + 
          dd*First[Last[data]]*Last[Last[data]]]]];
    data = Append[data, Last[data] + tr[[RandomChoice[
          {aa*First[Last[data]], 
            bb*First[Last[data]]*Last[Last[data]], 
            dd*First[Last[data]]*Last[Last[data]], 
            cc*Last[Last[data]]} -> {1, 2, 3, 4}]]]];];
Show[ListPlot[Transpose[First[Transpose[data]]],time],ListPlot[Transpose[Last[Transpose[data]]],time]]

So I start out with a vector t={0} and a vector data={{R,F}}, and i take a randomvariate from the exponential distribution, whose parameter is given by the sum of the individual parameters, add that to my previous time value, and append the new value to my vector t. I then append onto my list of lists one of the four: {1,0},{-1,0},{0,1},{0,-1} each with its own probability weighting.

Finally, I play around with ListPlot until Mathematica stops shouting at me and shows me something.

Now, I don't believe my model works. I have tried the parameters:

ssa[6000, .1, .001, 1, .1, 10, 4]

Sometimes, one species just goes extinct straight away. Sometimes, there is a huge amount of steps elapsed until a species comes back down again.

In any case, I can't seem to get beyond more than 0.5 seconds before Mathematica crashes...

One of the more common errors I get is:

"Parameter 0.` at position 1 in ExponentialDistribution[0.] is expected to be positive."

I don't know where these errors are coming from, and I also don't see where I've gone wrong in modelling the system. Any help would be greatly appreciated. Thanks!

Answer 1

Not acquainted with the SSA algorithm, but I tried to decipher it from your code. Please modify to fix it in your final work. If is not best practice in Mathematica to use For statements, so the following approach is taking advantage of the NestList Command.

    (*Define the function that will provide the vector state after the next event occured*)
ssa[{time_, aa_, bb_, cc_, dd_, r_, f_}] := Module[{t = RandomVariate[
     ExponentialDistribution[aa r +  bb r f + cc f + dd r f]], 
   v = {{1, 0}, {-1, 0}, {0, 1}, {0, -1} }, mc}, 
  mc = RandomChoice[{aa r,  bb r f , cc f , dd r f} -> {1, 2, 3, 
      4}]; {time + t, aa, bb, cc, dd, Max[0, r + v[[mc, 1]]], 
   Max[0, f + v[[mc, 2]]]}]

Simulate and graph one Run of the algorithm

sim = NestList[ssa, {0, .49, .01, .49, .01, 10, 10}, 100];
ListPlot[Evaluate[{sim[[All, {1, 6}]], sim[[All, {1, 7}]]}], 
 Joined -> True, PlotLegends -> {"R", "F"}]