Provision of essential resources as a persistence strategy in food webs

Michael Raatz (mraatz@evolbio.mpg.de)

Notebook #1 - calculating all but Fig. 3

ClearGlobal[]:=(ClearAll["Global`*"];Clear[Derivative];);RemoveGlobal[]:=(ClearGlobal[];Remove["Global`*"];);RemoveGlobal[]

SetOptions[Plot,BaseStyle->{FontFamily->"Helvetica",FontSize->14,FontWeight->Plain,FontColor->Black}];
SetOptions[Plot,{FrameTicksStyle->None,FrameTicks->Automatic,FrameStyle->Thickness[Medium],LabelStyle->Black, PlotStyle->Thickness[0.015],PlotPoints->40}];
SetOptions[LogPlot,BaseStyle->{FontFamily->"Helvetica",FontSize->14,FontWeight->Plain,FontColor->Black}];
SetOptions[LogPlot,{FrameTicksStyle->None,FrameTicks->Automatic,FrameStyle->Thickness[Medium],LabelStyle->Black}];
SetOptions[RegionPlot,BaseStyle->{FontFamily->"Helvetica",FontSize->14,FontWeight->Plain,FontColor->Black}];
SetOptions[RegionPlot,{FrameTicksStyle->None,FrameTicks->Automatic,FrameStyle->Thickness[Medium],LabelStyle->Black,AspectRatio->1/GoldenRatio}];
SetOptions[MatrixPlot,BaseStyle->{FontFamily->"Helvetica",FontSize->14,FontWeight->Plain,FontColor->Black}];
SetOptions[MatrixPlot,{FrameTicksStyle->None,FrameTicks->Automatic,FrameStyle->Thickness[Medium],LabelStyle->Black,AspectRatio->1}];

LaunchKernels[8]

Plotting=True;
SetDirectory[NotebookDirectory[]];

frame[legend_]:=Framed[legend,FrameStyle->Black,RoundingRadius->0,FrameMargins->0,Background->White]

dRdt=(R0-R)δ -r U[R] X1 - μX[X1]ϕ r U[R] X2
dX1dt=r U[R] X1-f[X1, X2] X1 Y -δ X1
dX2dt=ϵX[X1]μX[X1] ϕ r U[R] X2-μY[X1] p f[X1, X2]X2 Y -δ X2
dYdt= f[X1, X2] X1 Y+ϵY[X1] μY[X1] p f[X1, X2] X2 Y -δ Y

μX[0]=μX0;
μY[0]=μY0;
ϵX[0]=ϵX0;
ϵY[0]=ϵY0;

f[X1,X2]=g/(H+X1+p X2);

dYdt/.X1->0

solX2 =X2/.First[Solve[0==dYdt/.X1->0,X2]]

solY=Y/.First[Solve[0==dX2dt/.X1->0,Y]]/.X2->solX2

assume specific prey limitation function

U[R]=R/(K+R)

solR=R/.First[Solve[0==dRdt/.{X1->0, X2->solX2},R]]

invGRX1=Simplify[1/X1 dX1dt]/.Y->solY/.X2->solX2/.X1->0/.R->solR

Determine boundaries of different dynamics: for the resident state of no predator

ϵX0predmin=ϵX0/.First[Solve[0==solY/.R->solR, ϵX0]]
μX0predmin=μX0/.First[Solve[0==solY/.R->solR, μX0]]
ϵY0predmin=ϵY0/.First[Solve[0==solY/.R->solR, ϵY0]]
μY0predmin=μY0/.First[Solve[0==solY/.R->solR, μY0]]

now for the resident state of no predator

solRnopred=R/.First[Solve[0==dX2dt/.Y->0/.X1->0,R]]

solX2nopred=Simplify[X2/.First[Solve[0==dRdt/.X1->0/.R->solRnopred,X2]]]

invGRX1nopred=Simplify[1/X1 dX1dt/.Y->0/.R->solRnopred]

now for the resident state of no predator and no prey

ϵX0predminpreymin=ϵX0/.Last[Solve[0==solX2nopred,ϵX0]]

μX0predminpreymin=μX0/.Last[Solve[0==solX2nopred,μX0]]

solRnoprednoprey=R/.First[Solve[0==dRdt/.Y->0/.X1->0/.X2->0,R]]

invGRX1noprednoprey=Simplify[1/X1 dX1dt/.Y->0/.X2->1/.R->solRnoprednoprey]

