The Lotka - Volterra cannibalism time budget model

(d/(dt))R=r R (1-(R/K))-α R(∑_(j=1)^k)(1-(Subscript[x, j]))(Subscript[n, j])

(d/(dt))(Subscript[n, i])=ϵ α R (1-(Subscript[x, i]))(Subscript[n, i])-δ (1-(Subscript[x, i]))(Subscript[n, i])-((∑_(j=1)^k)β[(Subscript[x, j])](Subscript[x, j])(Subscript[n, j]))(1-(Subscript[x, i]))(Subscript[n, i])+

   γ β[(Subscript[x, i])](Subscript[x, i])(Subscript[n, i])(∑_(j=1)^k)(1-(Subscript[x, j]))(Subscript[n, j])-δ (Subscript[x, i])(Subscript[n, i])       for i=1,...,k

MONOMORPHIC RESIDENT POPULATION

Clear[α,β,β0,β1,γ,δ,ϵ,r,K,n,R]

Monomorphic resident population dynamics:

(R = resource density; n = resident population density)

dLogR=r-r R/K-α(1-x)n;
dLogn=ϵ α(1-x)R-δ+γ β[x]x(1-x)n-(1-x)β[x]x n;

Monomorphic resident population equilibrium:

Solve[{0==dLogR,0==dLogn},{R,n}]

(*copy and paste from above*)

R[x_]:=-((K (α δ-r x (-1+γ) β[x]))/(K (-1+x) (α^2) ϵ+r x (-1+γ) β[x]));

n[x_]:=-((r (δ+K (-1+x) α ϵ))/((-1+x) (K (-1+x) (α^2) ϵ+r x (-1+γ) β[x])));

Invasion fitness and its derivatives:

(Subscript[s, x_])[y_]:=ϵ α(1-y)R[x]-δ+γ β[y]y(1-x)n[x]-(1-y)β[x]x n[x];

ds[x_]:=((∂_y)(Subscript[s, x])[y])/.{y->x}; 

dds[x_]:=((∂_y)(∂_y)(Subscript[s, x])[y])/.{y->x};

Default parameter values and functions:

α=1;γ=0.2;δ=0.1;ϵ=0.05;r=1;K=10;

β0=0.;β1=1.5 ;p=.5;

β[x_]:=β0+β1 (x^p);

Plot[β[x],{x,0,1},AxesLabel->{"x","β[x]"},ImageSize->Medium]

Plots of R[x] and n[x]:

Plot[If[n[x]≥0,{R[x],n[x]}],{x,0,1},AxesLabel->{"x","n[x],R[x]"},ImageSize->Medium]

Pairwise invadability plot (PIP):

PIPint=RegionPlot[n[x]≥0&&(Subscript[s, x])[y]>0,{x,0,1},{y,0,1},PlotPoints->100];

nPos=ContourPlot[n[x],{x,0,1},{y,0,1},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->10];

Show[PIPint,nPos,ImageSize->Medium]

Bifurcation plot singular strategy x versus parameter α:

posEQ=RegionPlot[n[x]≥0,{α,0,2},{x,0,1}]; 

attESS=ContourPlot[If[n[x]>0&&ds'[x]<0&&dds[x]<0,ds[x]],{α,0,2},{x,0,1},Contours->{0},ContourStyle->{Green,Thick},ContourShading->False,PlotPoints->50];

repESS=ContourPlot[If[n[x]>0&&ds'[x]>0&&dds[x]<0,ds[x]],{α,0,2},{x,0,1},Contours->{0},ContourStyle->{Green,Thick,Dashed},ContourShading->False,PlotPoints->50];

attNESS=ContourPlot[If[n[x]>0&&ds'[x]<0&&dds[x]>0,ds[x]],{α,0,2},{x,0,1},Contours->{0},ContourStyle->{Red,Thick},ContourShading->False,PlotPoints->50];

repNESS=ContourPlot[If[n[x]>0&&ds'[x]>0&&dds[x]>0,ds[x]],{α,0,2},{x,0,1},Contours->{0},ContourStyle->{Red,Thick,Dashed},ContourShading->False,PlotPoints->50];

posEQbnd=ContourPlot[n[x],{α,0,2},{x,0,1},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->50];

Show[posEQ,posEQbnd,attESS,repESS,attNESS,repNESS,FrameLabel->{"α ⟶","(SuperStar[x]) ⟶"},ImageSize->Medium]

Bifurcation plot singular strategy x versus parameter p:

posEQ=RegionPlot[n[x]≥0,{p,0,2},{x,0,1}]; 

attESS=ContourPlot[If[n[x]>0&&ds'[x]<0&&dds[x]<0,ds[x]],{p,0,2},{x,0,1},Contours->{0},ContourStyle->{Green,Thick},ContourShading->False,PlotPoints->50];

repESS=ContourPlot[If[n[x]>0&&ds'[x]>0&&dds[x]<0,ds[x]],{p,0,2},{x,0,1},Contours->{0},ContourStyle->{Green,Thick,Dashed},ContourShading->False,PlotPoints->50];

attNESS=ContourPlot[If[n[x]>0&&ds'[x]<0&&dds[x]>0,ds[x]],{p,0,2},{x,0,1},Contours->{0},ContourStyle->{Red,Thick},ContourShading->False,PlotPoints->50];

repNESS=ContourPlot[If[n[x]>0&&ds'[x]>0&&dds[x]>0,ds[x]],{p,0,2},{x,0,1},Contours->{0},ContourStyle->{Red,Thick,Dashed},ContourShading->False,PlotPoints->50];

posEQbnd=ContourPlot[n[x],{p,0,2},{x,0,1},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->50];

Show[posEQ,posEQbnd,attESS,repESS,attNESS,repNESS,FrameLabel->{"p ⟶","(SuperStar[x]) ⟶"},ImageSize->Medium]

DIMORPHIC RESIDENT POPULATION

Reset:

Clear[α,β,γ,δ,ϵ,r,K];

Dimorphic resident population equilibrium:

eqn=
{0==r-r R/K-α((1-x1)n1+(1-x2)n2),0==ϵ α(1-x1)R-δ+γ β[x1]x1((1-x1)n1+(1-x2)n2)-(1-x1)(β[x1]x1 n1+β[x2]x2 n2),0==ϵ α(1-x2)R-δ+γ β[x2]x2((1-x1)n1+(1-x2)n2)-(1-x2)(β[x1]x1 n1+β[x2]x2 n2)};
var={R,n1,n2};
Solve[eqn, var]; 
Simplify[%]

R[x1_,x2_]:=
(K ((x1-x2) α δ+r x1 (-1+x2) γ β[x1]-r (-1+x1) x2 γ β[x2]))/(r γ (x1 (-1+x2) β[x1]-(-1+x1) x2 β[x2]));

n1[x1_,x2_]:=(K (x1-x2) (-1+x2) (α^2) δ ϵ+r x1 (-1+x2) γ (δ+K (-1+x2) α ϵ) β[x1]+r x2 (γ (δ-K α ϵ)+x2 (δ-γ δ+K α γ ϵ)-x1 (δ+K (-1+x2) α γ ϵ)) β[x2])/(r γ (((x1 (-1+x2) β[x1]-(-1+x1) x2 β[x2]))^2));

n2[x1_,x2_]:=(r x1 (-x2 δ+γ δ-K α γ ϵ+K x2 α γ ϵ+x1 (δ-γ δ-K (-1+x2) α γ ϵ)) β[x1]-(-1+x1) (K (x1-x2) (α^2) δ ϵ-r x2 γ (δ+K (-1+x1) α ϵ) β[x2]))/(r γ (((x1 (-1+x2) β[x1]-(-1+x1) x2 β[x2]))^2));

Dimorphic invasion fitness and derivatives:

(Subscript[s, x1_,x2_])[y_]:=ϵ α(1-y)R[x1,x2]-δ+γ β[y]y((1-x1)n1[x1,x2]+(1-x2)n2[x1,x2])-(1-y)(β[x1]x1 n1[x1,x2]+β[x2]x2 n2[x1,x2]);

x1ds[x1_,x2_]:=(∂_y)(Subscript[s, x1,x2])[y]/.{y->x1}; 
x2ds[x1_,x2_]:=(∂_y)(Subscript[s, x1,x2])[y]/.{y->x2}; 
x1dds[x1_,x2_]:=(∂_y)(∂_y)(Subscript[s, x1,x2])[y]/.{y->x1};
x2dds[x1_,x2_]:=(∂_y)(∂_y)(Subscript[s, x1,x2])[y]/.{y->x2};

Default parameter values and functions:

α=1;γ=0.2;δ=0.1;ϵ=0.05;r=1;K=10;

β0=0.;β1=1.5 ;p=.5;

β[x_]:=β0+β1 (x^p);

Coexistence plot (MIP):

MIP=RegionPlot[n1[x1,x2]>0∧n2[x1,x2]>0,{x1,0,1},{x2,0,1},PlotPoints->100];

antiMIP=RegionPlot[n1[x1,x2]<0∨n2[x1,x2]<0,{x1,0,1},{x2,0,1},PlotPoints->100,PlotStyle->White];

grad=VectorPlot[If[n1[x1,x2]>0&&n2[x1,x2]>0,{Sign[x1ds[x1,x2]],Sign[x2ds[x1,x2]]},{0,0}],{x1,0,1},{x2,0,1},VectorScale->Small];

Show[MIP,grad,antiMIP]