Prey-predator coevolution

(d/(dt))(Subscript[m, i])=r[(Subscript[x, i])](Subscript[m, i]) (1-(((∑_(j=1)^k)a[(Subscript[x, i]),(Subscript[x, j])](Subscript[m, j]))/(K[(Subscript[x, i])])))-(Subscript[m, i])(∑_(j=1)^l)b[(Subscript[x, i]),(Subscript[y, j])](Subscript[n, j])  (i=1,... ,k)  (prey)

(d/(dt))(Subscript[n, i])=(Subscript[n, i])(∑_(j=1)^k)c[(Subscript[x, j]),(Subscript[y, i])]b[(Subscript[x, j]),(Subscript[y, i])](Subscript[m, j])-d[(Subscript[y, i])](Subscript[n, i])           (i=1,... ,l) (predator)

MONOMORPHIC PREY AND PREDATOR POPULATIONS

Resident dynamic s

m=prey pop dens;n=pred pop dens;x=log body size prey;y=log body size pred

dLogm[m_,n_]:=r[x] (1-((m a[x,x])/(K[x])))-n b[x,y];
dLogn[m_,n_]:=m b[x,y] c[x,y]-d[y];

Population equilibrium (= time averages // explain!)

Solve[{dLogm[m,n]==0,dLogn[m,n]==0},{m,n}]

m[{x_,y_}]:=((d[y])/(b[x,y] c[x,y]));
n[{x_,y_}]:=((r[x])/(b[x,y]))-((d[y] r[x] a[x,x])/(K[x] ((b[x,y])^2) c[x,y]));

Parameter values and functions

r[x_]:=1;
K[x_]:=(E^(-(x^2)));
a[x1_,x2_]:=(E^(-α (((x1-x2))^2)));α=0.5;
b[x_,y_]:=β1 (E^(-β2 (((-p+x-y))^2)));β1=100;β2=0.2;p=Log[0.5];(* attack rate *)
c[x_,y_]:=γ (E^(x-y));γ=0.2;(* conversion coefficient *)
d[y_]:=δ1 (E^(y (-δ2)));δ1=1;δ2=1;

Row[{Show[Plot[b[x,0],{x,-5,5},AxesOrigin->{0,0},PlotStyle->Black],
Graphics[{Dashed,Line[{{p,0},{p,b[p,0]}}]}],
Graphics[Text["p",{p,12},Background->White]],Ticks->None,AxesLabel->{"x-y","b[x,y]"},ImageSize->Small],"    ",Show[Plot[c[x,0],{x,-5,5},AxesOrigin->{0,0},PlotStyle->Black],Ticks->None,AxesLabel->{"x-y","c[x,y]"},ImageSize->Small]}]

Coexistence

xmin=-1.5;xmax=2.5;ymin=-3.5;ymax=6;

coexBnd=ContourPlot[If[m[{x,y}]>0&&n[{x,y}]>0,1,-1],{x,xmin,xmax},{y,ymin,ymax},Contours->{0},ContourStyle->{Black},ContourShading->False,PlotPoints->60];

preyDens=DensityPlot[If[m[{x,y}]>0&&n[{x,y}]>0,m[{x,y}],k[x]],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,ColorFunction->"AvocadoColors",PlotPoints->60];

predDens=DensityPlot[If[m[{x,y}]>0&&n[{x,y}]>0,n[{x,y}]],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,ColorFunction->"SolarColors",PlotPoints->60];

optRatio=Plot[x-p,{x,xmin,xmax},PlotStyle->{White,Dashed}];

Row[{Show[preyDens,coexBnd,optRatio,FrameLabel->{"log prey size ⟶","log predator size ⟶"},PlotLabel->"Prey pop dens",ImageSize->Small],"      ",Show[predDens,coexBnd,optRatio,FrameLabel->{"log prey size ⟶","log predator size ⟶"},PlotLabel->"Predator pop dens",ImageSize->Small]}]

Invasion fitness and derivatives

Prey

(Subscript[sPrey, mo])[{x_,y_},xMut_]:=r[xMut] (1-((a[xMut,x] m[{x,y}])/(K[xMut])))-b[xMut,y] n[{x,y}];
(Subscript[dsPrey, mo])[{x_,y_}]:=(∂_xMut)(Subscript[sPrey, mo])[{x,y},xMut]/.{xMut->x};
(Subscript[ddsPrey, mo])[{x_,y_}]:=(∂_(xMut,xMut))(Subscript[sPrey, mo])[{x,y},xMut]/.{xMut->x};

Predator

(Subscript[sPred, mo])[{x_,y_},yMut_]:=b[x,yMut] c[x,yMut] m[{x,y}]-d[yMut];
(Subscript[dsPred, mo])[{x_,y_}]:=(∂_yMut)(Subscript[sPred, mo])[{x,y},yMut]/.{yMut->y};
(Subscript[ddsPred, mo])[{x_,y_}]:=(∂_(yMut,yMut))(Subscript[sPred, mo])[{x,y},yMut]/.{yMut->y};

Isoclines plot

Prey

solid=uninvadable (ES);dashed=invadable (NES)

preyES=ContourPlot[If[n[{x,y}]>0&&(Subscript[ddsPrey, mo])[{x,y}]≤0,(Subscript[dsPrey, mo])[{x,y}]],{x,xmin,xmax},{y,ymin,ymax},Contours->{0},ContourShading->False,ContourStyle->{Darker[Green],Thick},PlotPoints->30];
preyNES=ContourPlot[If[n[{x,y}]>0&&(Subscript[ddsPrey, mo])[{x,y}]>0,(Subscript[dsPrey, mo])[{x,y}]],{x,xmin,xmax},{y,ymin,ymax},Contours->{0},ContourShading->False,ContourStyle->{Darker[Green],Thick,Dashed},PlotPoints->30];

Predator

solid=uninvadable (ES);dashed=invadable (NES)

predES=ContourPlot[If[n[{x,y}]>0&&(Subscript[ddsPred, mo])[{x,y}]≤0,(Subscript[dsPred, mo])[{x,y}]],{x,xmin,xmax},{y,ymin,ymax},Contours->{0},ContourShading->False,ContourStyle->{Red,Thick},PlotPoints->30];
predNES=ContourPlot[If[n[{x,y}]>0&&(Subscript[ddsPred, mo])[{x,y}]>0,(Subscript[dsPred, mo])[{x,y}]],{x,xmin,xmax},{y,ymin,ymax},Contours->{0},ContourShading->False,ContourStyle->{Red,Thick,Dashed},PlotPoints->30];

Singularity

singPoint={x,y}/.FindRoot[{(Subscript[dsPrey, mo])[{x,y}]==0,(Subscript[dsPred, mo])[{x,y}]==0},{x,p},{y,0}];
sngPnt=Graphics[{PointSize->Large,Point[singPoint]}];

Isocline plot

green=prey isocline;red=predator isocline;
solid=uninvadable (ES);dashed=invadable (NES)

Show[preyES,preyNES,predES,predNES,coexBnd,sngPnt]

Stability

Jacobi matrix {{a11,a12},{a21,a22}}

a11[{x_,y_}]:=(∂_(X2,X2))(Subscript[sPrey, mo])[{x,y},X2]+(∂_(X1,X2))(Subscript[sPrey, mo])[{X1,y},X2]/.{X1->x,X2->x};
a12[{x_,y_}]:=(∂_(X,Y))(Subscript[sPrey, mo])[{x,Y},X]/.{X->x,Y->y};
a21[{x_,y_}]:=(∂_(X,Y))(Subscript[sPred, mo])[{X,y},Y]/.{X->x,Y->y};
a22[{x_,y_}]:=(∂_(Y2,Y2))(Subscript[sPred, mo])[{x,y},Y2]+(∂_(Y1,Y2))(Subscript[sPred, mo])[{x,Y1},Y2]/.{Y1->y,Y2->y};

Total stability

totStab=DensityPlot[If[n[{x,y}]>0&&a11[{x,y}]<0∧a22[{x,y}]<0∧a11[{x,y}] a22[{x,y}]>Abs[a12[{x,y}] a21[{x,y}]],Abs[(Subscript[dsPred, mo])[{x,y}]]+Abs[(Subscript[dsPrey, mo])[{x,y}]]],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,ColorFunction->"FuchsiaTones",PlotPoints->30];

Show[totStab,preyES,preyNES,predES,predNES,coexBnd,sngPnt]

Conclusion: the singular point is NOT ABSOLUTELY STABLE.

Strong

strStab=DensityPlot[If[n[{x,y}]>0∧a11[{x,y}]<0∧a22[{x,y}]<0∧a11[{x,y}] a22[{x,y}]>a12[{x,y}] a21[{x,y}],Abs[(Subscript[dsPred, mo])[{x,y}]]+Abs[(Subscript[dsPrey, mo])[{x,y}]]],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,ColorFunction->"CherryTones",PlotPoints->30];

Show[strStab,preyES,preyNES,predES,predNES,coexBnd,sngPnt]

Conclusion: the singular point is NOT STRONGLY STABLE.

stability Weak

wkStab=DensityPlot[If[n[{x,y}]>0∧(a11[{x,y}] <0∨a22[{x,y}] <0)∧a11[{x,y}] a22[{x,y}]>a12[{x,y}] a21[{x,y}],Abs[(Subscript[dsPred, mo])[{x,y}]]+Abs[(Subscript[dsPrey, mo])[{x,y}]]],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,ColorFunction->"CoffeeTones",PlotPoints->50];

Show[wkStab,preyES,preyNES,predES,predNES,coexBnd,sngPnt]

Conclusion:the singular point is WEAKLY STABLE

Canonical equation

Mutation probability per birth event(μ) and mutation variance (σ2)

μPrey=(10^(-6));σ2Prey=(10^(-6));
μPred=(10^(-6));σ2Pred=(10^(-6));

μPrey=1;σ2Prey=1;
μPred=1;σ2Pred=1;
(* for the shape of the orbits only the relative values matter *)

Deterministic drift

(Subscript[drift, mo])[{x_,y_}]:={(1/2) μPrey σ2Prey m[{x,y}] (Subscript[dsPrey, mo])[{x,y}],(1/2) μPred σ2Pred n[{x,y}] (Subscript[dsPred, mo])[{x,y}]};

CE plot Stream

CEstream=StreamPlot[If[n[{x,y}]>0,(Subscript[drift, mo])[{x,y}],{0,0}],{x,xmin,xmax},{y,ymin,ymax},StreamPoints->Fine];

Show[preyES,preyNES,predES,predNES,coexBnd,CEstream,sngPnt]

Deterministic orbit (Euler method)

v0={-1,-2};t0=0;t∞=20000;Δt=2;data={};
v=v0;t=t0;
While[t≤t∞&&n[v]>0,data=Join[data,{Append[v,t]}];v=Δt (Subscript[drift, mo])[v]+v;t=t+Δt;];

... projected in the (x,y)-plane

CEorbit=ListPlot[data[[All,{1,2}]],PlotStyle->{Gray,Thickness[.02]},Joined->True];

Show[coexBnd,CEorbit,preyES,preyNES,predES,predNES,CEstream,sngPnt]

... as evolutionary tree

xtree=ListPlot[data[[All,{1,3}]],PlotStyle->{Darker[Green],Thick},Joined->True];
ytree=ListPlot[data[[All,{2,3}]],PlotStyle->{Red,Thick},Joined->True];
Show[xtree,ytree,AxesLabel->{"log body size","evo time"}]

Different parameter set

μPrey=1;
σ2Prey=0.1;
μPred=1;
σ2Pred=1;

CE plot Stream

CEstream=StreamPlot[If[n[{x,y}]>0,(Subscript[drift, mo])[{x,y}],{0,0}],{x,xmin,xmax},{y,ymin,ymax},StreamPoints->Fine];

Show[preyES,preyNES,predES,predNES,coexBnd,CEstream,sngPnt]

Deterministic orbit (Euler method)

v0={2,4};t0=0;t∞=60000;Δt=10;data={};
v=v0;t=t0;
While[t≤t∞&&n[v]>0,data=Join[data,{Append[v,t]}];v=Δt (Subscript[drift, mo])[v]+v;t=t+Δt;];

... projected in the (x,y)-plane

CEorbit=ListPlot[data[[All,{1,2}]],PlotStyle->{Gray,Thickness[.02]},Joined->True];

Show[coexBnd,CEorbit,preyES,preyNES,predES,predNES,CEstream,sngPnt]

CEorbit=ListPlot[data[[All,{1,2}]],PlotStyle->{Gray,Thickness[.02]},Joined->True];

Show[coexBnd,CEorbit,preyES,preyNES,predES,predNES,CEstream,sngPnt,PlotRange->{{-.1,.1},{.2,1}}]

... as evolutionary tree

xtree=ListPlot[data[[All,{1,3}]],PlotStyle->{Darker[Green],Thick},Joined->True];
ytree=ListPlot[data[[All,{2,3}]],PlotStyle->{Red,Thick},Joined->True];
pLevel=Graphics[{Dashed,Line[{{-p,0},{-p,60000}}]}];

Show[xtree,ytree,pLevel,AxesLabel->{"log body size","evolutionary time"}]

Evolutionary branching in the prey?

Continuationwithlastdataset;
Noticethattheprey-isoclineneartheintersectionisnotevolutionarilystable(dashedisocline).Doesthatmeanthatthepreystrategywillbranch?

PIP for predator trait fixed at its singular value

PIPbnd=ContourPlot[If[n[{x,-p}]>0,(Subscript[sPrey, mo])[{x,-p},X]],{x,xmin,xmax},{X,xmin,xmax},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->60];
nPos=ContourPlot[n[{x,-p}],{x,xmin,xmax},{X,xmin,xmax},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->60];
PIPint=DensityPlot[If[n[{x,-p}]>0&&(Subscript[sPrey, mo])[{x,-p},X]>0,(Subscript[sPrey, mo])[{x,-p},X]],{x,xmin,xmax},{X,xmin,xmax},ColorFunction->"AvocadoColors",PlotPoints->60];
Show[PIPint,PIPbnd,nPos]

MIP near singular point

MIPbnd=ContourPlot[If[Abs[x-X]>0.001,(Subscript[sPrey, mo])[{x,-p},X] (Subscript[sPrey, mo])[{X,-p},x]],{x,-0.5,0.5},{X,-0.5,0.5},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->60];
MIPint=DensityPlot[If[(Subscript[sPrey, mo])[{x,-p},X]>0∧(Subscript[sPrey, mo])[{X,-p},x]>0,(Subscript[sPrey, mo])[{x,-p},X] (Subscript[sPrey, mo])[{X,-p},x]],{x,-0.5,0.5},{X,-0.5,0.5},ColorFunction->"AvocadoColors",PlotPoints->60];
Show[MIPint,MIPbnd]

Conclusion: evolutionary branching might be possible, but since the cone of mutual invadability of resident and mutant prey near the singularity is so very narrow, the dimorphic orbit might not stay inside it. This we have to investigate further.

DIMORPHIC PREY POPULATION AND MONOMORPHIC PREDATOR POPULATION

Reset

Clear[r,K,a,b,c,d];

Population equations

m1,m2 = prey pop dens; nn = pred pop dens; x1,x2 = prey log body size; y = pred log body size

dLogm1[m1_,m2_,nn_]:=r[x1] (1-((m1 a[x1,x1]+m2 a[x1,x2])/(K[x1])))-nn b[x1,y];
dLogm2[m1_,m2_,nn_]:=r[x2] (1-((m1 a[x2,x1]+m2 a[x2,x2])/(K[x2])))-nn b[x2,y];
dLognn[m1_,m2_,nn_]:=m1 b[x1,y] c[x1,y]+m2 b[x2,y] c[x2,y]-d[y];

Population equilibrium

Solve[{dLogm1[m1,m2,nn]==0,dLogm2[m1,m2,nn]==0,dLognn[m1,m2,nn]==0},{m1,m2,nn}]//Simplify

m1[{x1_,x2_,y_}]:=(-a[x1,x2] b[x2,y] d[y] K[x2] r[x1]+K[x1] (((b[x2,y])^2) c[x2,y] K[x2] r[x1]+a[x2,x2] b[x1,y] d[y] r[x2]-b[x1,y] b[x2,y] c[x2,y] K[x2] r[x2]))/(-a[x1,x2] b[x1,y] b[x2,y] c[x1,y] K[x2] r[x1]+a[x1,x1] ((b[x2,y])^2) c[x2,y] K[x2] r[x1]+b[x1,y] (a[x2,x2] b[x1,y] c[x1,y]-a[x2,x1] b[x2,y] c[x2,y]) K[x1] r[x2]);
m2[{x1_,x2_,y_}]:=(a[x1,x1] b[x2,y] d[y] K[x2] r[x1]+b[x1,y] K[x1] (-b[x2,y] c[x1,y] K[x2] r[x1]+(-a[x2,x1] d[y]+b[x1,y] c[x1,y] K[x2]) r[x2]))/(-a[x1,x2] b[x1,y] b[x2,y] c[x1,y] K[x2] r[x1]+a[x1,x1] ((b[x2,y])^2) c[x2,y] K[x2] r[x1]+b[x1,y] (a[x2,x2] b[x1,y] c[x1,y]-a[x2,x1] b[x2,y] c[x2,y]) K[x1] r[x2]);
nn[{x1_,x2_,y_}]:=(((a[x2,x2] b[x1,y] c[x1,y]-a[x2,x1] b[x2,y] c[x2,y]) K[x1]+a[x1,x2] (a[x2,x1] d[y]-b[x1,y] c[x1,y] K[x2])+a[x1,x1] (-a[x2,x2] d[y]+b[x2,y] c[x2,y] K[x2])) r[x1] r[x2])/(-a[x1,x2] b[x1,y] b[x2,y] c[x1,y] K[x2] r[x1]+a[x1,x1] ((b[x2,y])^2) c[x2,y] K[x2] r[x1]+b[x1,y] (a[x2,x2] b[x1,y] c[x1,y]-a[x2,x1] b[x2,y] c[x2,y]) K[x1] r[x2]);

Parameter values and functions

r[x_]:=1;
K[x_]:=(E^(-(x^2)));
a[x1_,x2_]:=(E^(-α (((x1-x2))^2)));α=0.5;
b[x_,y_]:=β1 (E^(-β2 (((x-y-p))^2)));β1=100;β2=0.2;p=Log[0.5];(* attack rate *)
c[x_,y_]:=γ (E^(x-y));γ=0.2;(* conversion coefficient *)
d[y_]:=δ1 (E^(-δ2 y));δ1=1;δ2=1;

Parameter set

μPrey=1;
σ2Prey=0.1;
μPred=1;
σ2Pred=1;

Invasion fitness, gradient and curvature

Prey

(Subscript[sPrey, di])[{x1_,x2_,y_},xMut_]:=r[xMut] (1-((a[xMut,x1] m1[{x1,x2,y}]+a[xMut,x2] m2[{x1,x2,y}])/(K[xMut])))-b[xMut,y] nn[{x1,x2,y}];

(Subscript[d1sPrey, di])[{x1_,x2_,y_}]:=(∂_xMut)(Subscript[sPrey, di])[{x1,x2,y},xMut]/.{xMut->x1};
(Subscript[d2sPrey, di])[{x1_,x2_,y_}]:=(∂_xMut)(Subscript[sPrey, di])[{x1,x2,y},xMut]/.{xMut->x2};

(Subscript[dd1sPrey, di])[{x1_,x2_,y_}]:=(∂_(xMut,xMut))(Subscript[sPrey, di])[{x1,x2,y},xMut]/.{xMut->x1};
(Subscript[dd2sPrey, di])[{x1_,x2_,y_}]:=(∂_(xMut,xMut))(Subscript[sPrey, di])[{x1,x2,y},xMut]/.{xMut->x2};

Predator invasion fitness and its derivatives:

(Subscript[sPred, di])[{x1_,x2_,y_},yMut_]:=b[x1,yMut] c[x1,yMut] m1[{x1,x2,y}]+b[x2,yMut] c[x2,yMut] m2[{x1,x2,y}]-d[yMut];

(Subscript[dsPred, di])[{x1_,x2_,y_}]:=(∂_yMut)(Subscript[sPred, di])[{x1,x2,y},yMut]/.{yMut->y};

(Subscript[ddsPred, di])[{x1_,x2_,y_}]:=(∂_(yMut,yMut))(Subscript[sPred, di])[{x1,x2,y},yMut]/.{yMut->y};

Canonical equation

Deterministic drift

(Subscript[drift, di])[{x1_,x2_,y_}]:={(1/2) μPrey σ2Prey m1[{x1,x2,y}] (Subscript[d1sPrey, di])[{x1,x2,y}],(1/2) μPrey σ2Prey m2[{x1,x2,y}] (Subscript[d2sPrey, di])[{x1,x2,y}],(1/2) μPred σ2Pred nn[{x1,x2,y}] (Subscript[dsPred, di])[{x1,x2,y}]};

Stream plot (note below that the boundaries of the mutual invadability cone are attracting, and so the cone is not forward invarant under the CE)

CEstream=StreamPlot[If[m1[{x1,x2,-p}]>0∧m2[{x1,x2,-p}]>0∧Abs[x1-x2]>0.01,(Subscript[drift, di])[{x1,x2,-p}][[1;;2]],{0,0}],{x1,-0.5,0.5},{x2,-0.5,0.5}]//Quiet;

Show[MIPint,MIPbnd,CEstream]

Deterministic evolutionary tree

(* monomorphic prey *)
v0={2,4};t0=0;t∞=100000;Δt=10;
x1data={};x2data={};ydata={};x1stab={};x2stab={};ystab={};
v=v0;t=t0;
While[t≤t∞∧n[v]>0∧Abs[v[[1]]]>(10^(-6))∧Abs[v[[2]]+p]>(10^(-6)),
x1data=Join[x1data,{{v[[1]],t}}];
x2data=Join[x2data,{{v[[1]],t}}];
ydata=Join[ydata,{{v[[2]],t}}];
x1stab=Join[x1stab,{{(Subscript[ddsPrey, mo])[v],t}}];
x2stab=Join[x2stab,{{(Subscript[ddsPrey, mo])[v],t}}];
ystab=Join[ystab,{{(Subscript[ddsPred, mo])[v],t}}];
v=v+Δt (Subscript[drift, mo])[v];
t=t+Δt];

(* dimorphic prey *)
v={v[[1]]-0.01 (Sqrt[σ2Prey]),v[[1]]+0.01 (Sqrt[σ2Prey]),v[[2]]};

While[t≤t∞&&m1[v]>0&&m2[v]>0&&nn[v]>0,
x1data=Join[x1data,{{v[[1]],t}}];
x2data=Join[x2data,{{v[[2]],t}}];
ydata=Join[ydata,{{v[[3]],t}}];
x1stab=Join[x1stab,{{(Subscript[dd1sPrey, di])[v],t}}];
x2stab=Join[x2stab,{{(Subscript[dd2sPrey, di])[v],t}}];
ystab=Join[ystab,{{(Subscript[ddsPred, di])[v],t}}];
v=Δt (Subscript[drift, di])[v]+v;
t=t+Δt];

Trait values

x1tree=ListPlot[x1data,PlotStyle->{Darker[Green],Thick},Joined->True];
x2tree=ListPlot[x2data,PlotStyle->{Darker[Green],Thick},Joined->True];
ytree=ListPlot[ydata,PlotStyle->{Red,Thick},Joined->True];
Show[x1tree,x2tree,ytree,AxesLabel->{"log body size","evolutionary time"}]

Uninvadability test

x1stree=ListPlot[x1stab,PlotStyle->{Black,Thick},Joined->True];
x2stree=ListPlot[x2stab,PlotStyle->{Black,Thick},Joined->True];
ystree=ListPlot[ystab,PlotStyle->{Red,Thick},Joined->True];

Show[x1stree,x2stree,ystree,AxesLabel->{"fitness curvature","evolutionary time"}]

Conclusion: after branching in the prey species, the population evoves to an evolutionarily stable point with two prey types and one predator type. However, as the mutual invadability cone is not forward invariant (see above) stochastic orbits may readily leave the cone so that branching may fail. For sufficiently small mutation step sizes, however, there is a positive probability of not leaving the cone, an so there is a positive probability of branching. Since after failed branching the population will in the neighborhood of the singular point, branching will occur sooner or later with probability one. The time till successful branching should be exponentially distributed and may be very long. This phenomenon has been dubbed "evolutionary loitering" (pc. Eva Kisdi). The trimorphism (2 prey + 1 pred) reached in the end is uninvadable, because the curvature of the fitness landscape is negative for both prey and for the predator (see last figure above).

Stability of trimorphism

Singular point

x1Sing=x1data[[-1,1]];
x2Sing=x2data[[-1,1]];
ySing=ydata[[-1,1]];

Jacobi matrix {{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}}

A=D[{(Subscript[d1sPrey, di])[{x1,x2,y}],(Subscript[d2sPrey, di])[{x1,x2,y}],(Subscript[dsPred, di])[{x1,x2,y}]},{{x1,x2,y}}]/.{x1->x1Sing,x2->x2Sing,y->ySing};

Total stability test

d1=1;(* just a matter of scaling *)

RegionPlot[ d1 A[[1,1]]+d2 A[[1,2]]+d3 A[[1,3]]<0∧d1 A[[2,1]]+d2 A[[2,2]]+d3 A[[2,3]]<0∧d1 A[[3,1]]+d2 A[[3,2]]+d3 A[[3,3]]<0,{d2,0,12},{d3,0,2},PlotPoints->100,FrameLabel->{"(Subscript[d, 2])⟶","(Subscript[d, 3])⟶"}]

Conclusion: The trimorphic singular point is totally stable (and hence also strongly and weakly stable);