Numerical determination of resident state is necessary for limit cycle resident dynamics, as there are many limit cycles

essentials={
μX[X1[t]]->μX0+(1-μX0)X1[t]/(b R0+X1[t]),
μY[X1[t]]->μY0+(1-μY0)X1[t]/(b R0+X1[t]),
ϵX[X1[t]]->ϵX0+(1-ϵX0)X1[t]/(b R0+X1[t]),
ϵY[X1[t]]->ϵY0+(1-ϵY0)X1[t]/(b R0+X1[t])}

invGR[paramsp_,μX0p_,μY0p_,ϵX0p_,ϵY0p_]:=
Module[{params=paramsp,μX0t=μX0p,μY0t=μY0p,ϵX0t=ϵX0p,ϵY0t=ϵY0p},
tf=2000;
tshift=tf-200;
stepsize=0.01;
s=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==1,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->μX0t, μY0->μY0t,ϵX0->ϵX0t,ϵY0->ϵY0t},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
convccheck1=Abs[D[Normal[LinearModelFit[Table[{t,Y[t]/.s[[1]]},{t,tshift,tf,stepsize}],x,x]],x]];
convccheck2=Abs[D[Fit[MovingMap[Variance,Table[Y[t]/.s[[1]],{t,tshift,tf,stepsize}],IntegerPart[5/stepsize]],{1,x},x],x]];
If[convccheck1>10^-3||convccheck2>10^-3,Print["Questionable convergence"]];
tt=FindPeaks[Table[Y[t]/.s[[1]],{t,tshift,tf,stepsize}]];
period=If[Length[tt]>3,Mean[Differences[tt[[2;;Length[tt]-1,1]]]]*stepsize,10];
Rfunc=R/.s[[1]];
X2func=X2/.s[[1]];
Yfunc=Y/.s[[1]];
ttt=Transpose[Table[{Rfunc[t], X2func[t],Yfunc[t]},{t,tshift, tshift+period, stepsize}]];
{Mean[Table[Simplify[1/X1 dX1dt]/.X1->0/.params/.{R-> Rfunc[t], X2-> X2func[t], Y-> Yfunc[t]},{t,tshift, tshift+period, stepsize}]],
Table[MinMax[ttt[[i]]],{i,1,3}]}
]

Plotting function

InvGRPlot :=Grid[{{,
Graphics[Style[Text["           Uptake rate"],16], AspectRatio->1/10, ImageSize->300],
 Graphics[Style[Text["          Conversion efficiency"],16], AspectRatio->1/10, ImageSize->300]},
{Graphics[Rotate[Style[Text["Essential for\ncompetitor"],16],0 Degree], AspectRatio->1/5, ImageSize->100],
Show[RegionPlot[(e<1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1})&&(e>1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}),{e,0, 1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->50],
RegionPlot[e<1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµX,-10},{bifurcµX,10}}]}],
Plot[invGRX1noprednoprey/.params/.{μY0->1,ϵX0->1,ϵY0->1}/.μX0->1-e,{e,1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1},1}, PlotStyle->{Gray,Dashed, Thickness[0.011]}],
Plot[invGRX1nopred/.params/.{μY0->1,ϵX0->1,ϵY0->1}/.μX0->1-e,
{e,1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1},
1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}],
Plot[invGRX1/.params/.{μY0->1,ϵX0->1,ϵY0->1}/.μX0->1-e,{e,0,1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}], 
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[4]],Thickness[0.015]}, Joined->True],
 PlotRange->{{0,1},{-1,1}},FrameLabel->{"","Invasion growth rate"}, ImageSize->300,
