Using RegionPlot and NDSolve to plot areas of species coexistence across 2 parameter space

Question

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!

Answer 1

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]