Continuous-time matching-alleles coevolution

Here we explore host-parasite coevolution between a host and free-living pathogen.  By assuming host and parasite population sizes remain constant over time we are able to focus solely on the effects of coevolution (change in allele frequencies).  Using a Moran model design every host (pathogen) birth is coupled with a host (pathogen) death.  A schematic diagram showing the host-pathogen life cycle is given by figure 2 of the main text.

By running the following line and then compiling this Mathematica notebook will generate the plots found in the main text.

(*CreateDirectory[NotebookDirectory[]<>"Plots"]*)

Deterministic Model

Here we model the dynamics of host and pathogen in the deterministic limit, assuming that host and pathogen population population sizes are infinite.

Differential Equations

We model host-parasite coevolution in a symmetric matching-alleles model as given by the substitution list βSub. 
Although, it is most straightforward to derive the differential equations in terms the number of hosts of type i, nH[i,t], and the number of pathogens of type j, nP[j,t], we transform into the host and pathogen allele frequencies using the substitution pSub.

βSub={β[1,1]->(X/κ),β[1,2]->(Y/κ),β[2,1]->(Y/κ),β[2,2]->(X/κ)};
pSub={nH[1,t]->κ pH[t],nH[2,t]->κ (1-pH[t]),nP[1,t]->κ pP[t],nP[2,t]->κ (1-pP[t])};
Denom=((Sum[β[i,j]α nH[i,t]nP[j,t],{j,1,2},{i,1,2}]+Sum[δ nH[i,t],{i,1,2}])/κ);

Differential equation for host allele frequency (d/dt)(Subscript[p, H])(t)=gH(t)

gH[t_]:=(1/κ)(Sum[((β[i,j]α nH[i,t]nP[j,t])/Denom) ((-Mod[i,2])+((nH[1,t])/κ)),{j,1,2},{i,1,2}]+Sum[((δ nH[i,t])/Denom)((-Mod[i,2])+((nH[1,t])/κ)),{i,1,2}])/.βSub/.pSub//FullSimplify

gH[t]//FullSimplify

Differential equation for pathogen allele frequency (d/dt)(Subscript[p, H])(t)=gH(t)

gP[t_]:=(1/κ)(((Sum[β[i,j]α nH[i,t]nP[j,t]((Mod[j,2] M)-M((nP[1,t])/κ)),{j,1,2},{i,1,2}])/Denom)+Sum[((γ nP[j,t])/Denom) (Mod[j,2]M-M((nP[1,t])/κ)),{j,1,2}])/.M->1/.βSub/.pSub//FullSimplify

gP[t]

Equilibria

equi=Solve[{gH[t]==0,gP[t]==0},{pH[t],pP[t]}]

Stability of the polymorphic equilibrium

jMtrxMoran = {{D[gH[t], pH[t]], D[gH[t], pP[t]]}, {D[gP[t], pH[t]], D[gP[t], pP[t]]}};

Eigenvalues[jMtrxMoran/.equi[[3]]]//FullSimplify

These eigenvalues are purely imaginary (0 real part) and hence the endemic equilibrium is characterized by neutrally stable cycles.  Not that it is straightforward to show that the remaining four equilibria are all unstable.

Numerical Integration

Numerical integration of host and pathogen allele frequency dynamics.

Clear[nsolMoran]
nsolMoran[pars_,inits_]:=nsolMoran[pars,inits]=NDSolve[{D[pH[t],t]==gH[t]/.pars,D[pP[t],t]==gP[t]/.pars,pH[0]==pH0/.inits,pP[0]==pP0/.inits},{pH[t],pP[t]},{t,0,tMax/.pars}]//Flatten
(*Solution for allele frequencies*)
pHsol[τ_,pars_,inits_]:=pH[t]/.nsolMoran[pars,inits]/.{t->τ}
pPsol[τ_,pars_,inits_]:=pP[t]/.nsolMoran[pars,inits]/.{t->τ}

Test parameters and initial conditions

TestPars={X->0.8,Y->0.4,κ->100,δ->1,α->0.2,γ->1,(*M->3,*)tMax->1000};
TestInits0={pH0->0.55,pP0->0.45};
TestInits={pH0->0.85,pP0->0.15};
TestInits2={pH0->0.7,pP0->0.3};

Allele frequency dynamics.  Host=Black Pathogen=Red.

Plot[{pHsol[τ,TestPars,TestInits],pPsol[τ,TestPars,TestInits]},{t,1,500},PlotStyle->{{Black,Thick},{Red,Thick}},Frame->True,FrameTicks->{True,False,True,False},ImageSize->Medium,PlotRange->{0,1},AxesLabel->{"Time","Allele Frequency"}]

Phase plane of allele frequency dynamics.

DetDynAF[pars_,inits_,color_]:=Show[ListLinePlot[Table[{pHsol[τ,pars,inits],pPsol[τ,pars,inits]},{τ,0,1000}],PlotRange->{{0,1},{0,1}},AspectRatio->1,PlotStyle->color,Frame->True,FrameTicksStyle->Directive[Black,Bold,14](*,FrameLabel->{"Host Allele Freq.","Pathogen Allele Freq."}, LabelStyle->Directive[Black,Bold, 18]*)]/.Line[x_]:>{Arrowheads[Table[.05,{r,{2}}]],Arrow[x]},ListPlot[{{pH0,pP0}/.inits},PlotStyle->Red,PlotMarkers->{●,10}],ImageSize->{Automatic,300}]

Show[(*DetDynAF[TestPars,TestInits0,Black],DetDynAF[TestPars,TestInits2,Blue],*)DetDynAF[TestPars,TestInits,Black(*Darker[Green]*)],Background->White]
Export[NotebookDirectory[]<>"Plots/Moran_Det_AF.pdf",%,ImageResolution->300];

Plot of heterozygosity as a function of time.

DetDynHet[pars_,inits_,color_]:=Show[ListLinePlot[Table[2pHsol[τ,pars,inits](1-pHsol[τ,pars,inits]),{τ,0,1000}],PlotRange->{0,0.5},AspectRatio->0.75,PlotStyle->color,Frame->True,FrameTicksStyle->Directive[Black,Bold,14](*,FrameLabel->{"Time","Host Heterozygosity"}, LabelStyle->Directive[Black,Bold, 18]*)],ListPlot[{{0,2pH0(1-pH0)}/.inits},PlotStyle->Red,PlotMarkers->{●,10}],ImageSize->{Automatic,300}]