Epilog-> Text[Style["(a)",20],Scaled[{0.08,0.9}]]],
Show[RegionPlot[(e<1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1})&& e>1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->40],
RegionPlot[e<1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcϵX,-10},{bifurcϵX,10}}]}],
Plot[invGRX1noprednoprey/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,{e,1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1},1}, PlotStyle->{Gray,Dashed, Thickness[0.011]}],
Plot[invGRX1nopred/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
{e,1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1},
1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}],
Plot[invGRX1/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,{e,0,1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}], 
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[4]],Thickness[0.015]}, Joined->True],
PlotRange->{{0,1},{-1,1}}, FrameLabel->{"",""}, ImageSize->300,
Epilog-> Text[Style["(b)",20],Scaled[{0.08,0.9}]]]},
{Graphics[Rotate[Style[Text["Essential for\npredator"],16],0 Degree], AspectRatio->1/5, ImageSize->100],Show[RegionPlot[e>1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµY,-10},{bifurcµY,10}}]}],
Plot[invGRX1nopred/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
{e,1,1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}],
Plot[invGRX1/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,{e,0,1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}], 
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,1]]}]], PlotStyle->{ColorData[97,"ColorList"][[4]],Thickness[0.015]}, Joined->True,PlotRange->All],
PlotRange->{{0,1},{-1,1}}, FrameLabel->{"Essentiality e","Invasion growth rate"}, ImageSize->300,
Epilog-> Text[Style["(c)",20],Scaled[{0.08,0.9}]]],
Show[RegionPlot[e>1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcϵY,-10},{bifurcϵY,10}}]}],
Plot[invGRX1nopred/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
{e,1,1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}],
Plot[invGRX1/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,{e,0,1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1}}, PlotStyle->{Gray,Dashed, Thickness[0.011]}], 
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,1]]}]],  PlotStyle->{ColorData[97,"ColorList"][[4]],Thickness[0.015]}, Joined->True],
PlotRange->{{0,1},{-1,1}}, FrameLabel->{"Essentiality e",""}, ImageSize->300,
Epilog-> Text[Style["(d)",20],Scaled[{0.08,0.9}]]]}}, Spacings->0]

Put in parameters

First food web scenario

params={R0->1,δ->1,r->2,ϕ->1.05,K->0.1,g->1.5,p->0.8,H->0.1, b->0.5}

invGRX1µX=ParallelTable[{0.01i, invGR[params,0.01i,1,1,1]},{i,0,100}];

invGRX1µY=ParallelTable[{0.01i, invGR[params,1,0.01i,1,1]},{i,0,100}];

invGRX1ϵX=ParallelTable[{0.01i, invGR[params,1,1,0.01i,1]},{i,0,100}];

invGRX1ϵY=ParallelTable[{0.01i, invGR[params,1,1,1,0.01i]},{i,0,100}];

bifurcµX=1-invGRX1µX[[First[FirstPosition[Table[First[Differences[invGRX1µX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µX]}],True]],1]];

bifurcµY=1-invGRX1µY[[First[FirstPosition[Table[First[Differences[invGRX1µY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µY]}],True]],1]];

bifurcϵX=1-invGRX1ϵX[[First[FirstPosition[Table[First[Differences[invGRX1ϵX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵX]}],True]],1]];

bifurcϵY=1-invGRX1ϵY[[First[FirstPosition[Table[First[Differences[invGRX1ϵY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵY]}],True]],1]];

pinvgr1=InvGRPlot

If[Plotting==True,Export[StringJoin["Fig_InvGr_1.pdf"],pinvgr1]]
If[Plotting==True,Export[StringJoin["Fig_InvGr_1.png"],pinvgr1, ImageResolution->300]]

pstates1=Grid[{{,
Graphics[Style[Text["           Uptake rate"],16], AspectRatio->1/10, ImageSize->300],
 Graphics[Style[Text["          Conversion efficiency"],16], AspectRatio->1/10, ImageSize->300]},
{Graphics[Rotate[Style[Text["Essential for\ncompetitor"],16],0 Degree], AspectRatio->1/5, ImageSize->100],
Show[
RegionPlot[(e<1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1})&&(e>1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}),{e,0, 1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->40],
RegionPlot[e<1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµX,-10},{bifurcµX,10}}]}],
Plot[{solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solX2/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solY/.R->solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e},{e,0,1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All,
PlotLegends->Placed[
LineLegend[{ColorData[97,"ColorList"][[1]],
ColorData[97,"ColorList"][[2]],
ColorData[97,"ColorList"][[3]]},
{"Resource","Competitor","Predator"},
LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.8,0.5}]],
 PlotRange->{{0,1},{0,1.1}},FrameLabel->{"","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(a)",20],Scaled[{0.08,0.9}]]],
