I have a predator prey food chain model set up as follows, with x representing the resource, y the consumer, and z the predator:

dxdt[x_, y_, z_] := r x (1 - x/k) - a12 x y - a13 x z;
dydt[x_, y_, z_] := e12 a12 x y - a23 y z - m2 y;
dzdt[x_, y_, z_] := e23 a23 y z + e13 a13 x z - m3 z;

model[{x0_, y0_, z0_}] := {
x'[t] == dxdt[x[t], y[t], z[t]],
y'[t] == dydt[x[t], y[t], z[t]],
z'[t] == dzdt[x[t], y[t], z[t]],
x[0] == x0,
y[0] == y0,
z[0] == z0

I have the following code to plot the equlibrium biomasses of the resource predator and consumer over time and view how they change depending on the values of the parameters in the model using the manipulate function:

 Module[{ sol, tend = 1000},
  sol =
    model[{1.0, 0.5, 0.1}] /. {a12 -> a12in, a23 -> a23in, 
       a13 -> a13in, m2 -> min, m3 -> min} /. param, {x, y, z}, {t, 0,
      tend}, perturb];
  Plot[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tend}, 
   PlotLegends -> {"R", "C", "P"}, Frame -> True, 
   PlotRange -> {All, {0, 2}}]
 {{a12in, 1.0}, 0.01, 1.5, Appearance -> "Labeled"},
 {{a23in, 1.0}, 0.01, 1.5, Appearance -> "Labeled"},
 {{a13in, 0.0}, 0.0, 1.0, Appearance -> "Labeled"},
 {{rin, 3.0}, 1.0, 3.0, Appearance -> "Labeled"},
 {{e12in, 0.8}, 0.0, 3.0, Appearance -> "Labeled"},
 {{e23in, 0.8}, 0.0, 3.0, Appearance -> "Labeled"},
 {{min, 0.7}, 0.0, 2.0, Appearance -> "Labeled"}

I want to add a perturbation to the model, at a specified time interval after the biomasses reach equilibrium. I want the perturbation to be something like, z[t] -> 0.9z[t], where 10%, or another percentage, of the predator population is lopped off. I then want to be able to measure and record the time and distance required for the predator population to return to equilibrium, as well as any population overshoot that occurs, in order to calculate stability metrics for the system.

I'm aware that I could likely use something like a WhenEvent function.

I would also like to be able to add a periodic increase in resource density at a specific recurring time interval. Something that would allow me to increase resource density by x% every 50th or so time step.

  • $\begingroup$ Welcome to M.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Chris K
    Nov 28, 2017 at 23:46
  • $\begingroup$ I couldn't get your code as posted to run. Could you have a look at it? Rather than Manipulate, can you give a set of default parameter values? Also, when you say you want to increase the resource density, do you mean x[t] or its carrying capacity k? $\endgroup$
    – Chris K
    Nov 28, 2017 at 23:49


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy