Provision of essential resources as a persistence strategy in food webs

Michael Raatz (mraatz@evolbio.mpg.de)

Notebook #2 - calculating Fig. 3

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

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);

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

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_,ϕp_,pp_]:=
Module[{params=paramsp,μX0t=μX0p,μY0t=μY0p,ϵX0t=ϵX0p,ϵY0t=ϵY0p,ϕt=ϕp,pt=pp},
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/.{ϕ->ϕt,p->pt}/.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/.{ϕ->ϕt,p->pt}/.params/.{R-> Rfunc[t], X2-> X2func[t], Y-> Yfunc[t]},{t,tshift, tshift+period, stepsize}]]
]

compEss[invGRList_]:=
Module[{X=invGRList},
npos=Total[Boole[Table[X[[i,2]]>10^-6,{i,Length[X]}]]];
ifunc=Interpolation[X];
If[npos==Length[X],
ec=0,
If[npos==0,
ec=1,
ec=ess/.FindRoot[ifunc[ess]==0,{ess,0.2,0,1}]]];
{ec,ifunc'[ec]}]

options[A_,B_]:=Module[{},If[A&&B,Return[1]];If[!A&&!B,Return[2]];If[!A&&B,Return[3]];If[A&&!B,Return[4]]]

decideEffect[params_,ϕ_,p_]:=
Module[{paramst=params,ϕt=ϕ,pt=p},
(*targetting competitor's uptake rate*)
IG0=invGR[params,1,1,1,1,ϕt,pt]>=0;
IG1μX=invGR[params,0,1,1,1,ϕt,pt]>=0;
IG1ϵX=invGR[params,1,1,0,1,ϕt,pt]>=0;
IG1μY=invGR[params,1,0,1,1,ϕt,pt]>=0;
IG1ϵY=invGR[params,1,1,1,0,ϕt,pt]>=0;
{options[IG0,IG1μX],options[IG0,IG1ϵX],options[IG0,IG1μY],options[IG0,IG1ϵY]}]

potential outcomes:
	A - persistence independent of essentiality (IG(e=0)>0 & IG(e=1)>0)
	B - extinction independent of essentiality (IG(e=0)<0 & IG(e=1)<0)
	C - persistence for high essentiality, extinction for low essentiality (promoting effect) (IG(e=0)<0 & IG(e=1)>0)
	D - persistence for low essentiality, extinction for high essentiality (detrimental effect) (IG(e=0)>0 & IG(e=1)<0)

Scanning (for different ϕ and p values determine where IGR crosses axis)

params={R0->1,δ->1,r->2,K->0.1,g->1.5,H->0.1, b->0.5};

xrow=Table[x,{x,0.75,1.25,0.025}];
ycolumn=Table[x,{x,0.75,1.25,0.025}];

yticks=MapIndexed[{-4+5 #2[[1]],#}&,xrow[[1;; ;;5]]];
xticks=MapIndexed[{-4+5#2[[1]],#}&,ycolumn[[1;; ;;5]]];
frame[legend_] := 
 Framed[legend, FrameStyle -> Black, RoundingRadius -> 10, 
  FrameMargins -> 0, Background -> LightGray]

EffectsTable=ParallelTable[decideEffect[params,ϕ,p],{ϕ,xrow},{p,ycolumn}]

Below is the result of the above computation for the parameters given in the manuscript.

(*EffectsTable={{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{3,3,3,3},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{3,3,3,3},{3,3,3,3},{3,3,3,3},{3,3,3,3},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{3,3,3,3},{3,3,3,3},{3,3,3,3},{3,3,3,3},{3,3,3,3},{3,3,3,3},{3,3,3,3},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}},{{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{3,3,2,3},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1},{1,1,4,1}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{1,1,4,4},{1,1,4,4},{1,1,4,4},{1,1,4,4},{1,1,4,4},{1,1,4,4},{1,1,4,4}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{1,1,4,4},{1,1,4,4},{1,1,4,4}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}},{{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2},{3,3,2,2}}};*)

EffectsμX=EffectsTable[[;;,;;,1]];
EffectsϵX=EffectsTable[[;;,;;,2]];
EffectsμY=EffectsTable[[;;,;;,3]];
EffectsϵY=EffectsTable[[;;,;;,4]];

pParaScan=GraphicsGrid[{{Show[{MatrixPlot[EffectsμX, FrameTicks->{{yticks,None},{xticks,None}}, DataReversed->True,ColorFunction->(ColorData["GrayYellowTones"][Rescale[#,{1,4}]]&), ColorFunctionScaling->False,
PlotLabel->"Essentiality to competitor's uptake rate",
FrameLabel->{"Resource competitiveness ϕ",""},
Epilog-> {Text[Style["(a)",20],Scaled[{0.08,0.92}]],
Text[Style["I",20],Scaled[{0.175,0.63}]],
Text[Style["II",20],Scaled[{0.175,0.45}]],
Text[Style["IV",20],Scaled[{0.92,0.63}]],
Text[Style["III",20],Scaled[{0.92,0.45}]]}],Graphics[{PointSize[0.03],Point[({{0.8,1.05},{0.8,0.95},{1.2,0.95},{1.2,1.05}}-Min[xrow])/Mean[Differences[xrow]]+0.5]}]}],
Show[{MatrixPlot[EffectsϵX, FrameTicks->{{yticks,None},{xticks,None}}, DataReversed->True, ColorFunction->(ColorData["GrayYellowTones"][Rescale[#,{1,4}]]&), ColorFunctionScaling->False,
PlotLabel->"Essentiality to competitor's conversion efficiency",
FrameLabel->{"",""},
Epilog-> {Text[Style["(b)",20],Scaled[{0.08,0.92}]],
Text[Style["I",20],Scaled[{0.175,0.63}]],
Text[Style["II",20],Scaled[{0.175,0.45}]],
Text[Style["IV",20],Scaled[{0.92,0.63}]],
Text[Style["III",20],Scaled[{0.92,0.45}]]}],Graphics[{PointSize[0.03],Point[({{0.8,1.05},{0.8,0.95},{1.2,0.95},{1.2,1.05}}-Min[xrow])/Mean[Differences[xrow]]+0.5]}]}]},
{Show[{MatrixPlot[EffectsμY, FrameTicks->{{yticks,None},{xticks,None}}, DataReversed->True, ColorFunction->(ColorData["GrayYellowTones"][Rescale[#,{1,4}]]&), ColorFunctionScaling->False,
PlotLabel->"Essentiality to predator's uptake rate",
FrameLabel->{"Resource competitiveness ϕ","Vulnerability to predation p"},Epilog-> {Text[Style["(c)",20],Scaled[{0.08,0.92}]],
Text[Style["I",20],Scaled[{0.175,0.63}]],
Text[Style["II",20],Scaled[{0.175,0.45}]],
Text[Style["IV",20],Scaled[{0.92,0.63}]],
Text[Style["III",20],Scaled[{0.92,0.45}]]}],Graphics[{PointSize[0.03],Point[({{0.8,1.05},{0.8,0.95},{1.2,0.95},{1.2,1.05}}-Min[xrow])/Mean[Differences[xrow]]+0.5]}]}],
Show[{MatrixPlot[EffectsϵY, FrameTicks->{{yticks,None},{xticks,None}}, DataReversed->True,ColorFunction->(ColorData["GrayYellowTones"][Rescale[#,{1,4}]]&), ColorFunctionScaling->False,
 PlotLegends->Placed[
SwatchLegend[Array[ColorData["GrayYellowTones"],4,{0,1}],{"None, indep. persisting", "None, indep. extinct", "Promoting", "Detrimental" },LegendLabel->"Effect of essentiality on persistence",LegendFunction->frame,LabelStyle->Directive[14]],
{0.6,0.2}],PlotLabel->"Essentiality to predator's conversion efficiency",
FrameLabel->{"","Vulnerability to predation p"},
Epilog-> {Text[Style["(d)",20],Scaled[{0.08,0.92}]], 
Text[Style["I",20],Scaled[{0.175,0.63}]],
Text[Style["II",20],Scaled[{0.175,0.45}]],
Text[Style["IV",20],Scaled[{0.92,0.63}]],
Text[Style["III",20],Scaled[{0.92,0.45}]]}],Graphics[{PointSize[0.03],Point[({{0.8,1.05},{0.8,0.95},{1.2,0.95},{1.2,1.05}}-Min[xrow])/Mean[Differences[xrow]]+0.5]}]}]}},ImageSize->800, Spacings->{-30,-5}]

If[Plotting==True,Export[StringJoin["Fig_ParaScan.pdf"],pParaScan]]

potential outcomes:
	A - persistence independent of essentiality (IG(e=0)>0 & IG(e=1)>0)
	B - extinction independent of essentiality (IG(e=0)<0 & IG(e=1)<0)
	C - persistence for high essentiality, extinction for low essentiality (promoting effect) (IG(e=0)<0 & IG(e=1)>0)
	D - persistence for low essentiality, extinction for high essentiality (detrimental effect) (IG(e=0)>0 & IG(e=1)<0)