Show[
RegionPlot[(e<1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1})&&
e>1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcϵX,-10},{bifurcϵX,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e},{e,0,1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"",""}, ImageSize->300,
Epilog-> Text[Style["(b)",20],Scaled[{0.08,0.9}]]]},
{Graphics[Rotate[Style[Text["Essential for\npredator"],16],0 Degree], AspectRatio->1/5, ImageSize->100],Show[RegionPlot[e>1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµY,-10},{bifurcµY,10}}]}],
Plot[{solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solX2/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solY/.R->solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e},{e,0,1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(c)",20],Scaled[{0.08,0.9}]]],
Show[RegionPlot[e<1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
RegionPlot[e>1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
Graphics[{Dotted,Line[{{bifurcϵY,-10},{bifurcϵY,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e},{e,0,1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,2]]}]],
PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e",""}, ImageSize->300,
Epilog-> Text[Style["(d)",20],Scaled[{0.08,0.9}]]]}}, Spacings->0]

If[Plotting==True,Export[StringJoin["Fig_States_1.pdf"],pstates1]]
If[Plotting==True,Export[StringJoin["Fig_States_1.png"],pstates1, ImageResolution->300]]

tf=1000;
tshift=tf-150;
stepsize=0.01;
s01=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->0.9, μY0->1,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s1=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,2]],
Y[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,3]]}/.essentials/.params/.{μX0->0.9, μY0->1,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s02=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->0.7, μY0->1,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,2]],
Y[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,3]]}/.essentials/.params/.{μX0->0.7, μY0->1,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];

pdyn1=GraphicsRow[{Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s1],{t,0,100},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time","Biomass"}, 
Prolog-> Text[Style["(a)",20],Scaled[{0.07,0.91}]],
Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.1"]],Background->White], Scaled[{0.865,0.11}]],
AspectRatio->1/GoldenRatio, ImageSize->350] ,
Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2],{t,0,100},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time",""},
Prolog-> Text[Style["(b)",20],Scaled[{0.07,0.91}]],Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.3"]],Background->White], Scaled[{0.89,0.11}]],PlotLegends->Placed[LineLegend[{"Resource","Competitor","Predator", "Focal Prey"},LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.84,0.73}],
AspectRatio->1/GoldenRatio,ImageSize->350]}, ImageSize->700, Spacings->{0,0}]

If[Plotting==True,Export[StringJoin["Fig_Dyn_1.pdf"],pdyn1]]
If[Plotting==True,Export[StringJoin["Fig_Dyn_1.png"],pdyn1, ImageResolution->300]]

pMuX1=Plot[{μX[X1[t]]/. essentials/.μX0->0.7/.params,
μX[X1[t]]/. essentials/.μX0->0.9/.params},{X1[t],0,1},PlotRange->{Automatic,{-0.05,1.05}},PlotStyle->{{Thickness[0.015],Dashing[0.03],ColorData[97,"ColorList"][[1]]},{Thickness[0.015],ColorData[97,"ColorList"][[1]]}}, Frame->True, FrameLabel->{ToString[Subscript["X","1"],StandardForm],ToString[Subscript["μ","X"],StandardForm]}, 
PlotLegends->Placed[{"e = 0.1","e = 0.3"},{0.8,0.5}]]

If[Plotting==True,Export[StringJoin["Fig_MuX_1.png"],pMuX1, ImageResolution->300]]

Second food web scenario

params={R0->1,δ->1,r->2.0,ϕ-> 0.95, K->0.1,g->1.5,p->0.8,H->0.1, b->0.5}

invGRX1µX=ParallelTable[{0.01i, invGR[params,0.01i,1,1,1]},{i,0,100}];

invGRX1µY=ParallelTable[{0.01i, invGR[params,1,0.01i,1,1]},{i,0,100}];

invGRX1ϵX=ParallelTable[{0.01i, invGR[params,1,1,0.01i,1]},{i,0,100}];

invGRX1ϵY=ParallelTable[{0.01i, invGR[params,1,1,1,0.01i]},{i,0,100}];

bifurcµX=1-invGRX1µX[[First[FirstPosition[Table[First[Differences[invGRX1µX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µX]}],True]],1]];

bifurcµY=1-invGRX1µY[[First[FirstPosition[Table[First[Differences[invGRX1µY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µY]}],True]],1]];

bifurcϵX=1-invGRX1ϵX[[First[FirstPosition[Table[First[Differences[invGRX1ϵX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵX]}],True]],1]];

bifurcϵY=1-invGRX1ϵY[[First[FirstPosition[Table[First[Differences[invGRX1ϵY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵY]}],True]],1]];

pinvgr2=InvGRPlot