Show[(*DetDynHet[TestPars,TestInits0,Black],DetDynHet[TestPars,TestInits2,Blue],*)DetDynHet[TestPars,TestInits,Black(*Darker[Green]*)],Background->White]
Export[NotebookDirectory[]<>"Plots/Moran_Det_Het.pdf",%,ImageResolution->300];

Neutral Expectation

For the neutral Moran model, heterozygosity after k events (1/2)(((1-(2/(κ^2))))^k)  given the initial heterozygosity of (Subscript[H, 0])=(1/2).  Probability of k events by time (((((κ t))^k)Exp[-κ t])/(k!))

Clear[neutExp]
neutExp[κ_,t_]:=If[t==0,(1/2),Sum[(((((κ t))^k)Exp[-κ t])/(k!))(1/2)(((1-(2/(κ^2))))^k),{k,0,∞}]]

This simplifies to:

Sum[(((((κ t))^k)Exp[-κ t])/(k!))(1/2)(((1-(2/(κ^2))))^k),{k,0,∞}]

Ensemble Moments

Substitutions

Let δX be the deviation of variable X from its value at the deterministic equilibrium.  
Let μδX be the expectation of deviation of X from equilibrium over all realizations of the stochastic process.
Let CδXY, CδXYZ,...  be the correlation, skew and higher order centralized moments of the variables X,Y,Z,...

Moment Substitutions

Substitution lists that allows us to take the expectation of by substituting higher order moments in terms of lower order moments and the corresponding centralized moment.

first={δH->μδH,δP->μδP};

second=Block[{out},out={};
AppendTo[out,Solve[CδHH==Expand[(((δH-μδH))^2)]/.((δH)^2)->X,X][[1]]/.first/.X->((δH)^2)];
AppendTo[out,Solve[CδHP==Expand[(δH-μδH)(δP-μδP)]/.δH δP->X,X][[1]]/.first/.X->δH δP];
AppendTo[out,Solve[CδPP==Expand[(((δP-μδP))^2)]/.((δP)^2)->X,X][[1]]/.first/.X->((δP)^2)];
Flatten[out]
];

third=Block[{out},out={};
AppendTo[out,Solve[CδHHH==Expand[(((δH-μδH))^3)]/.((δH)^3)->X,X][[1]]/.second/.first/.X->((δH)^3)];
AppendTo[out,Solve[CδHHP==Expand[(((δH-μδH))^2)(δP-μδP)]/.((δH)^2) δP->X,X][[1]]/.second/.first/.X->((δH )^2)δP];
AppendTo[out,Solve[CδHPP==Expand[(δH-μδH)(((δP-μδP))^2)]/.δH ((δP)^2)->X,X][[1]]/.second/.first/.X->δH ((δP)^2)];
AppendTo[out,Solve[CδPP==Expand[(((δP-μδP))^3)]/.((δP)^3)->X,X][[1]]/.second/.first/.X->((δP)^3)];
Flatten[out]
];

fourth=Block[{out},out={};
AppendTo[out,Solve[CδHHHH==Expand[(((δH-μδH))^4)]/.((δH)^4)->X,X][[1]]/.third/.second/.first/.X->((δH)^4)];
AppendTo[out,Solve[CδHHHP==Expand[(((δH-μδH))^3)(δP-μδP)]/.((δH)^3) δP->X,X][[1]]/.third/.second/.first/.X->((δH )^3)δP];
AppendTo[out,Solve[CδHHPP==Expand[(((δH-μδH))^2)(((δP-μδP))^2)]/.((δH)^2) ((δP)^2)->X,X][[1]]/.third/.second/.first/.X->((δH )^2)((δP)^2)];
AppendTo[out,Solve[CδHPPP==Expand[(δH-μδH)(((δP-μδP))^3)]/.δH ((δP)^3)->X,X][[1]]/.third/.second/.first/.X->δH ((δP)^3)];
AppendTo[out,Solve[CδPPPP==Expand[(((δP-μδP))^4)]/.((δP)^4)->X,X][[1]]/.third/.second/.first/.X->((δP)^4)];
Flatten[out]
];

Simplification

Substitution of notation for the counts.

nSub={nH[1]->H,nH[2]->κ-H,nP[1]->P,nP[2]->κ-P};

Substituting counts into (and back from) deviations from the equilibrium.

Let δH=H-(OverHat[H]) and δP=P-(OverHat[P]).  We will assume that the deviations from equilibrium, δH and δP, are small.

Note that the covariance between the deviations from equilibrium are equal to the the covariance in the original variables.

eSub={H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP};

μeSubRev={μδH->μH[τ]-(κ/2),μδP->μP[τ]-(κ/2),CδHH-> CHH[τ],CδHP->CHP[τ],CδPP->CPP[τ],CδHHH->CHHH[τ],CδHHP->CHHP[τ],CδHPP-> CHPP[τ],CδPPP->CPPP[τ]};

Moment Closure

assump3={CδHHHH->0,CδHHHP->0,CδHHPP->0,CδHPPP->0,CδPPPP->0};

Rates f[x] and change in x

d[H_,P_]:=((Sum[α β[h,p]nH[h]nP[p],{h,1,2},{p,1,2}]+Sum[δ nH[h],{h,1,2}])/κ)/.nSub/.βSub//FullSimplify

f[H_,P_]:=Join[(*Host Free-living death*)Flatten[Table[((δ nH[h1])/(d[H,P]))((nH[h2])/κ)/.nSub,{h1,1,2},{h2,1,2}]],(*Infectious Death of h1 by p1, h2 born,p2 born*)Flatten[Table[((α β[h1,p1])/(d[H,P]))nH[h1]nP[p1]((nH[h2])/κ)((nP[p2])/κ)/.nSub/.βSub,{h1,1,2},{p1,1,2},{h2,1,2},{p2,1,2}]],(*Pathogen Free-living death*)Flatten[Table[((γ nP[p1])/(d[H,P]))((nP[p2])/κ)/.nSub,{p1,1,2},{p2,1,2}]]]//FullSimplify

