Simulate adaptation on a fitness landscape

Author: Alexandre Blanckaert, Claudia Bank
Version: 17 April 2018

Instruction: Evaluate the notebook first by choosing “Evaluation->Evaluate Notebook” in the Menu. To the system dialog asking whether you want to evaluate initialization cells, you reply “yes”. You can then try out different configurations of the model at the end of this notebook.

Simulation code (can be ignored)

Fisher[x0_,y0_,σ_]:=Exp[-((x-x0)^2+((y-y0)^2))/σ] (* Fisher geometric model; creates a peak around {x0, y0}, width of the pic is controlled by σ *)

OneStep[{x0_,y0_},var_,func_,Nind_]:=Module[{x1,y1,t}, (* create a random step and decide whether it is selected for or not*)
x1=RandomVariate[NormalDistribution[0,var]]; (*change in the x axis*)
y1=RandomVariate[NormalDistribution[0,var]]; (*change in the y axis*)
t=(func/.{x->(x0+x1),y->(y0+y1)})-(func/.{x->x0,y->y0}); (*calculate the fitness of the mutant*)
If[(1-Exp[-2*t])/(1-Exp[-2*Nind*t])>RandomVariate[UniformDistribution[]],Return[{x0+x1,y0+y1}],Return[{x0,y0}]]] (* see if the mutation escapes drift, using Kimura formula for fixation proba*)

OneTraj[{x0_,y0_},var_,func_,Nind_,nPoints_,nrej_]:=Module[{traj,u,reject},traj={{x0,y0}};reject=0; (* calculate a single trajectories along the fitness landscape*)
While[Length[traj]<nPoints&&reject<nrej, (* we lookk for the next until the trajectories includes nPoints, or we made nrej unsuccessful attempt (to avoid forever calculation) *)
u=OneStep[traj[[Length[traj]]],var,func,Nind];
reject=reject+1;
If[u!=traj[[Length[traj]]],reject=0;AppendTo[traj,u]];
];
Return[traj];
]

MakeHotPlates[nbPic_,mutsize_,nrep_,Nind_,origin_]:=Module[{h,zz,trajzz,veccol,var,c0,nPoints,nrej,start,start3D,pnt},
Which[mutsize=="Very small",var=.005,mutsize=="Small",var=0.02,mutsize=="Large",var=.05]; (*choose mutation size*)
Which[nbPic==1,h=Fisher[-.3,-.3,10]+Fisher[.75,-.75,.2],nbPic==2,h=Fisher[-.3,-.3,10]+Fisher[.42,.42,.2]*.8+Fisher[.75,-.75,.2],nbPic==3,h=Fisher[-.3,-.3,10]+Fisher[.42,.42,.2]*.8+Fisher[-.2,-.8,.1]+Fisher[.75,-.75,.2],nbPic==4,h=Fisher[-.3,-.3,10]+Fisher[.42,.42,.2]*.8+Fisher[-.2,-.8,.1]+Fisher[.75,-.75,.2]+Fisher[-.8,0,.05]*.9,nbPic==5,h=Fisher[-.3,-.3,10]+Fisher[.42,.42,.2]*.8+Fisher[-.2,-.8,.1]+Fisher[.75,-.75,.2]+Fisher[-.8,0,.05]*.9+Fisher[-.5,.5,.01]]; (*define the number of peaks*)
Which[origin=="Pnt 1",c0={0,-.2},origin=="Pnt 2",c0={-.318,0},origin=="Pnt 3",c0={-.9,-.5},origin=="Pnt 4",c0={-.8,.8}]; (* define potential starting point *)
If[nrep==0,Which[origin=="Pnt 1",pnt=1,origin=="Pnt 2",pnt=2,origin=="Pnt 3",pnt=3,origin=="Pnt 4",pnt=4]]; (* to highlight the possible staring position*)
start={{0,-.2},{-.318,0},{-.9,-.5},{-.8,.8}};
start3D=Transpose[Flatten[{Transpose[start],{(h/.{x->start[[;;,1]],y->start[[;;,2]]})+.01}},1]];
Which[Nind> 100||var>.1, nPoints=100,Nind>10||var>.01,nPoints=200,True,nPoints=1000]; (* chose number of points depending on parameters of the simulations*)
Which[Nind< 100||var>.1, nrej=1000,Nind<10000||var>.01,nrej=5000,True,nrej=10000]; (* define how hard we look for the next step depending on paremeters of the simulations*)
zz=Table[OneTraj[c0,var,h,Nind,nPoints,nrej],{nrep}]; (* calculate all trajectories*)
trajzz=Table[Transpose[Flatten[{Transpose[zz[[k]]],{(h/.{x->zz[[k]][[;;,1]],y->zz[[k]][[;;,2]]})+.01}},1]],{k,1,nrep}]; (* calculate fitness associated to the trajectories*)
veccol={Blue,Cyan,Green,Red,Magenta,Pink,Purple,Gray,Black,White};
Show[Plot3D[h,{x,-1,1},{y,-1,1},PlotRange->All,ImageSize->600,ColorFunction->"FallColors",Boxed->False,Axes->{False, False, False}],If[nrep==0,Show[ListPointPlot3D[Tooltip[{start3D[[1]]},"Pnt 1"],PlotStyle->{{If[pnt==1,Red,Black],PointSize[.02]}}],ListPointPlot3D[Tooltip[{start3D[[2]]},"Pnt 2"],PlotStyle->{{If[pnt==2,Red,Black],PointSize[.02]}}],ListPointPlot3D[Tooltip[{start3D[[3]]},"Pnt 3"],PlotStyle->{{If[pnt==3,Red,Black],PointSize[.02]}}],ListPointPlot3D[Tooltip[{start3D[[4]]},"Pnt 4"],PlotStyle->{{If[pnt==4,Red,Black],PointSize[.02]}}]],Table[ListPointPlot3D[trajzz[[k]],PlotStyle->{{veccol[[k]],PointSize[Large]}}],{k,1,nrep}]] (* making the plot*)
]]

Interactive visualization of adaptive walks

Some tips for getting started
0) You change the parameters of the model by directly clicking on the choices displayed below.
1) Initially, the “Number of replicates” is set to “0”, thus, no adaptive walk is present. Get the simulation started by setting it to >0.
2) “Position of the population” is indicated in red in the landscape. It is the point at which the adaptive walks will begin.
3) The time to compute a new scenario increases with the population size and with the number of replicates.
4) “Size of mutation” refers to the mean distance a single mutation can traverse in the fitness landscape.
5) Try rotating the fitness landscape plot and looking at the walks from different perspectives.

Manipulate[MakeHotPlates[nPic,var,nrep,Nind,origin],{{nPic,3,"Number of fitness peaks"},{1,2,3,4,5}},{{var,"Small","Size of mutation"},{"Very small","Small","Large"}},{{nrep,0,"Number of replicates"},{0,1,2,3,4,5,6,7,8,9,10}},{{Nind,100,"Population size"},{10,100,500,1000,10000,1000000}},{{origin,"Pnt 1","Position of the population"},{"Pnt 1","Pnt 2","Pnt 3","Pnt 4"}},ControlType->Setter,ContinuousAction->False]