If[Plotting==True,Export[StringJoin["Fig_InvGr_2.pdf"],pinvgr2]]
If[Plotting==True,Export[StringJoin["Fig_InvGr_2.png"],pinvgr2, ImageResolution->300]]

pstates2=Grid[{{,
Graphics[Style[Text["           Uptake rate"],16], AspectRatio->1/10, ImageSize->300],
 Graphics[Style[Text["          Conversion efficiency"],16], AspectRatio->1/10, ImageSize->300]},
{Graphics[Rotate[Style[Text["Essential for\ncompetitor"],16],0 Degree], AspectRatio->1/5, ImageSize->100],
Show[
RegionPlot[(e<1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1})&&(e>1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}),{e,0, 1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->40],
RegionPlot[e<1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµX,-10},{bifurcµX,10}}]}],
Plot[{solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solX2/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solY/.R->solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e},{e,0,1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All,
PlotLegends->Placed[
LineLegend[{ColorData[97,"ColorList"][[1]],
ColorData[97,"ColorList"][[2]],
ColorData[97,"ColorList"][[3]]},
{"Resource","Competitor","Predator"},
LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.8,0.5}]],
 PlotRange->{{0,1},{0,1.1}},FrameLabel->{"","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(a)",20],Scaled[{0.08,0.9}]]],
Show[
RegionPlot[(e<1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1})&&
e>1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->40],
RegionPlot[e<1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True,PlotPoints->40],
Graphics[{Dotted,Line[{{bifurcϵX,-10},{bifurcϵX,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e},{e,0,1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"",""}, ImageSize->300,
Epilog-> Text[Style["(b)",20],Scaled[{0.08,0.9}]]]},
{Graphics[Rotate[Style[Text["Essential for\npredator"],16],0 Degree], AspectRatio->1/5, ImageSize->100],Show[RegionPlot[e>1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµY,-10},{bifurcµY,10}}]}],
Plot[{solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solX2/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solY/.R->solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e},{e,0,1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(c)",20],Scaled[{0.08,0.9}]]],
Show[RegionPlot[e<1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
RegionPlot[e>1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
Graphics[{Dotted,Line[{{bifurcϵY,-10},{bifurcϵY,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e},{e,0,1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,2]]}]],
PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e",""}, ImageSize->300,
Epilog-> Text[Style["(d)",20],Scaled[{0.08,0.9}]]]}}, Spacings->0]

If[Plotting==True,Export[StringJoin["Fig_States_2.pdf"],pstates2]]
If[Plotting==True,Export[StringJoin["Fig_States_2.png"],pstates2, ImageResolution->300]]

tf=1000;
tshift=tf-150;
stepsize=0.01;
s01=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->1, μY0->0.9,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s1=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,2]],
Y[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,3]]}/.essentials/.params/.{μX0->1, μY0->0.9,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s02=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2a=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,2]],
Y[0]==0}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2b=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,2]],
Y[0]==10^-3}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];

pdyn2=GraphicsRow[{Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s1],{t,0,200},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time","Biomass"}, 
Prolog-> Text[Style["(a)",20],Scaled[{0.07,0.91}]],
Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.1"]],Background->White], Scaled[{0.865,0.11}]],
PlotLegends->Placed[LineLegend[{"Resource","Competitor","Predator", "Focal Prey"},LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.84,0.73}],
AspectRatio->1/GoldenRatio, ImageSize->350] ,
Show[Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2a],{t,0,300},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time",""},
Prolog-> Text[Style["(b)",20],Scaled[{0.07,0.91}]],Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.3"]],Background->White], Scaled[{0.89,0.3}]],
AspectRatio->1/GoldenRatio, ImageSize->350],
Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2b],{t,0,1000}, PlotStyle->{{Thickness[0.012],Dashing[0.025,0.5,"Round"]}}]]}, ImageSize->700, Spacings->{0,0}]

If[Plotting==True,Export[StringJoin["Fig_Dyn_2.pdf"],pdyn2]]
If[Plotting==True,Export[StringJoin["Fig_Dyn_2.png"],pdyn2, ImageResolution->300]]

Third food web scenario

params={R0->1,δ->1,r->2.0,ϕ-> 0.95,K->0.1,g->1.5,p->1.2 ,H->0.1, b->0.5}

invGRX1µX=ParallelTable[{0.01i, invGR[params,0.01i,1,1,1]},{i,0,100}];