ΔH=Join[(*Host Free-living death*)Flatten[Table[-Mod[h1,2]+Mod[h2,2],{h1,1,2},{h2,1,2}],1],(*Infectious Death of h1 by p1, h2 born,p2 dies*)Flatten[Table[-Mod[h1,2]+Mod[h2,2],{h1,1,2},{p1,1,2},{h2,1,2},{p2,1,2}]],(*Pathogen Free-living death*)Flatten[Table[0,{p1,1,2},{p2,1,2}]]]

ΔP=Join[(*Natural death of h1 birth of h2*)Flatten[Table[0,{h1,1,2},{h2,1,2}],1],(*Infectious Death of h1 by p1, h2 born,p2 dies*)Flatten[Table[Mod[p1,2]-Mod[p2,2],{h1,1,2},{p1,1,2},{h2,1,2},{p2,1,2}]],(*Natural death of p1 birth of p2*)Flatten[Table[-Mod[p1,2]+Mod[p2,2],{p1,1,2},{p2,1,2}]]]

Differential Equations

Means

Clear[dμHdt]
dμHdt[t_]:=dμHdt[t]=FullSimplify[Collect[Normal[Series[Simplify[f[H,P].ΔH]/.eSub,{ϵ,0,2}]],{ϵ,δH,δP},FullSimplify]/.third/.second/.first/.{ϵ->1}]/.μeSubRev/.{τ->t}//FullSimplify

Clear[dμPdt]
dμPdt[t_]:=dμPdt[t]=FullSimplify[Collect[Normal[Series[Simplify[f[H,P].ΔP]/.eSub,{ϵ,0,2}]],{ϵ,δH,δP},FullSimplify]/.second/.first/.{ϵ->1}]/.μeSubRev/.{τ->t}//FullSimplify

Variances

Clear[dHHdt]
dHHdt[t_]:=dHHdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(((H+ΔH))^2)-(H^2)]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,3}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}(*FullSimplify[Collect[Normal[Series[Simplify[f[H,P].Simplify[(((H+ΔH))^2)-(H^2)]/.eSub],{ϵ,0,3}]],{ϵ,δH,δP},FullSimplify]/.third/.second/.first/.{ϵ->1}]/.μeSubRev/.{τ->t}//FullSimplify*)

Clear[dHPdt]
dHPdt[t_]:=dHPdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(H+ΔH)(P+ΔP)-H P]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,3}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}

Clear[dPPdt]
dPPdt[t_]:=dPPdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(((P+ΔP))^2)-(P^2)]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,3}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}

Clear[dCHHdt,dCHPdt,dCPPdt]
dCHHdt[t_]:=(*dCHHdt[i,j,M,t]=*)dHHdt[t]-μH[t]dμHdt[t]-μH[t]dμHdt[t]
dCHPdt[t_]:=(*dCHPdt[i,j,M,t]=*)dHPdt[t]-μH[t]dμPdt[t]-μP[t]dμHdt[t]
dCPPdt[t_]:=(*dCPPdt[i,j,M,t]=*)dPPdt[t]-μP[t]dμPdt[t]-μP[t]dμPdt[t]

Series[dμPdt[t]/.κ->(1/ϵ),{ϵ,0,-1}]//Simplify

Series[dCHHdt[t]/.κ->(1/ϵ),{ϵ,0,-1}]//Simplify

Series[dPPdt[t]/.κ->(1/ϵ),{ϵ,0,-1}]//Simplify

Skews

Clear[dHHHdt]
dHHHdt[t_]:=dHHHdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(((H+ΔH))^3)-(H^3)]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,4}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}

Clear[dHHPdt]
dHHPdt[t_]:=dHHPdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(((H+ΔH))^2)(P+ΔP)-(H^2) P]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,4}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}

Clear[dHPPdt]
dHPPdt[t_]:=dHPPdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(H+ΔH)(((P+ΔP))^2)-H (P^2)]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,4}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}

Clear[dPPPdt]
dPPPdt[t_]:=dPPPdt[t]=Expand[Normal[Series[Simplify[f[H,P].Simplify[(((P+ΔP))^3)-(P^3)]/.{H->(κ/2)+ϵ δH,P->(κ/2)+ϵ δP}],{ϵ,0,4}]]]/.fourth/.third/.second/.first/.{ϵ->1}/.assump3/.μeSubRev/.{τ->t}

Clear[dCHHHdt]
dCHHHdt[t_]:=dCHHHdt[t]=dHHHdt[t]-μH[t]dCHHdt[t]-μH[t]dCHHdt[t]-μH[t]dCHHdt[t]-dμHdt[t]CHH[t]-dμHdt[t]CHH[t]-dμHdt[t]CHH[t]-μH[t]μH[t]dμHdt[t]-μH[t]μH[t]dμHdt[t]-μH[t]μH[t]dμHdt[t]
Clear[dCHHPdt]
dCHHPdt[t_]:=dCHHPdt[t]=dHHPdt[t]-μH[t]dCHPdt[t]-μH[t]dCHPdt[t]-μP[t]dCHHdt[t]-dμHdt[t]CHP[t]-dμHdt[t]CHP[t]-dμPdt[t]CHH[t]-μH[t]μH[t]dμPdt[t]-μH[t]μP[t]dμHdt[t]-μH[t]μP[t]dμHdt[t]
Clear[dCHPPdt]
dCHPPdt[t_]:=dCHPPdt[t]=dHPPdt[t]-μH[t]dCPPdt[t]-μP[t]dCHPdt[t]-μP[t]dCHPdt[t]-dμHdt[t]CPP[t]-dμPdt[t]CHP[t]-dμPdt[t]CHP[t]-μH[t]μP[t]dμPdt[t]-μH[t]μP[t]dμPdt[t]-μP[t]μP[t]dμHdt[t]
Clear[dCPPPdt]
dCPPPdt[t_]:=dCPPPdt[t]=dPPPdt[t]-μP[t]dCPPdt[t]-μP[t]dCPPdt[t]-μP[t]dCPPdt[t]-dμPdt[t]CPP[t]-dμPdt[t]CPP[t]-dμPdt[t]CPP[t]-μP[t]μP[t]dμPdt[t]-μP[t]μP[t]dμPdt[t]-μP[t]μP[t]dμPdt[t]

Genetic Variance

D[(2/κ)μH[t]-(2/(κ^2))(((μH[t])^2)+CHH[t]),t]

