Competition colonization trade-off

Model description

Population dynamics

((Subscript[dp, i])/dt)=(Subscript[β, i]) (Subscript[p, i])((∑_(∀j))(Subscript[γ, j,i])(Subscript[p, j]))-(Subscript[p, i])((∑_(∀j))(Subscript[γ, i,j])(Subscript[β, j])(Subscript[p, j]))-(Subscript[δ, i])(Subscript[p, i]);
                        ⟵---⟶            ⟵-----⟶
per offspring          rate of   
                     probability of        displacement
                     colonization           from a site
with

(Subscript[p, i]) = proportion of sites occupied by 'species' i

and

(Subscript[β, i]) = per capita birth rate for type i
(Subscript[δ, i]) = per capita death rate for type i
(Subscript[γ, i,j])= prob. of site owner of type i being displaced by intruder of type j

Rewrite as

((Subscript[dp, i])/dt)=(Subscript[r, i])(Subscript[p, i])(1-(((∑_(∀j))(Subscript[a, i,j])(Subscript[p, j]))/(Subscript[K, i])))

with

(Subscript[r, i])=(Subscript[β, i])-(Subscript[δ, i])

(Subscript[K, i])=1-((Subscript[δ, i])/(Subscript[β, i]))
(Subscript[a, i,j])=1-(Subscript[γ, j,i])+(Subscript[γ, i,j])((Subscript[β, j])/(Subscript[β, i])) 

Invasion fitness

(Subscript[s, x])[y]=r[y](1-((a[y,x]K[x])/(K[y])))

(Subscript[s, (Subscript[x, 1]),(Subscript[x, 2])])[y]=r[y](1-((a[y,(Subscript[x, 1])](Subscript[n, 1])[(Subscript[x, 1]),(Subscript[x, 2])]+a[y,(Subscript[x, 2])](Subscript[n, 2])[(Subscript[x, 1]),(Subscript[x, 2])])/(K[y])))

with

(Subscript[n, 1])[(Subscript[x, 1]),(Subscript[x, 2])]=((K[(Subscript[x, 1])]-a[(Subscript[x, 1]),(Subscript[x, 2])]K[(Subscript[x, 2])])/(1-a[(Subscript[x, 1]),(Subscript[x, 2])]a[(Subscript[x, 2]),(Subscript[x, 1])]))

(Subscript[n, 2])[(Subscript[x, 1]),(Subscript[x, 2])]=((K[(Subscript[x, 2])]-a[(Subscript[x, 2]),(Subscript[x, 1])]K[(Subscript[x, 1])])/(1-a[(Subscript[x, 1]),(Subscript[x, 2])]a[(Subscript[x, 2]),(Subscript[x, 1])]))

Implementation

Ecology

r[x_]:=β[x]-δ[x];

K[x_]:=1-((δ[x])/(β[x]));
a[x_,y_]:=1-γ[y,x]+γ[x,y]((β[y])/(β[x]));

Invasion fitness (monomorphic) and fitness gradient & curvatures

(Subscript[s, x_])[y_]:=r[y](1-((a[y,x]K[x])/(K[y]))); 

grad[x_]:=(∂_η)(Subscript[s, ξ])[η]/.{ξ->x,η->x};
xCurv[x_]:=(∂_(ξ,ξ))(Subscript[s, ξ])[η]/.{ξ->x,η->x};
yCurv[x_]:=(∂_(η,η))(Subscript[s, ξ])[η]/.{ξ->x,η->x};

Invasion fitness & population densities (dimorphic) and fitness gradient & curvatures

(Subscript[s, x1_,x2_])[y_]:=r[y](1-((a[y,x1]n1[x1,x2]+a[y,x2]n2[x1,x2])/(K[y])));

n1[x1_,x2_]:=((K[x1]-a[x1,x2]K[x2])/(1-a[x1,x2]a[x2,x1]));

n2[x1_,x2_]:=((K[x2]-a[x2,x1]K[x1])/(1-a[x1,x2]a[x2,x1]));

grad1[x1_,x2_]:=(∂_η)(Subscript[s, ξ1,ξ2])[η]/.{ξ1->x1,ξ2->x2,η->x1};
grad2[x1_,x2_]:=(∂_η)(Subscript[s, ξ1,ξ2])[η]/.{ξ1->x1,ξ2->x2,η->x2};
curv1[x1_,x2_]:=(∂_(η,η))(Subscript[s, ξ1,ξ2])[η]/.{ξ1->x1,ξ2->x2,η->x1};
curv2[x1_,x2_]:=(∂_(η,η))(Subscript[s, ξ1,ξ2])[η]/.{ξ1->x1,ξ2->x2,η->x2};

First run

Case: γ[(Subscript[x, i]),(Subscript[x, j])]=1-γ[(Subscript[x, j]),(Subscript[x, i])]  (i.e., no effect of ownership)

Functions and parameters

β[x_]:=1/x;
δ[x_]:=1;
γ[x_,y_]:=c[y]/(c[x]+c[y]);
c[x_]:=(E^(α x)); (* α is a measure for the asymmetry of competition *)

Strategy range

xMin=0;
xMax=1;

Pairwise invadability plots (PIP)

Do[
PIP[α]=
RegionPlot[(Subscript[s, x])[y]>0,{x,0,1},{y,0,1}],
{α,{5,5.82,7,8.4,9.3,10}}
];