invGRX1µY=ParallelTable[{0.01i, invGR[params,1,0.01i,1,1]},{i,0,100}];

invGRX1ϵX=ParallelTable[{0.01i, invGR[params,1,1,0.01i,1]},{i,0,100}];

invGRX1ϵY=ParallelTable[{0.01i, invGR[params,1,1,1,0.01i]},{i,0,100}];

bifurcµX=1-invGRX1µX[[First[FirstPosition[Table[First[Differences[invGRX1µX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µX]}],True]],1]];

bifurcµY=1-invGRX1µY[[First[FirstPosition[Table[First[Differences[invGRX1µY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µY]}],True]],1]];

bifurcϵX=1-invGRX1ϵX[[First[FirstPosition[Table[First[Differences[invGRX1ϵX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵX]}],True]],1]];

bifurcϵY=1-invGRX1ϵY[[First[FirstPosition[Table[First[Differences[invGRX1ϵY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵY]}],True]],1]];

pinvgr3=InvGRPlot

If[Plotting==True,Export[StringJoin["Fig_InvGr_3.pdf"],pinvgr3]]
If[Plotting==True,Export[StringJoin["Fig_InvGr_3.png"],pinvgr3, ImageResolution->300]]

pstates3=Grid[{{,
Graphics[Style[Text["           Uptake rate"],16], AspectRatio->1/10, ImageSize->300],
 Graphics[Style[Text["          Conversion efficiency"],16], AspectRatio->1/10, ImageSize->300]},
{Graphics[Rotate[Style[Text["Essential for\ncompetitor"],16],0 Degree], AspectRatio->1/5, ImageSize->100],
Show[
RegionPlot[(e<1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1})&&(e>1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}),{e,0, 1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->50],
RegionPlot[e<1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµX,-10},{bifurcµX,10}}]}],
Plot[{solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solX2/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solY/.R->solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e},{e,0,1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All,
PlotLegends->Placed[
LineLegend[{ColorData[97,"ColorList"][[1]],
ColorData[97,"ColorList"][[2]],
ColorData[97,"ColorList"][[3]]},
{"Resource","Competitor","Predator"},
LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.8,0.5}]],
 PlotRange->{{0,1},{0,1.1}},FrameLabel->{"","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(a)",20],Scaled[{0.08,0.9}]]],
Show[
RegionPlot[(e<1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1})&&
e>1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->40],
RegionPlot[e<1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True,PlotPoints->40],
Graphics[{Dotted,Line[{{bifurcϵX,-10},{bifurcϵX,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e},{e,0,1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"",""}, ImageSize->300,
Epilog-> Text[Style["(b)",20],Scaled[{0.08,0.9}]]]},
{Graphics[Rotate[Style[Text["Essential for\npredator"],16],0 Degree], AspectRatio->1/5, ImageSize->100],Show[RegionPlot[e>1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµY,-10},{bifurcµY,10}}]}],
Plot[{solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solX2/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solY/.R->solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e},{e,0,1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(c)",20],Scaled[{0.08,0.9}]]],
Show[RegionPlot[e<1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
RegionPlot[e>1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
Graphics[{Dotted,Line[{{bifurcϵY,-10},{bifurcϵY,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e},{e,0,1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,2]]}]],
PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e",""}, ImageSize->300,
Epilog-> Text[Style["(d)",20],Scaled[{0.08,0.9}]]]}}, Spacings->0]

If[Plotting==True,Export[StringJoin["Fig_States_3.pdf"],pstates3]]
If[Plotting==True,Export[StringJoin["Fig_States_3.png"],pstates3, ImageResolution->300]]

tf=1000;
tshift=tf-150;
stepsize=0.01;
s01=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->1, μY0->0.9,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s1=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,2]],
Y[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,3]]}/.essentials/.params/.{μX0->1, μY0->0.9,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s02=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2a=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,2]],
Y[0]==0}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2b=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,2]],
Y[0]==10^-3}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];

