1
$\begingroup$

I have four ODEs corresponding to four species:

dRdt = g*R (1 - R/k) - c1*M*R - c2*P*R;
dMdt = e1*c1*M*R - c3*Z*M*(1 - s) - M*d1;
dPdt = e2*c2*P*R - c5*Z*P*s - u1*P*(M + P) + u2*R*Z - P*d2;
dZdt = e3*c3*Z*M*(1 - s) + e5*c5*Z*P*s + u1*P*(M + P) - u2*R*Z - Z*d3;

I check that everything works by first simulating the dynamics over time:

sim = {R -> R[t], M -> M[t], P -> P[t], Z -> Z[t]};

dR = {R'[t] == First@{dRdt /. sim}, R[0] == 1};
dM = {M'[t] == First@{dMdt /. sim}, M[0] == 1};
dP = {P'[t] == First@{dPdt /. sim}, P[0] == 1};
dZ = {Z'[t] == First@{dZdt /. sim}, Z[0] == 1};

Assigning my parameters and solving with NDSolve:

pars = {g -> 5, k -> 25, c1 -> 1, c2 -> 0.5, c3 -> 0.5, 
   c5 -> 0.5, e1 -> 0.5, e2 -> 0.5, e3 -> 0.5, e5 -> 0.5, d1 -> 0.5, 
   d2 -> 0.5, d3 -> 0.5, u1 -> 0.5, u2 -> 0.5, s -> 0.7};
sol = NDSolve[{{dR, dM, dP, dZ} /. pars}, {R, M, P, Z}, {t, 0, 
    1000}];

And finally plotting the simulations across my two variables of interest:

Manipulate[
 sol = NDSolve[{{dR, dM, dP, dZ} /. {g -> 5, k -> var1, c1 -> 0.5, 
      c2 -> 0.5, c3 -> 0.5, c5 -> 0.5, e1 -> 0.5, e2 -> 0.5, 
      e3 -> 0.5, e5 -> 0.5, d1 -> 0.5, d2 -> 0.5, d3 -> 0.5, 
      u1 -> 0.5, u2 -> 0.5, s -> var2}}, {R, M, P, Z}, {t, 0, 
    10000}];
 
 Plot[Evaluate[{R[t], M[t], P[t], Z[t]} /.
    First[sol]], {t, 10, 
   10000}, PlotStyle -> cols, PlotLegends -> Automatic, 
  PlotRange -> {{0, 1000}, {0, 10}}],
 
 {{var1, 50}, 0, 100, 1},
 {{var2, 0.7}, 0, 1, 0.01}
 ]

This is great because I can use the toggle to visualize when certain species go extinct vs. when they all coexist. Now, I want to show exactly that using RegionPlot. Where I'm stuck is that I'm not sure how to ask RegionPlot to evaluate the areas of all species having equilibrium greater than zero as being True. I've tried the following but I know this is wrong:

Table[
 sol = NDSolve[{{dR, dM, dP, dZ} /. {g -> 5, k -> var1, c1 -> 0.5, 
      c2 -> 0.5, c3 -> 0.5, c5 -> 0.5, e1 -> 0.5, e2 -> 0.5, 
      e3 -> 0.5, e5 -> 0.5, d1 -> 0.5, d2 -> 0.5, d3 -> 0.5, 
      u1 -> 0.5, u2 -> 0.5, s -> var2}}, {R, M, P, Z}, {t, 0, 
    10000}];
 
 RegionPlot[
  Evaluate[{R[t], M[t], P[t], Z[t]} /.
     sol] > 0, {0, 100, 1}, {0,
    1, 0.01}],
 
 {{var1, 50}, 0, 100, 1},
 {{var2, 0.7}, 0, 1, 0.01}
 ]

Any help is greatly appreciated, thank you!

$\endgroup$

1 Answer 1

3
$\begingroup$

We use NDSolve to simulate from t=0 to t=10000, and the value at time t=10000 is used as a proxy for the equilibrium value. For the plot I use ListContourPlot:

(* OPs code *)
dRdt=g*R (1-R/k)-c1*M*R-c2*P*R;
dMdt=e1*c1*M*R-c3*Z*M*(1-s)-M*d1;
dPdt=e2*c2*P*R-c5*Z*P*s-u1*P*(M+P)+u2*R*Z-P*d2;
dZdt=e3*c3*Z*M*(1-s)+e5*c5*Z*P*s+u1*P*(M+P)-u2*R*Z-Z*d3;
sim={R->R[t],M->M[t],P->P[t],Z->Z[t]};
dR={R'[t]==First@{dRdt/. sim},R[0]==1};
dM={M'[t]==First@{dMdt/. sim},M[0]==1};
dP={P'[t]==First@{dPdt/. sim},P[0]==1};
dZ={Z'[t]==First@{dZdt/. sim},Z[0]==1};

(* adaptation of OPs code *)
aux={g->5,c1->0.5,c2->0.5,c3->0.5,c5->0.5,e1->0.5,e2->0.5,
     e3->0.5,e5->0.5,d1->0.5,d2->0.5,d3->0.5,u1->0.5,u2->0.5};
equilibrium[var1_,var2_,tend_:10000]:=NDSolveValue[{{dR,dM,dP,dZ}/.aux/. {k->var1,s->var2}},             
                                         {R[t],M[t],P[t],Z[t]},{t,0,tend}]/.{t->tend};
                                         
(* simulate, this takes a moment *)
var1s=Range[1,100,0.5];
var2s=Range[0.01,1,0.005];
sim=Outer[{#1,#2,equilibrium[#1,#2]}&,var1s,var2s];

(*plot*)
plotdata[i_]:=Flatten[Map[{#[[1]],#[[2]],#[[3,i]]}&,sim,{2}],1];
plotdata["all"]:=Flatten[Map[{#[[1]],#[[2]],Times@@Chop[Abs[#[[3]]]]^(1/4)}&,sim,{2}],1];
plot[data_,species_,color_]:=ListContourPlot[data,Contours->{0.0001},PlotLabel->"Species "<>species,ContourShading->{None,Opacity[0.2,color]},InterpolationOrder->3,FrameLabel->{"var1","var2"}];
GraphicsGrid[{{plot[plotdata[1],"R",Blue],plot[plotdata[2],"M",Green]},
             {plot[plotdata[3],"P",Yellow],plot[plotdata[4],"Z",Magenta]}}]
plot[plotdata["all"],"ALL",Cyan]

enter image description here

enter image description here

$\endgroup$
3
  • $\begingroup$ Very neat, thank you! Is there any way to overlay these plots onto one another so that there is one single plot with a shaded region where all species have positive equilibrium? $\endgroup$
    – Clara
    Jul 29 at 16:28
  • $\begingroup$ @Clara I added another plot. I understand that the contours do not look great, one could work on it. Btw, I hacked this together pretty quickly, and have not checked carefully if the result is plausible, that will be much easier for you to judge. Good luck. $\endgroup$
    – user293787
    Jul 29 at 16:48
  • $\begingroup$ That is perfect, thank you so much for your help! $\endgroup$
    – Clara
    Jul 29 at 19:14

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.