dHetdt[t_]:=((2 dμHdt[t])/κ)-((2 (dCHHdt[t]+2 μH[t] dμHdt[t]))/(κ^2))

nearEquSub={μH[t]->(κ/2)+ϵ μdH[t],μP[t]->(κ/2)+ϵ μdH[t],CHH[t]->0+ϵ CδHH[t],CHP[t]->0+ϵ CδHP[t],CPP[t]->0+ϵ CδPP[t],CHHH[t]->0+ϵ CδHHH[t],CHHP[t]->0+ϵ CδHHP[t],CHPP[t]->0+ϵ CδHPP[t],CPPP[t]->0+ϵ CδPPP[t]};

If we are very near equilibrium then the ensemble moment dynamics of heterozygosity are equal to that of the neutral expectation

Normal[Series[dHetdt[t]/.nearEquSub/.{κ->(1/ϵ)},{ϵ,0,1}]]/.{ϵ->(1/κ)}

DSolve[{D[HEM[t],t]==-(1/κ),HEM[0]==(1/2)},{HEM[t]},t]//Simplify

Normal[Series[D[(1/2) (E^(-((2 t)/κ))),t]/.{t->ϵ τ,κ->(1/ϵ)},{ϵ,0,1}]]/.{ϵ->(1/κ)}

Numerical Integration

Functions

Clear[odesEM]

TestPars={X->0.9,Y->0.7,α->0.2,γ->0.5,δ->0.8,κ->1000};

initsSubEM[pars_]:={μH[0]->(κ/2),μP[0]->(κ/2),CHH[0]->0,CHP[0]->0,CPP[0]->0,CHHH[0]->0,CHHP[0]->0,CHPP[0]->0,CPPP[0]->0}/.pars;

initsEM[pars_]:={μH[0]==(κ/2),μP[0]==(κ/2),CHH[0]==0,CHP[0]==0,CPP[0]==0,CHHH[0]==0,CHHP[0]==0,CHPP[0]==0,CPPP[0]==0}/.pars;

odesEM[pars_]:={D[μH[t],t]==dμHdt[t],D[μP[t],t]==dμPdt[t],D[CHH[t],t]==dCHHdt[t],D[CHP[t],t]==dCHPdt[t],D[CPP[t],t]==dCPPdt[t],D[CHHH[t],t]==dCHHHdt[t],D[CHHP[t],t]==dCHHPdt[t],D[CHPP[t],t]==dCHPPdt[t],D[CPPP[t],t]==dCPPPdt[t]}/.pars;

varsEM={μH[t],μP[t],CHH[t],CHP[t],CPP[t],CHHH[t],CHHP[t],CHPP[t],CPPP[t]};

Clear[nsolEM]
nsolEM[pars_,tMax_]:=nsolEM[pars,tMax]=NDSolve[Join[odesEM[pars],initsEM[pars]],varsEM,{t,0,tMax}][[1]]
nHet[pars_,tMax_,τ_]:=(2/κ)μH[t]-(2/(κ^2))(((μH[t])^2)+CHH[t])/.nsolEM[pars,tMax]/.pars/.{t->τ}

Plotting

Effect of deviation

Parsγ[γin_,δin_]:={X->0.9,Y->0.7,α->0.2,γ->γin,δ->δin,κ->100000};

Dequ=((Sum[β[h,p]α ((((κ/2)))^2),{h,1,2},{p,1,2}]+δ κ)/κ)/.βSub;

ΔHEM[pars_,τ_]:=nHet[pars,τ,τ]-neutExp[κ/.pars,τ]

γDeivEM=Block[{Xr,Yr,αr,γr,δr,out,ct,pars},out={};
For[ct=1,ct≤200,ct++,
Xr=RandomReal[];Yr=RandomReal[{0,Xr}]; αr=RandomReal[{0,0.5}];
γr=RandomReal[{0.1,1}];δr=RandomReal[{0.1,1}];
pars={X->Xr,Y->Yr,α->αr,γ->γr,δ->δr,κ->100000};
(*AppendTo[out,{γr,δr,Dequ/.pars,ΔHEM[pars,50]}]*)
AppendTo[out,{γr,δr,Dequ/.Parsγ[γr,δr],ΔHEM[Parsγ[γr,δr],50]}]
];
out
];

Figure 3 in the main text. Note that due to random sampling of parameters, plot may not be 100% identical.

γDeivEMPlot=Module[{data},data=Transpose[{γDeivEM[[;;,1]]-γDeivEM[[;;,2]],γDeivEM[[;;,4]]}];
(*lmEM=LinearModelFit[data,x,x];*)
lmEM=NonlinearModelFit[data,y0+(L/(1+Exp[-k(x-x0)])),{y0,L,k,x0},x];
Show[ListPlot[data,Frame->True,PlotStyle->Black],Plot[lmEM[x],{x,Min[data[[;;,1]]],Max[data[[;;,1]]]},PlotStyle->Gray],PlotRange->{{-1,1},Automatic},Epilog->Text[Style["(R^2)= "<>ToString[Round[lmEM["RSquared"],0.001]],14,Black],Scaled[{0.85,0.75}]],ImageSize->{Automatic,300},FrameTicksStyle->Directive[Black,Bold,14]]
]
Export[NotebookDirectory[]<>"Plots/gammaDeivEM.pdf",γDeivEMPlot,ImageResolution->500];

Individual Based Simulations (C++)

Making Input files

Choosing Parameters

randPars[npts_(*,μIn_*)]:=Block[{out,ct,a,x,y,d,(*m,*)g},out={};
For[ct=1,ct≤npts,ct++,
d=RandomReal[{0.1,1}];a=RandomReal[{0,0.5}];x=RandomReal[];y=RandomReal[{0,x}];g=RandomReal[{0.1,1}];(*m=RandomInteger[{1,5}];*)
AppendTo[out,{X->x,Y->y,α->a,δ->d,γ->g(*,M->m,μ->μIn*)}];
];
out
]

pars1=randPars[100(*,0.001*)];

Making Input File

file[pars_]:="input file for Moran MAM simulations