GraphicsRow[{Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s1],{t,0,200},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time","Biomass"}, 
Prolog-> Text[Style["(a)",20],Scaled[{0.07,0.91}]],
Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.1"]],Background->White], Scaled[{0.865,0.11}]],
AspectRatio->1/GoldenRatio, ImageSize->350] ,
Show[Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2a],{t,0,200},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time",""},
Prolog-> Text[Style["(b)",20],Scaled[{0.07,0.91}]],Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.3"]],Background->White], Scaled[{0.89,0.11}]],PlotLegends->Placed[LineLegend[{"Resource","Competitor","Predator", "Focal Prey"},LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.84,0.73}],
AspectRatio->1/GoldenRatio, ImageSize->350],
Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2b],{t,0,1000}, PlotStyle->{{Thickness[0.01],Dashing[0.025,0.5,"Round"]}}]]}, ImageSize->700, Spacings->{0,0}]

Fourth food web scenario

params={R0->1,δ->1,r->2.0,ϕ-> 1.05, ϕ2->1,K->0.1,g->1.5,p->1.2,H->0.1, b->0.5}

invGRX1µX=ParallelTable[{0.01i, invGR[params,0.01i,1,1,1]},{i,0,100}];

invGRX1µY=ParallelTable[{0.01i, invGR[params,1,0.01i,1,1]},{i,0,100}];

invGRX1ϵX=ParallelTable[{0.01i, invGR[params,1,1,0.01i,1]},{i,0,100}];

invGRX1ϵY=ParallelTable[{0.01i, invGR[params,1,1,1,0.01i]},{i,0,100}];

bifurcµX=1-invGRX1µX[[First[FirstPosition[Table[First[Differences[invGRX1µX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µX]}],True]],1]];

bifurcµY=1-invGRX1µY[[First[FirstPosition[Table[First[Differences[invGRX1µY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1µY]}],True]],1]];