Row[
Table[
Show[PIP[α],PlotLabel->α,ImageSize->Small],
{α,{5,5.82,7,8.4,9.3,10}}]
]

Bifurcation Plot (x vs α)

ContourPlot[grad[x],{α,0,20},{x,xMin,xMax},Contours->{0},ContourStyle->{Black,Thick},ContourShading->False,PlotPoints->30,ImageSize->Small]

Show[
(* ES attractor *)
ContourPlot[If[yCurv[x]<Min[xCurv[x],0],grad[x]],{α,0,20},{x,xMin,xMax},Contours->{0},ContourStyle->{Black,Thickness[0.02]},ContourShading->False,PlotPoints->30],
(* ES repeller *)
ContourPlot[If[xCurv[x]<yCurv[x]<0,grad[x]],{α,0,20},{x,xMin,xMax},Contours->{0},ContourStyle->{Black,Thickness[.005]},ContourShading->False,PlotPoints->30],
(* non-ES repeller *)
ContourPlot[If[Max[xCurv[x],0]<yCurv[x],grad[x]],{α,0,20},{x,xMin,xMax},Contours->{0},ContourStyle->{Red,Thickness[.005]},ContourShading->False,PlotPoints->30],
(* BP *)
ContourPlot[If[0<yCurv[x]<xCurv[x],grad[x]],{α,0,20},{x,xMin,xMax},Contours->{0},ContourStyle->{Red,Thickness[0.02]},ContourShading->False,PlotPoints->30],
ImageSize->Small
]//Quiet

Mutual invadability plots (MIP)

Do[
MIP[α]=
RegionPlot[(Subscript[s, x1])[x2]>0 ∧ (Subscript[s, x2])[x1]>0,{x1,0,1},{x2,0,1}],
{α,{5,5.82,7,8.4,9.3,10}}
]//Quiet;

Row[
Table[
Show[MIP[α],PlotLabel->α,ImageSize->Small],
{α,{5,5.82,7,8.4,9.3,10}}
]
]

"Anti-MIP" (to be used later as a mask to cover what happens outside the oexistence region)

Do[
antiMIP[α]=
RegionPlot[
(Subscript[s, x])[y]≤0 ∨ (Subscript[s, y])[x]≤0,{x,xMin,xMax},{y,xMin,xMax},
PlotPoints->100
],
{α,{5,5.82,7,8.4,9.3,10}}
]

Evolutionary isoclines

(* SLOW PROCEDURE - do not run in class *)

Do[
IPes1[α]=ContourPlot[If[(0<n1[x1,x2])∧(0<n2[x1,x2])∧0>curv1[x1,x2]∧(Abs[x1-x2]>.01(xMax-xMin)),grad1[x1,x2]],{x1,xMin,xMax},{x2,xMin,xMax},Contours->{0},ContourStyle->{Darker[Green],Thick},ContourShading->False,PlotPoints->10];
IPbp1[α]=ContourPlot[If[(0<n1[x1,x2])∧(0<n2[x1,x2])∧0<curv1[x1,x2]∧(Abs[x1-x2]>.01(xMax-xMin)),grad1[x1,x2]],{x1,xMin,xMax},{x2,xMin,xMax},Contours->{0},ContourStyle->{Darker[Red],Thick},ContourShading->False,PlotPoints->10];
IPes2[α]=ContourPlot[If[(0<n1[x1,x2])∧(0<n2[x1,x2])∧0>curv2[x1,x2]∧(Abs[x1-x2]>.01(xMax-xMin)),grad2[x1,x2]],{x1,xMin,xMax},{x2,xMin,xMax},Contours->{0},ContourStyle->{Darker[Green],Dashed,Thick},ContourShading->False,PlotPoints->10];
IPbp2[α]=ContourPlot[If[(0<n1[x1,x2])∧(0<n2[x1,x2])∧0<curv2[x1,x2]∧(Abs[x1-x2]>.01(xMax-xMin)),grad2[x1,x2]],{x1,xMin,xMax},{x2,xMin,xMax},Contours->{0},ContourStyle->{Darker[Red],Dashed,Thick},ContourShading->False,PlotPoints->10],
{α,{5,5.82,7,8.4,9.3,10}}
]//Quiet;

Row[
Table[
Show[IPes1[α],IPbp1[α],IPes2[α],IPbp2[α],antiMIP[α],
PlotLabel->α,ImageSize->Small],
{α,{5,5.82,7,8.4,9.3,10}}
]
]

Canonical equation (CE)

drift[{x1_,x2_}]={.5μ σ2[x1] n1[x1,x2]grad1[x1,x2],.5μ σ2[x2] n2[x1,x2]grad2[x1,x2]};

μ=1;
σ2[x_]:=.001(Sqrt[x(1-x)]);

Do[
stream[α]=StreamPlot[drift[{x1,x2}],{x1,xMin,xMax},{x2,xMin,xMax},StreamPoints->Medium,StreamScale->Large,StreamStyle->LightGray],
{α,{5,5.82,7,8.4,9.3,10}}
]//Quiet;

Row[
Table[
Show[stream[α],IPes1[α],IPbp1[α],IPes2[α],IPbp2[α],antiMIP[α],
PlotLabel->α,ImageSize->Small],
{α,{5,5.82,7,8.4,9.3,10}}
]
]