number of reps (nRep): 7
number of demes (nDemes): 1000
run time (tmax): 500
sampling frequency (sampFreq): 10
matching infection rate (X): "<>ToString[X/.pars]<>"
mismatching infection rate (Y): "<>ToString[Y/.pars]<>"
death probability (alpha): "<>ToString[α/.pars]<>"
host natural death rate (delta): "<>ToString[δ/.pars]<>"
pathogen natural death rate (gamma): "<>ToString[γ/.pars](*<>"
pathogen births per host (M): "<>ToString[M/.pars]<>"
pathogen mutation rate (mu): "<>ToString[μ/.pars]*);

exportFile[run_]:=Export["/Users/ailenemacpherson/Documents/VisualStudio/ConstantMAM/inputConst_"<>ToString[run]<>".txt",file[pars1[[run]]]]

Do[exportFile[run],{run,1,100}]

pars1[[1]]

Importing Data

Importing Data

File paths to simulation files.

file=NotebookDirectory[]<>"Sims/Moran/runConst_";
file2=NotebookDirectory[]<>"Sims/Moran_10_15/runConst_";

Clear[pars]
pars[file_,run_]:=pars[file,run]=Import[file<>ToString[run]<>"/out.csv"]

Clear[nDemes,nRep,tmax,δ,X,Y,α,κ,γ]
nRep[file_,run_]:=pars[file,run][[2,2]]
nDemes[file_,run_]:=pars[file,run][[3,2]]
tmax[file_,run_]:=pars[file,run][[4,2]]
sampFreq[file_,run_]:=pars[file,run][[5,2]]
X[file_,run_]:=pars[file,run][[6,2]]
Y[file_,run_]:=pars[file,run][[7,2]]
α[file_,run_]:=pars[file,run][[8,2]]
δ[file_,run_]:=pars[file,run][[9,2]]
γ[file_,run_]:=pars[file,run][[10,2]]
(*M[file_,run_]:=pars[file,run][[11,2]]
μ[file_,run_]:=pars[file,run][[12,2]]*)
κ[file_,run_,rep_]:=pars[file,run][[12+rep,2]];

simPars[file,run,rep] gives a substitution list of the simulation parameters for simulations in the folder “file”, for parameter condition specified by “run” (run ranges between 1 and 100)., and growth limit κ specified by rep (rep ranges between 0 and 6).

simPars[file_,run_,rep_]:={X->X[file,run],Y->Y[file,run],α->α[file,run],δ->δ[file,run],γ->γ[file,run],(*M->M[file,run],*)κ->κ[file,run,rep]};

While no-longer strictly necessary fileExist is a useful troubleshooting tool that allows us to exclude cases where the simulation file doesn’t exist.

fileExist[file_,run_,rep_,t_]:=FileExistsQ[file<>ToString[run]<>"/rep_"<>ToString[rep]<>"/state_"<>ToString[t]<>".csv"]

repData is a {1000,8} matrix.  each row specifies a deme.  Columns 1 through 4 given the number of hosts of type 1, number of hosts of type 2, number of pathogens of type 1, and number of pathogens of 2 under host-parasite coevolution. Columns 5-8 are the analogous values under neutral genetic drift.

Clear[repData]
repData[file_,run_,rep_,t_]:=repData[file,run,rep,t]=Block[{name},
name=file<>ToString[run]<>"/rep_"<>ToString[rep]<>"/state_"<>ToString[t]<>".csv";
If[FileExistsQ[name]==True,Import[name]]
]

Processing

For a given parameter condition (run and rep) the following three functions determine in which demes coevolution is ongoing (CoevSel), in which directional selection occurs (DirSel), and in which the host is fixed (fixed).

