Evo-eco dynamics

Example: LV symmentric competition model

Model definition

Partitioning strategy space

xMin=-1(*min trait value*);
xMax=1(*max trait value*);
iMax=50 (*number of bins*);
dx=N[(xMax-xMin)/iMax] (*bin width*); 
Do[(Subscript[x, i])=xMin+i dx,{i,0,iMax}](*discretized trait values*);

Functions and parameters

Do[(Subscript[r, i])=.5Exp[-((Subscript[x, i])^2)],{i,0,iMax}];
Do[(Subscript[K, i])=Exp[-((Subscript[x, i])^2)],{i,0,iMax}];
Do[(Subscript[a, i,j])=Exp[-2.((((Subscript[x, i])-(Subscript[x, j])))^2)],{i,0,iMax},{j,0,iMax}];

Population model

Do[(Subscript[dn, i])[t_]=(Subscript[r, i])(Subscript[n, i])[t](1-(((∑_(j=1)^iMax)(Subscript[a, i,j])(Subscript[n, j])[t])/(Subscript[K, i]))),{i,0,iMax}];

Equations and solve variables for NDSolve

eqs=Flatten[Table[{(Subscript[n, i])'[t]==(Subscript[dn, i])[t],(Subscript[n, i])[t0]==(Subscript[n0, i])},{i,1,iMax}]];
var=Table[(Subscript[n, i]),{i,1,iMax}];

Trait substitutions (ecological time scale)

Simulation parameters

xIni=xMin+2dx(*initial trait value*);
nIni=nMin(*initial population density of initial trait value*);
nMin=.001(*extinction threshold = initial mututant density*);
mpb=.05(*mutation probability per unit of time and per individual *);
nMax=1.(*max dens -- affects plots only*);
dt=50(*update time for pruning and mutating*);
tMax=2000(*max simulation time*);

Initialize and run

(*initialize*)
Do[(Subscript[n0, i])=0,{i,0,iMax}](*all population densities set to zero*);
iIni=Round[(xIni-xMin)/dx](*except for the initial trait vallue xIni*);
(Subscript[n0, iIni])=nIni;
t0=0;
dat={}(*plot output list*);

(*run*)
While[t0≤tMax,
sol=NDSolve[eqs,var,{t,t0,t0+dt}](*solve*);
dat=Join[dat,{Plot[Evaluate[Table[(Subscript[n, i])[t],{i,0,iMax}]/.sol],{t,t0,t0+dt},PlotRange->{{0,tMax},{0,nMax}},Filling->Bottom]}](*store plot*);
Do[(Subscript[n0, i])=Chop[Evaluate[(Subscript[n, i])[t0+dt]/.sol[[1]]],nMin],{i,0,iMax}](*prune*);
Do[(Subscript[dn0, i])=0,{i,0,iMax}](*reset change*);
Do[(*mutate? loop*)
If[RandomReal[]<mpb (Subscript[n0, i])dt (*yep, mutate*),
If[RandomReal[]>.5(*up or down?*),
{(Subscript[dn0, i-1]),(Subscript[dn0, i])}={(Subscript[dn0, i-1])+nMin,(Subscript[dn0, i])-nMin}(*down*),
{(Subscript[dn0, i+1]),(Subscript[dn0, i])}={(Subscript[dn0, i+1])+nMin,(Subscript[dn0, i])-nMin}(*up*)
]
],
{i,1,iMax-1}(*note that the boundaries of the strategy space are absorbing*)
](*end mutate? loop*);
Do[(Subscript[n0, i])=(Subscript[n0, i])+(Subscript[dn0, i]),{i,0,iMax}](*update population densities*);
t0=t0+dt(*update time*);
];

Population density vs time

Show[dat,AxesLabel->{"t","(Subscript[n, i])"}](*N.B. the small-scale dips in curves are due to mutants being added or taken away*)

Evolutionary tree (evolutionary time scale)

Simulation parameters

xIni=xMin+2dx(*initial trait value*);
nMin=.005(*extinction threshold = initial mututant density*);
mpb=.05(*mutation probability per unit of time and per individual *);
dt=10(*update time for pruning and mutating*);
Δt=50(*output time, i.e., updating the output*);
tMax=20000(*max simulation time*);

Initialize and run

(*initialize*)
Do[(Subscript[n0, i])=0,{i,0,iMax}];
iIni=Round[(xIni-xMin)/dx];
(Subscript[n0, iIni])=(Subscript[K, iIni]);
t0=0;
Clear[dat];
dat={{t0,(Subscript[x, iIni-1]),0},{t0,(Subscript[x, iIni]),(Subscript[n0, iIni])},{t0,(Subscript[x, iIni+1]),0}};

(*counters*)
cntMax=Round[Δt/dt](*update loop*);
CntMax=Round[tMax/Δt](*output loop*);

(*run*)
Do[(*output loop*)
Do[(*update loop*)
sol=NDSolve[eqs,var,{t,t0,t0+dt}](*solve*);
Do[(Subscript[n0, i])=Chop[Evaluate[(Subscript[n, i])[t0+dt]/.sol[[1]]],nMin],{i,0,iMax}](*prune*);
Do[(Subscript[dn0, i])=0,{i,0,iMax}](*reset change*);
Do[(*mutate? loop*)
If[RandomReal[]<mpb (Subscript[n0, i])dt (*mutate*),
If[RandomReal[]>.5(*up or down?*),
{(Subscript[dn0, i-1]),(Subscript[dn0, i])}={(Subscript[dn0, i-1])+nMin,(Subscript[dn0, i])-nMin}(*down*),
{(Subscript[dn0, i+1]),(Subscript[dn0, i])}={(Subscript[dn0, i+1])+nMin,(Subscript[dn0, i])-nMin}(*up*)
]
],
{i,1,iMax-1}(*note that the boundaries of the strategy space are absorbing*)
](*end mutate? loop*);
Do[(Subscript[n0, i])=(Subscript[n0, i])+(Subscript[dn0, i]),{i,0,iMax}](*update population densities*);
t0=t0+dt(*update time*),
{cnt,1,cntMax}
](*end update loop*);
Do[If[(Subscript[n0, i-1])≥nMin∨(Subscript[n0, i])≥nMin∨(Subscript[n0, i+1])≥nMin,dat=Join[dat,{{t0,(Subscript[x, i]),(Subscript[n0, i])}}]],{i,1,iMax-1}](*store data*),
{Cnt,1,CntMax}
](*end output loop*);

Length[dat]

Show[
{DensityPlot[0,{t,0,tMax},{x,xMin,xMax},ColorFunction->"AvocadoColors",AspectRatio->.7],

ListDensityPlot[dat,InterpolationOrder->0,PlotRange->All,ColorFunction->"AvocadoColors"]},
PlotLabel->"Evolutionary dynamics"
]

Show[
Plot3D[nMin,{t,0,tMax},{x,xMin,xMax},PlotRange->{nMin,.8},PlotStyle->Black],

ListPlot3D[dat,InterpolationOrder->0,PlotRange->{0,.8},Mesh->None,BoundaryStyle->None,PlotStyle->Green,Filling->Bottom,FillingStyle->Darker[Green]],

AxesLabel->{"t","x","n"},PlotLabel->"Eco-evo dynamics"
]

Show[
Plot3D[0,{t,0,tMax},{x,xMin,xMax},PlotRange->All,ColorFunction->"AvocadoColors"],

ListPlot3D[dat,InterpolationOrder->1,PlotRange->All,ColorFunction->"AvocadoColors",
AxesLabel->{"t","x","n"}]
]

Show[
Plot3D[0,{t,0,tMax},{x,xMin,xMax},PlotRange->All,ColorFunction->"AvocadoColors"],

ListPlot3D[dat,InterpolationOrder->1,PlotRange->{0,10nMin},Mesh->None,ColorFunction->"AvocadoColors",
AxesLabel->{"t","x","n"}]
]

Show[
Plot3D[0,{t,0,tMax},{x,xMin,xMax},PlotRange->All],

ListPlot3D[dat,InterpolationOrder->1,PlotRange->All]
]

Show[
Plot3D[0,{t,0,tMax},{x,xMin,xMax},PlotRange->All],

ListPlot3D[dat,InterpolationOrder->1,PlotRange->All,Mesh->None]
]