bifurcϵX=1-invGRX1ϵX[[First[FirstPosition[Table[First[Differences[invGRX1ϵX[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵX]}],True]],1]];

bifurcϵY=1-invGRX1ϵY[[First[FirstPosition[Table[First[Differences[invGRX1ϵY[[i,2,2,3]]]]>10^-5,{i,1,Length[invGRX1ϵY]}],True]],1]];

pinvgr4=InvGRPlot

If[Plotting==True,Export[StringJoin["Fig_InvGr_4.pdf"],pinvgr4]]
If[Plotting==True,Export[StringJoin["Fig_InvGr_4.png"],pinvgr4, ImageResolution->300]]

pstates4=Grid[{{,
Graphics[Style[Text["           Uptake rate"],16], AspectRatio->1/10, ImageSize->300],
 Graphics[Style[Text["          Conversion efficiency"],16], AspectRatio->1/10, ImageSize->300]},
{Graphics[Rotate[Style[Text["Essential for\ncompetitor"],16],0 Degree], AspectRatio->1/5, ImageSize->100],
Show[
RegionPlot[(e<1-μX0predminpreymin/.params/.{μY0->1,ϵX0->1,ϵY0->1})&&(e>1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}),{e,0, 1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->50],
RegionPlot[e<1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµX,-10},{bifurcµX,10}}]}],
Plot[{solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solX2/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e,
solY/.R->solR/.params/.{ϵX0->1,μY0->1,ϵY0->1}/.μX0->1-e},{e,0,1-μX0predmin/.params/.{μY0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µX[[;;,1]]},{invGRX1µX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All,
PlotLegends->Placed[
LineLegend[{ColorData[97,"ColorList"][[1]],
ColorData[97,"ColorList"][[2]],
ColorData[97,"ColorList"][[3]]},
{"Resource","Competitor","Predator"},
LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.8,0.5}]],
 PlotRange->{{0,1},{0,1.1}},FrameLabel->{"","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(a)",20],Scaled[{0.08,0.9}]]],
Show[
RegionPlot[(e<1-ϵX0predminpreymin/.params/.{μX0->1,μY0->1,ϵY0->1})&&
e>1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True,PlotPoints->40],
RegionPlot[e<1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True,PlotPoints->40],
Graphics[{Dotted,Line[{{bifurcϵX,-10},{bifurcϵX,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵY0->1}/.ϵX0->1-e},{e,0,1-ϵX0predmin/.params/.{μX0->1,μY0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵX[[;;,1]]},{invGRX1ϵX[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"",""}, ImageSize->300,
Epilog-> Text[Style["(b)",20],Scaled[{0.08,0.9}]]]},
{Graphics[Rotate[Style[Text["Essential for\npredator"],16],0 Degree], AspectRatio->1/5, ImageSize->100],Show[RegionPlot[e>1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
RegionPlot[e<1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
Graphics[{Dotted,Line[{{bifurcµY,-10},{bifurcµY,10}}]}],
Plot[{solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solX2/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e,
solY/.R->solR/.params/.{μX0->1,ϵX0->1,ϵY0->1}/.μY0->1-e},{e,0,1-μY0predmin/.params/.{μX0->1,ϵX0->1,ϵY0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1µY[[;;,1]]},{invGRX1µY[[;;,2,2,3,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e","Biomass"}, ImageSize->300,
Epilog-> Text[Style["(c)",20],Scaled[{0.08,0.9}]]],
Show[RegionPlot[e<1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->Lighter[Gray],Axes->True],
RegionPlot[e>1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1},{e,0,1},{x,-10,10},BoundaryStyle ->None,PlotStyle->LightGray,Axes->True],
Graphics[{Dotted,Line[{{bifurcϵY,-10},{bifurcϵY,10}}]}],
Plot[{solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solX2/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e,
solY/.R->solR/.params/.{μX0->1,μY0->1,ϵX0->1}/.ϵY0->1-e},{e,0,1-ϵY0predmin/.params/.{μX0->1,μY0->1,ϵX0->1}}, PlotStyle->Dashed],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,1,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[1]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,2,2]]}]],PlotStyle->{ColorData[97,"ColorList"][[2]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,1]]}]],PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
ListPlot[Transpose[Join[{1-invGRX1ϵY[[;;,1]]},{invGRX1ϵY[[;;,2,2,3,2]]}]],
PlotStyle->{ColorData[97,"ColorList"][[3]],Thickness[0.015]}, Joined->True, PlotRange->All],
PlotRange->{{0,1},{0,1.1}}, FrameLabel->{"Essentiality e",""}, ImageSize->300,
Epilog-> Text[Style["(d)",20],Scaled[{0.08,0.9}]]]}}, Spacings->0]

If[Plotting==True,Export[StringJoin["Fig_States_4.pdf"],pstates4]]
If[Plotting==True,Export[StringJoin["Fig_States_4.png"],pstates4, ImageResolution->300]]

tf=1000;
tshift=tf-150;
stepsize=0.01;
s01=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->1, μY0->0.9,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s1=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,2]],
Y[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s01]/.t->tf)[[1,3]]}/.essentials/.params/.{μX0->1, μY0->0.9,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s02=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==0.8,
X1[0]==0,
X2[0]==10^-1,
Y[0]==10^-1}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2a=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,2]],
Y[0]==0}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];
s2b=NDSolve[{R'[t]==(dRdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X1'[t]==(dX1dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
X2'[t]==(dX2dt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
Y'[t]==(dYdt/.{R->R[t],X1->X1[t],X2->X2[t],Y->Y[t]}),
R[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,1]],
X1[0]==10^-3,
X2[0]==(Evaluate[{R[t],X2[t],Y[t]}/.s02]/.t->tf)[[1,2]],
Y[0]==10^-3}/.essentials/.params/.{μX0->1, μY0->0.7,ϵX0->1,ϵY0->1},{R,X1,X2,Y},{t,0,tf},MaxStepSize->stepsize];

GraphicsRow[{Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s1],{t,0,200},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time","Biomass"}, 
Prolog-> Text[Style["(a)",20],Scaled[{0.07,0.91}]],
Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.1"]],Background->White], Scaled[{0.865,0.11}]],
AspectRatio->1/GoldenRatio, ImageSize->350] ,
Show[Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2a],{t,0,200},PlotRange->{Automatic,{-0.05,1.19}}, Frame->True, FrameLabel->{"Time",""},
Prolog-> Text[Style["(b)",20],Scaled[{0.07,0.91}]],Epilog->Inset[Framed[Style[StringJoin[ToString[Subscript["e",""],StandardForm]," = 0.3"]],Background->White], Scaled[{0.89,0.11}]],PlotLegends->Placed[LineLegend[{"Resource","Competitor","Predator", "Focal Prey"},LabelStyle->12,LegendMarkerSize->{10,3},LegendFunction->frame,LegendLayout->{"Column",1}],{0.84,0.73}],
AspectRatio->1/GoldenRatio, ImageSize->350],
Plot[Evaluate[{R[t],X2[t],Y[t],X1[t]}/.s2b],{t,0,200}, PlotStyle->{{Thickness[0.01],Dashing[0.025,0.5,"Round"]}}]]}, ImageSize->700, Spacings->{0,0}]