Clear[CoevSel]
CoevSel[file_,run_,rep_,t_]:=CoevSel[file,run,rep,t]=If[fileExist[file,run,rep,t],Quiet[Flatten[Position[repData[file,run,rep,t],_?((#[[1]]>0&&#[[2]]>0)&&(#[[3]]>0&&#[[4]]>0)&)]]]]

Clear[DirSel]
DirSel[file_,run_,rep_,t_]:=DirSel[file,run,rep,t]=If[fileExist[file,run,rep,t],Quiet[Flatten[Position[repData[file,run,rep,t],_?((#[[1]]>0&&#[[2]]>0)&&Xor[#[[3]]>0,#[[4]]>0]&)]]]]

Clear[fixed]
fixed[file_,run_,rep_,t_]:=fixed[file,run,rep,t]=If[fileExist[file,run,rep,t],Quiet[Flatten[Position[repData[file,run,rep,t],_?(Xor[#[[1]]>0.,#[[2]]>0.]&)]]]]

Functions

Equilibrium

Clear[Dequ]
Dequ[file_,run_,rep_]:=((Sum[β[h,p]α ((((κ/2)))^2),{h,1,2},{p,1,2}]+δ κ)/κ)/.βSub/.simPars[file,run](*/.{κ-> 2+rep 0.25}*)

Bootstrap Functions

Takes nBoot bootstrap samples of data and calculates the mean of each.

BootstrapMean[data_,nBoot_]:=BootstrapMean[data,nBoot]=Block[{out,b},out={};
For[b=1,b≤nBoot,b++,
AppendTo[out,Mean[RandomChoice[data,Length[data]]]]
];
out]

Returns a bootstrap sample of data indexed by a dummy variable.

Bootstrap[data_,dummy_]:=Bootstrap[data,dummy]=RandomChoice[data,Length[data]]

Allele Frequency

Host allele frequency under coevolution

Clear[pH]
pH[file_,run_,rep_,t_]:=pH[file,run,rep,t]=If[fileExist[file,run,rep,t],((repData[file,run,rep,t][[;;,1]])/(Total[Transpose[repData[file,run,rep,t][[;;,{1,2}]]]]))//N]

Host allele frequency under neutral drift

Clear[pHn]
pHn[file_,run_,rep_,t_]:=pHn[file,run,rep,t]=If[fileExist[file,run,rep,t],((repData[file,run,rep,t][[;;,5]])/(Total[Transpose[repData[file,run,rep,t][[;;,{5,6}]]]]))//N]

Pathogen allele frequency under coevolution

Clear[pP]
pP[file_,run_,rep_,t_]:=pP[file,run,rep,t]=If[fileExist[file,run,rep,t],((repData[file,run,rep,t][[;;,3]])/(Total[Transpose[repData[file,run,rep,t][[;;,{3,4}]]]]))//N ]

Heterozygosity

Genetic variation (measured as heterozygosity) at time t simulated under coevolutionary selection.

Clear[GVar]
GVar[file_,run_,rep_,t_]:=GVar[file,run,rep,t]=If[fileExist[file,run,rep,t],2pH[file,run,rep,t](1-pH[file,run,rep,t])]

Genetic variation (measured as heterozygosity) at time t simulated under neutral genetic drift

Clear[GVarn]
GVarn[file_,run_,rep_,t_]:=GVarn[file,run,rep,t]=If[fileExist[file,run,rep,t],2pHn[file,run,rep,t](1-pHn[file,run,rep,t])]

Analytical expectation of genetic variation (measured as heterozygosity) at time t.

Clear[GVarn2]
GVarn2[file_,run_,rep_,t_]:=GVarn2[file,run,rep,t]=neutExp[κ/.simPars[file,run,rep],t]//N

Clear[ΔHet]
ΔHet[file_,run_,rep_,t_]:=If[fileExist[file,run,rep,t],Mean[GVar[file,run,rep,t]]-GVarn2[file,run,rep,t](*Mean[GVarn[file,run,rep,t]]*)]

Plotting

Genetic Variation

Dynamics

Plot of the dynamics for a given parameter condition (specified by file,run, and rep) up to time tMax.  Note that you have to run the Ensemble moment dynamics to calculate nHet, the dynamics of heterozygosity expected under that approximation.

dynPlot[file_,run_,rep_,tMax_]:=Block[{Δt,tList},Δt=sampFreq[file,run];tList=Join[Table[t,{t,0(*10*),tMax,Δt}]];
Show[
(*Coevolutionary heterozygosity*)
ListLinePlot[{Table[{t,Mean[GVar[file,run,rep,t]]},{t,tList}],Table[{t,Mean[GVar[file,run,rep,t]]+2((StandardDeviation[GVar[file,run,rep,t]])/(Sqrt[Length[GVar[file,run,rep,t]]]))},{t,tList}],Table[{t,Mean[GVar[file,run,rep,t]]-2((StandardDeviation[GVar[file,run,rep,t]])/(Sqrt[Length[GVar[file,run,rep,t]]]))},{t,tList}]},Filling->2->{3},PlotStyle->{Red,Pink,Pink},PlotRange->{0,0.5}],
(*Neutral heterozygosity from simulations*)ListLinePlot[{Table[{t,Mean[GVarn[file,run,rep,t]]},{t,tList}],Table[{t,Mean[GVarn[file,run,rep,t]]+2((StandardDeviation[GVarn[file,run,rep,t]])/(Sqrt[Length[GVarn[file,run,rep,t]]]))},{t,tList}],Table[{t,Mean[GVarn[file,run,rep,t]]-2((StandardDeviation[GVarn[file,run,rep,t]])/(Sqrt[Length[GVarn[file,run,rep,t]]]))},{t,tList}]},Filling->2->{3},PlotStyle->{Directive[Thick,Black],Gray,Gray},PlotRange->{0,0.5}],
(*Analytical neutral expectation*)
ListLinePlot[Table[{t,neutExp[κ[file,run,rep],t]},{t,tList}],PlotRange->{0,0.5},PlotStyle->Directive[Thick,Blue]],
(*Ensemble Moment Dynamics*)
Plot[nHet[simPars[file,run,rep],60,τ],{τ,0,60},PlotStyle->Directive[{Thickness[0.008],Darker[Green]}]],PlotRange->{{0,tMax},{Floor[neutExp[κ[file,run,rep],tMax],0.1],0.5}},Frame->True,FrameTicks->{{{0,0.1,0.2,0.3,0.4,0.5},False},{True,False}},ImageSize->{400,Automatic}(*,Epilog->Inset[barplot[run,rep,tMax,2],Scaled[{0.49,0.12}]]*)]]

Figure 4B (and inset) in the main text.

Show[dynPlot[file2,1,3,60],ImageSize->{Automatic,300},Frame->True,FrameTicksStyle->Directive[Black,Bold,14],FrameTicks->{{True,False},{True,False}},Background->White]
Export[NotebookDirectory[]<>"Plots/EgDynSub.pdf",%];

Show[dynPlot[file2,1,3,500],ImageSize->{Automatic,300},Frame->True,FrameTicksStyle->Directive[Black,Bold,14],FrameTicks->{{True,False},{True,False}},Background->White]
Export[NotebookDirectory[]<>"Plots/EgDynFull.pdf",%];

simPars[file2,1,3]

(10^2.75)

ΔH Distribution

An alternative way of exploring the effect of coevolution on host heterozygosity by focusing on the distributions in ΔH, the change in heterozygosity between the coevolutionary and neutral simulations.

ΔHList[file_,rep_,t_]:=Select[Table[ΔHet[file,run,rep,t],{run,1,60}],#≠"Null"&]

HetBW[file_,nReps_,t_]:=Block[{sigData,rep,temp,point},sigData={};
point=Graphics[{EdgeForm[Black],White,Disk[Scaled[{0.5,0.5}],0.02]}];
(*Which reps are signifcantly DIFFERENT From 0*)
For[rep=0,rep≤nReps,rep++,
If[Mean[BootstrapMean[ΔHList[file,rep,t],1000]]+1.6 StandardDeviation[BootstrapMean[ΔHList[file,rep,t],1000]]<0 (*||Mean[BootstrapMean[ΔHList[file,rep,t],10000]]-2 StandardDeviation[BootstrapMean[ΔHList[file,rep,t],10000]]>0*),
AppendTo[sigData,{{rep+1,0.13}}],If[Mean[BootstrapMean[ΔHList[file,rep,t],1000]]-1.6 StandardDeviation[BootstrapMean[ΔHList[file,rep,t],1000]]>0,AppendTo[sigData,{{rep+1,-0.45}}],AppendTo[sigData,{{rep+1,-0.6}}]]]
];
Show[
(*Adds stars for all distributions where the mean differs significantly from 0*)
ListPlot[sigData,PlotStyle->Table[ColorData["CMYKColors"][(rep/nReps)],{rep,0,nReps}],PlotMarkers->✶],
(*Violins*)
DistributionChart[Table[ΔHList[file,rep,t],{rep,0,nReps}],ChartStyle->Table[ColorData["CMYKColors"][(rep/nReps)],{rep,0,nReps}]],
(*Neutral line*)
ListLinePlot[{{0.2,0},{8,0}},PlotStyle->Directive[Black]],
(*Data points*)ListPlot[Table[Transpose[{RandomReal[{1+rep-0.1,1+rep+0.1},Length[ΔHList[file,rep,t]]],ΔHList[file,rep,t]}],{rep,0,6}],PlotMarkers->point],
(*Plot Options*)
PlotRange->{-0.35,Automatic},AxesOrigin->{Automatic,-0.5},Frame->True,FrameLabel->{"(Subscript[Log, 10])[κ]","ΔH"},FrameTicks->{{True,False},{Table[{rep+1,2+rep*0.25},{rep,0,nReps}],False}},LabelStyle->Directive[Black,Bold, 18],FrameTicksStyle->Directive[14],ImageSize->Large]
]

HetBW[file,6,500]
Export[NotebookDirectory[]<>"Plots/Moran_HetBW.gif",%,ImageResolution->300];

Mean Heterozygosity vs. Neutral Expectation

LOESS Smoothing

Project data onto the X=Y line.

Clear[dataNew]
dataNew[data_]:=dataNew[data]=Transpose[{Cos[Abs[ArcTan[((data[[;;,2]])/(data[[;;,1]]))]-(π/4)]]*Sqrt[Total[Transpose[(data^2)]]],Sin[ArcTan[((data[[;;,2]])/(data[[;;,1]]))]-(π/4)]*Sqrt[Total[Transpose[(data^2)]]]}]

A gaussian weighting of the data  near the point xNew along the X=Y line

Clear[weights]
weights[xNew_,αW_,data_]:=weights[xNew,αW,data]=Exp[-((((dataNew[data][[;;,1]]-xNew))^2)/(2 ((αW)^2)))];

Find the weighted mean and standard deviation and standard error.

Clear[meanNew]
meanNew[xNew_,αW_,data_]:=meanNew[xNew,αW,data]=Total[weights[xNew,αW,data]dataNew[data][[;;,2]]]/Total[weights[xNew,αW,data]]

SDNew[xNew_,αW_,data_]:=Sqrt[Total[weights[xNew,αW,data]*(((dataNew[data][[;;,2]]-meanNew[xNew,αW,data]))^2)]/Total[weights[xNew,αW,data]]]

SENew[xNew_,αW_,data_]:=((SDNew[xNew,αW,data])/(Sqrt[Total[weights[xNew,αW,data]]]))

The rolling average fit along the X+Y line (rotated plot)

Clear[RollFit]
RollFit[αW_,npts_,data_]:=RollFit[αW,npts,data]=Table[{x,meanNew[x,αW,data]},{x,0.001(*((Sqrt[2*(0.5^2)])/(npts+2))*),Sqrt[2*(0.5^2)](*-((Sqrt[2*(0.5^2)])/(npts+2))*),((Sqrt[2*(0.5^2)]-0.001)/(npts+2))}];

Projecting back into the original coordinates.

RollFitOld[αW_,npts_,data_]:=Transpose[{Cos[ArcTan[((RollFit[αW,npts,data][[;;,2]])/(RollFit[αW,npts,data][[;;,1]]))]+(π/4)]Sqrt[((RollFit[αW,npts,data][[;;,2]])^2)+((RollFit[αW,npts,data][[;;,1]])^2)],Sin[ArcTan[((RollFit[αW,npts,data][[;;,2]])/(RollFit[αW,npts,data][[;;,1]]))]+(π/4)]Sqrt[((RollFit[αW,npts,data][[;;,2]])^2)+((RollFit[αW,npts,data][[;;,1]])^2)]}]

Plotting

GVarList[file_,t_]:=Block[{run,rep,out},out={};
For[rep=0,rep≤6,rep++,
AppendTo[out,{}];
For[run=1,run≤60,run++,
If[fileExist[file,run,rep,t],AppendTo[out[[-1]],{GVarn2[file,run,rep,t](*Mean[GVarn[file,run,rep,t]]*),Mean[GVar[file,run,rep,t]]}]]
]
];
out
]

(*Export[NotebookDirectory[]<>"MoranMean.csv",RollFitOld[0.1,40,Flatten[test[file,250],1]]];*)

Identifying which data points differ significantly from the neutral expectation (given the Bonferroni correction).

ΔHrep[file_,run_,rep_,t_]:=GVar[file,run,rep,t]-GVarn[file,run,rep,t]
Clear[sigΔH]
sigΔH[file_,t_]:=sigΔH[file,t]=Block[{outSig,outNotSig,run,rep,CIU,CIL},outSig={};outNotSig={};
For[rep=0,rep≤6,rep++,
AppendTo[outSig,{}];AppendTo[outNotSig,{}];
For[run=1,run≤100,run++,
If[fileExist[file,run,rep,t],
CIU=Mean[ΔHrep[file,run,rep,t]]+(*1.96*)2.64 StandardDeviation[BootstrapMean[ΔHrep[file,run,rep,t],1000]];
CIL=Mean[ΔHrep[file,run,rep,t]]-(*1.96*)2.64 StandardDeviation[BootstrapMean[ΔHrep[file,run,rep,t],1000]];
If[CIU>0&&CIL<0,AppendTo[outNotSig[[-1]],{(*GVarn2[file,run,rep,t]*)Mean[GVarn[file,run,rep,t]],Mean[GVar[file,run,rep,t]]}],AppendTo[outSig[[-1]],{(*GVarn2[file,run,rep,t]*)Mean[GVarn[file,run,rep,t]],Mean[GVar[file,run,rep,t]]}]]]
]
];
{outNotSig,outSig}
]

Percent of data points that did not differ significantly from the neutral expectation.

PercentNeut[file_,rep_,t_]:=Length[sigΔH[file,t][[1,rep+1]]]

Figure 4A in the main text.

MeanFit[file_,t_]:=Show[ListLinePlot[Table[RollFitOld[0.1,40,Bootstrap[Flatten[GVarList[file,t],1],x]],{x,1,1000}],PlotStyle->Directive[Pink,Opacity[0.02]]],Plot[x,{x,0,0.5},PlotStyle->Directive[Black,Dashed]],
ListLinePlot[RollFitOld[0.1,40,Flatten[GVarList[file,t],1]],PlotStyle->Red],ListPlot[Select[sigΔH[file,t][[1]],#≠{}&],PlotStyle->Table[Lighter[ColorData["CMYKColors"][(rep/6)]],{rep,0,6}],PlotMarkers->○],ListPlot[Select[sigΔH[file,t][[2]],#≠{}&],PlotStyle->Table[Directive[Darker[ColorData["CMYKColors"][(rep/6)]]],{rep,0,6}],PlotMarkers->●],AspectRatio->1,PlotRange->{{0,0.5},{0,0.5}},ImageSize->{Automatic,500},Frame->True,FrameTicksStyle->Directive[Black,Bold,16],FrameTicks->{{True,False},{True,False}},Epilog->Table[Text[Style[ToString[PercentNeut[file,rep,t]],14,ColorData["CMYKColors"][(rep/6)]],{neutExp[(10^(2+0.25 rep)),t],neutExp[(10^(2+0.25 rep)),t]+0.07}(*Scaled[{0.1,0.1+((0.6-0.1)/7)(6-rep)}]*)],{rep,0,6}] ]

MeanFit[file2,250]
Export[NotebookDirectory[]<>"Plots/NeutExp.pdf",%,ImageResolution->300];

κLeg=PointLegend[Table[ColorData["CMYKColors"][(rep/6)],{rep,0,6}],Table[ToString[Round[(10^(2+0.25rep))]],{rep,0,6}],LegendFunction->"Frame",LegendLabel->Text[Style["κ",14,Bold]],LegendMarkerSize->15]
Export[NotebookDirectory[]<>"Plots/kappaLeg.pdf",κLeg];

lineLeg1=LineLegend[{{Black,Dashed},Red,{Pink,Opacity[0.5]}},{"Neutral Expectation","IBS Mean","Bootstrap Means"},LegendFunction->"Frame"]
Export[NotebookDirectory[]<>"Plots/lineLeg1.pdf",lineLeg1];

Effect of γ on heterozygosity

Although not presented in the main text we can observe the effect of perturbations on the deviation from the neutral expectation as expected from the ensemble moment approximation.

Clear[ΔHγList]

ΔHγList[file_,rep_,t_]:=Select[Table[{((γ[file,run])/(Dequ[file,run,rep])),ΔH[file,run,rep,t]},{run,1,60}] ,#[[2]]≠"Null"&]

γPlot[file_,rep_,t_]:=Block[{data,lm},
data=Select[ΔHγList[file,rep,t],(*#[[2]]<0&&*)#[[1]]<9&];
lm=LinearModelFit[data,x,x];
(*lm=NonlinearModelFit[data,a Exp[b x],{a,b},x];*)
R2[rep]=lm["RSquared"];
Show[ListPlot[data,PlotStyle->ColorData["CMYKColors"][(rep/6)],Frame->True,(*FrameLabel->{"(γ/D): Rate of Pathogen Deaths","ΔH"},*)ImageSize->Large],Plot[lm[x],{x,Min[data[[;;,1]]],Max[data[[;;,1]]]},PlotStyle->ColorData["CMYKColors"][(rep/6)]],PlotRange->All]
]

Show[Table[γPlot[file,rep,50],{rep,0,6}],Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[R2[rep],0.001]],14,ColorData["CMYKColors"][(rep/6)]],Scaled[{0.1,0.1+((0.6-0.1)/7)(6-rep)}]],{rep,0,6}],FrameTicksStyle->Directive[Black,Bold,14],ImageSize->{Automatic,300}(*,LabelStyle->Directive[Black,Bold, 18],FrameLabel->{"γ/(OverHat[D])",ΔH}*)]
Export[NotebookDirectory[]<>"Plots/gammaDeivIBS.pdf",%,ImageResolution->300];

Allele Frequency Dynamics

Clear[traj]
traj[file_,run_,rep_,tLim_]:=traj[file,run,rep,tLim]=Block[{outC,outD,outN,t,d},outC={};outD={};outN={};
For[d=1,d<500(*nDemes[run]*),d++,
AppendTo[outC,{}];AppendTo[outD,{}];AppendTo[outN,{}];
For[t=0,t≤tLim,t=t+sampFreq[file,run],
(*Print[{t,pH[run,rep,d,t],pP[run,rep,d,t]} ];*)
If[0<pH[file,run,rep,t][[d]]<1,If[0<pP[file,run,rep,t][[d]]<1,
(*Ongoing Coevolution*)AppendTo[outC[[-1]],{t,pH[file,run,rep,t][[d]]}],
(*Monomorphic pathogen*)AppendTo[outD[[-1]],{t,pH[file,run,rep,t][[d]]}]],
(*Fixed Host*)AppendTo[outN[[-1]],{t,pH[file,run,rep,t][[d]]}]
]
]
];
{outC,Select[outD,#≠{}&],Select[outN,#≠{}&]}//N
]

AFDyn[file_,run_,rep_,tLim_]:=Show[{ListLinePlot[traj[file,run,rep,tLim][[1]],PlotRange->{0,1},PlotStyle->Directive[Gray,Opacity[0.1]]],ListLinePlot[traj[file,run,rep,tLim][[2]],PlotRange->{0,1},PlotStyle->Directive[Red,Opacity[0.2]]],ListLinePlot[traj[file,run,rep,tLim][[3]],PlotRange->{0,1},PlotStyle->Directive[Blue,Opacity[0.2]]]},ImageSize->{Automatic,300},Frame->True,(*FrameLabel->{"Time","(Subscript[p, H])"},LabelStyle->Directive[Bold, Black,20],*)FrameTicksStyle->Directive[Black,Bold,14],FrameTicks->{{{0,0.25,0.5,0.75,1},False},{Table[100*i,{i,0,5}],False}},PlotRange->{-0.05,1.05}]

Figure 4C in the main text.

AFDyn[file,2,5,500]
Export[NotebookDirectory[]<>"Plots/Moran_AF_Dyn.pdf",%,ImageResolution->300];

simPars[file,2,5]

LineLegend[{Gray,Red,Blue},{"Coevolution","Directional Selection","Fixed"}]
Export[NotebookDirectory[]<>"Plots/AFDynLeg.pdf",%,ImageResolution->200];