A Wright-Fisher Model of host-pathogen coevolution

The deterministic limit (difference equations)

WH[1,t_]:=1-αH(X pP[t]+Y (1-pP[t]))
WH[2,t_]:=1-αH(X (1-pP[t])+Y pP[t])
WP[1,t_]:=1+αP(X pH[t]+Y (1-pH[t]))
WP[2,t_]:=1+αP(X (1-pH[t])+Y pH[t])

gH[t_]:=pH[t](((WH[1,t])/(WH[1,t]pH[t]+WH[2,t](1-pH[t]))))//FullSimplify
gP[t_]:=pP[t](((WP[1,t])/(WP[1,t]pP[t]+WP[2,t](1-pP[t]))))//FullSimplify

gH[t]-pH[t]//FullSimplify

gP[t]-pP[t]//FullSimplify

Equilibrium and Stability

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

jMtrx={{D[gH[t],pH[t]],D[gH[t],pP[t]]},{D[gP[t],pH[t]],D[gP[t],pP[t]]}}/.equi[[3]]//Simplify;

λs=FullSimplify[Eigenvalues[jMtrx]]

This can be rewritten as:
λ=1±(((X-Y) (Sqrt[-1])(Sqrt[(Subscript[α, H])(Subscript[α, P])(2-(X+Y)(Subscript[α, H]))(2+(X+Y) (Subscript[α, P]))]))/((2-(X+Y)(Subscript[α, H]))(2+(X+Y) (Subscript[α, P]))))

The magnitude of these is always greater than one. The eigenvalues are imaginary and hence the equilibrium characterized by unstable cycles.

The rate at which the cycles deviate from equilibrium is given by the magnitude of the eigenvalue.

λLen=Sqrt[(1^2)+((((((X-Y) (Sqrt[-1])(Sqrt[(Subscript[α, H])(Subscript[α, P])(2-(X+Y)(Subscript[α, H]))(2+(X+Y) (Subscript[α, P]))]))/((2-(X+Y)(Subscript[α, H]))(2+(X+Y) (Subscript[α, P]))))))^2)];

λLen

Numerical Solution

Clear[nSol]
nSol[tMax_,inits_,pars_]:=nSol[tMax,inits,pars]=Block[{out,τ},out={{pH0,pP0}/.inits};
For[τ=1,τ≤tMax,τ++,
AppendTo[out,{gH[t],gP[t]}/.pars/.{pH[t]->out[[-1,1]],pP[t]->out[[-1,2]]}]
];
out
]

TestPars={X->0.8,Y->0.4,α->0.4};
TestInits={pH0->0.53,pP0->0.47};

DetPlotAF[tMax_,inits_,pars_]:=Show[ListLinePlot[nSol[tMax,inits,pars],AspectRatio->1,PlotRange->{{0,1},{0,1}},PlotStyle->Black,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}]

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

DetPlotAF[1000,TestInits,TestPars]
Export[NotebookDirectory[]<>"Plots/WF_Det_AF.pdf",%,ImageResolution->300];

Show[DetPlotHet[1000,TestInits,TestPars](*,ListLinePlot[{{500,0},{500,0.5}},PlotStyle->Blue]*)]
Export[NotebookDirectory[]<>"Plots/WF_Det_Het.pdf",%,ImageResolution->300];

λLen

ContourPlot[λLen/.{α->0.2},{X,0,1},{Y,0,X},ColorFunction->GrayLevel,Frame->True,FrameLabel->{"X","Y"},FrameTicksStyle->Directive[Black,Bold,14], LabelStyle->Directive[Black,Bold, Large],ImageSize->Large]

Neutral Expectation

Hneut[n_,t_]:=0.5 (((1-(1/n)))^t)

Simulations

Simulations

Drawing Parameters

Clear[randPars]
randPars[nRun_,μIn_,tMaxIn_,nDemesIn_]:=randPars[nRun,μIn,tMaxIn,nDemesIn]=Block[{x,y,c,a,out,rep},out={};
For[c=1,c≤nRun,c++,
x=RandomReal[];y=RandomReal[{0,x}];a=RandomReal[{0,0.5}];
AppendTo[out,{}];
For[rep=0,rep<7,rep++,
AppendTo[out[[-1]],{X->x,Y->y,α->a,tMax->tMaxIn,nDemes->nDemesIn,κ->Round[(10^(2+rep*0.5)),2],μ->μIn}];
];
];
out
]

parsList=randPars[60,0.001,500,1000];
temp=Flatten[Table[{run,rep,X,Y,α,tMax,nDemes,κ,μ}/.parsList[[run,rep]],{run,1,60},{rep,1,7}],1];
parsTable=PrependTo[temp,{"run","rep","X","Y","alpha","tMax","nDemes","kappa","mu"}];
Export[NotebookDirectory[]<>"ConstSims/pars.csv",parsTable];

Running Simulations

export[run_,rep_]:=Block[{hostTab,pathTab},hostTab=simWF[parsList[[run,rep]]][[;;,;;,2]];pathTab=simWF[parsList[[run,rep]]][[;;,;;,3]];
 Export[NotebookDirectory[]<>"ConstSims/Host_run_"<>ToString[run]<>"_rep_"<>ToString[rep]<>".csv",hostTab];Export[NotebookDirectory[]<>"ConstSims/Path_run_"<>ToString[run]<>"_rep_"<>ToString[rep]<>".csv",pathTab];]

Clear[simWF]
simWF[pars_]:=simWF[pars]=Block[{out,d,t,WH,WP,nH1,nH2,nP1,nP2},out={};
WH[nPX_,nPY_]:=1-α(X (nPX/κ)+Y (nPY/κ))/.pars; 
WP[nHX_,nHY_]:=1+α(X (nHX/κ)+Y (nHY/κ))/.pars; 
For[d=1,d≤nDemes/.pars,d++,
AppendTo[out,{{0,Round[(κ/2)/.pars],Round[(κ/2)/.pars]}}];
For[t=1,t≤tMax/.pars,t++,
nH1=out[[-1,-1,2]];nH2=κ-out[[-1,-1,2]]/.pars;
nP1=out[[-1,-1,3]];nP2=κ-out[[-1,-1,3]]/.pars;
AppendTo[out[[-1]],{t,RandomVariate[BinomialDistribution[κ/.pars,((WH[nP1,nP2]nH1)/(WH[nP1,nP2]nH1+WH[nP2,nP1]nH2))]],RandomVariate[BinomialDistribution[κ/.pars,((WP[nH1,nH2]nP1 (1-μ/.pars)+WP[nH2,nH1]nP2(μ/.pars))/(WP[nH1,nH2]nP1+WP[nH2,nH1]nP2))]]}];
];
];
out//N
]

Clear[simWFAll]
simWFAll=Do[simWF[parsList[[run,rep]]];export[run,rep];Print[{run,rep}];,{run,1,60},{rep,1,7}]

Explaining variation

TestPars={X->0.9,Y->0.89,tMax->200,nDemes->1000,κ->100};
TestPars2={X->0.1,Y->0.09,tMax->200,nDemes->1000,κ->100};
TestPars3={X->0.9,Y->0.8955757084820807`,tMax->200,nDemes->1000,κ->100};

simWF[TestPars]

simWF[TestPars2]

simWF[TestPars3]

{λLen/.TestPars,TestPars}

{λLen/.TestPars2,TestPars2}

{λLen/.TestPars3,TestPars3}

ListLinePlot[{Mean[Hdyn[TestPars]],Mean[Hdyn[TestPars2]],Mean[Hdyn[TestPars3]]}]

GraphicsRow[{NPlot[TestPars],NPlot[TestPars2],NPlot[TestPars3]},ImageSize->Full]

GraphicsRow[{ListLinePlot[nSol[1000,{pH0->0.51,pP0->0.49},TestPars],AspectRatio->1,PlotRange->{{0.45,0.55},{0.45,0.55}}],ListLinePlot[nSol[1000,{pH0->0.51,pP0->0.49},TestPars2],AspectRatio->1,PlotRange->{{0.45,0.55},{0.45,0.55}}],ListLinePlot[nSol[1000,{pH0->0.51,pP0->0.49},TestPars3],AspectRatio->1,PlotRange->{{0.45,0.55},{0.45,0.55}}]},ImageSize->Full]

Select[Table[simWF[TestPars2][[deme,DirSel[TestPars2,deme],{1,2}]],{deme,51,60}],#≠{}&][[;;,1]]//ListPlot

ListLinePlot[Select[Table[simWF[TestPars2][[deme,DirSel[TestPars2,deme],{1,2}]],{deme,51,60}],#≠{}&],PlotStyle->Directive[Red,Opacity[0.5]]]

ListLinePlot[Select[Table[simWF[TestPars3][[deme,DirSel[TestPars3,deme],{1,2}]],{deme,101,120}],#≠{}&],PlotStyle->Directive[Red,Opacity[0.5]]]

Importing Runs

noMu=NotebookDirectory[]<>"Sims/WF_no_mutation/";
Mu=NotebookDirectory[]<>"Sims/WF_mutation/";

(*{"run","rep","X","Y","alpha","tMax","nDemes","kappa"}*)
Clear[pars]
pars[file_]:=pars[file]=Import[file<>"pars.csv"][[2;;]]

simPars[file_,run_,rep_]:=Block[{temp},temp=pars[file][[7*(run-1)+rep]];
{X->temp[[3]],Y->temp[[4]],α->temp[[5]],κ->temp[[8]](*,μ->temp[[-1]]*)}
]

Clear[repData]
repData[file_,run_,rep_]:=repData[file,run,rep]=Import[file<>"Host_run_"<>ToString[run]<>"_rep_"<>ToString[rep]<>".csv"]

repDataP[file_,run_,rep_]:=repDataP[file,run,rep]=Import[file<>"Path_run_"<>ToString[run]<>"_rep_"<>ToString[rep]<>".csv"]

Processing

(*Both host and pathogen are polymorphic*)
Clear[CoevSel]
CoevSel[file_,run_,rep_,deme_]:=CoevSel[file,run,rep,deme]=Block[{data},
data=Transpose[{repData[file,run,rep][[deme]],repDataP[file,run,rep][[deme]]}];
Quiet[Flatten[Position[data,_?(#[[1]]≠0&&#[[1]]≠(κ/.simPars[file,run,rep])&&#[[2]]≠0&&#[[2]]≠(κ/.simPars[file,run,rep])&)],1]]
]

(*Host Polymorpic, Pathogen Fixed*)
Clear[DirSel]
DirSel[file_,run_,rep_,deme_]:=DirSel[file,run,rep,deme]=Block[{data},
data=Transpose[{repData[file,run,rep][[deme]],repDataP[file,run,rep][[deme]]}];
Quiet[Flatten[Position[data,_?(#[[1]]≠0&&#[[1]]≠(κ/.simPars[file,run,rep])&&(#[[2]]==0||#[[2]]==(κ/.simPars[file,run,rep]))&)],1]]
]

(*Host Fixed*)
Clear[Fixed]
Fixed[file_,run_,rep_,deme_]:=Fixed[file,run,rep,deme]=Block[{data},
data=Transpose[{repData[file,run,rep][[deme]],repDataP[file,run,rep][[deme]]}];
Quiet[Flatten[Position[data,_?((#[[1]]==0||#[[1]]==(κ/.simPars[file,run,rep]))&)],1]]
]

Clear[Hdyn]
Hdyn[file_,run_,rep_]:=Hdyn[file,run,rep]=2((repData[file,run,rep])/(κ/.simPars[file,run,rep]))(1-((repData[file,run,rep])/(κ/.simPars[file,run,rep])))

Hneut[file_,run_,rep_,t_]:=Hneut[κ/.simPars[file,run,rep],t]

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

Plotting

Heterozygosity

HDynPlot[file_,run_,rep_]:=ListLinePlot[Mean[Hdyn[file,run,rep]]];

HDynPlot2[file_,rep_]:=ListLinePlot[Table[Mean[Hdyn[file,run,rep]],{run,1,15(*20*)}],PlotStyle->ColorData["CMYKColors"][((rep-1)/6)],PlotRange->{0,0.5}];

HDynPlot3[file_,rep_]:=Show[HDynPlot2[file,rep],Plot[Hneut[κ/.simPars[file,1,rep],t],{t,0,500},PlotStyle->Directive[Black,Thick],PlotRange->{0,0.5}],PlotRange->All,Frame->True,FrameLabel->{"Time","(Subscript[H, coev])"},LabelStyle->Directive[Black,Bold, Large],FrameTicksStyle->Directive[14],ImageSize->Large]

Box Plot

HetBW[file_,t_]:=Block[{data,sigData,rep,point},sigData={};
point=Graphics[{EdgeForm[Black],White,Disk[Scaled[{0.5,0.5}],0.02]}];
data=Table[Mean[Hdyn[file,run,rep]][[t]]-Hneut[file,run,rep,t],{rep,1,7},{run,1,60}];
(*Which reps are signifcantly below 0*)
For[rep=1,rep≤7,rep++,
If[Mean[BootstrapMean[data[[rep]],10000]]+1.64 StandardDeviation[BootstrapMean[data[[rep]],10000]]<0,
AppendTo[sigData,{{rep,0.1}}],AppendTo[sigData,{{rep,0}}]]
];
Show[DistributionChart[data,ChartStyle->Table[ColorData["CMYKColors"][(rep/6)],{rep,0,6}]],ListPlot[Table[Transpose[{RandomReal[{rep-0.1,rep+0.1},60],data[[rep]]}],{rep,1,7}],PlotMarkers->point],ListPlot[sigData,PlotStyle->Table[ColorData["CMYKColors"][(rep/6)],{rep,0,6}],PlotMarkers->✶],ListLinePlot[{{0.2,0},{8,0}},PlotStyle->Directive[Black]],PlotRange->{{0.5,7.5},{-0.5,0.1}},Frame->True,(*FrameLabel->{"(Subscript[Log, 10])[κ]","ΔH"},*)FrameTicks->{{True,False},{Table[{rep,2+(rep-1)*0.5},{rep,1,7}],False}},LabelStyle->Directive[Black,Bold, 16],FrameTicksStyle->Directive[14],ImageSize->Large]
]

HetBW[noMu,0,500]

Export[NotebookDirectory[]<>"Plots/WF_HetBW_Mu.gif",HetBW[Mu,500],ImageResolution->300];

Mean Fitting Function

Functions

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,((Sqrt[2*(0.5^2)])/(npts+2)),Sqrt[2*(0.5^2)]-((Sqrt[2*(0.5^2)])/(npts+2)),((Sqrt[2*(0.5^2)])/(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)]}]

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

Plotting

Clear[data]
data[file_,t_]:=data[file,t]=Table[{Hneut[file,run,rep,t],Mean[Hdyn[file,run,rep]][[t]]},{rep,1,7},{run,1,60}]

data[noMu,500][[;;,;;,2]]-data[noMu,500][[;;,;;,1]]//BoxWhiskerChart

ΔHrep[file_,run_,rep_,t_]:=Hdyn[file,run,rep][[t]]-Hneut[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=1,rep≤7,rep++,
AppendTo[outSig,{}];AppendTo[outNotSig,{}];
For[run=1,run≤60,run++,
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]],{Hneut[file,run,rep,t],Mean[Hdyn[file,run,rep]][[t]]}],AppendTo[outSig[[-1]],{Hneut[file,run,rep,t],Mean[Hdyn[file,run,rep]][[t]]}]]
];
If[outSig[[rep]]=={},AppendTo[outSig[[rep]],{-0.01,-0.01}]];
If[outNotSig[[rep]]=={},AppendTo[outNotSig[[rep]],{-0.01,-0.01}]];
];
{outNotSig,outSig}
]

moranMean=Import[NotebookDirectory[]<>"MoranMean.csv"];

Clear[ModelComparison]

ModelComparison[file_,t_]:=Show[Plot[x,{x,0,0.5},PlotStyle->Directive[Black,Dashed]],ListLinePlot[moranMean,PlotStyle->Red],ListPlot[sigΔH[file,t][[1]],PlotStyle->Table[Lighter[ColorData["CMYKColors"][(rep/6)]],{rep,0,6}],PlotMarkers->○],ListPlot[sigΔH[file,t][[2]],PlotStyle->Table[Directive[Darker[ColorData["CMYKColors"][(rep/6)]]],{rep,0,6}],PlotMarkers->●],ListLinePlot[RollFitOld[0.1,40,Flatten[data[file,t],1]],PlotStyle->Directive[Darker[Green],Thick]],ListLinePlot[Table[RollFitOld[0.1,40,Bootstrap[Flatten[data[file,t],1],x]],{x,0,1000}],PlotStyle->Directive[Green,Opacity[0.01]]],PlotRange->All,AspectRatio->1,Frame->True,PlotRange->{{0,0.5},{0,0.5}},AspectRatio->1,FrameTicks->{True,True,False,False},FrameTicksStyle->Directive[Black,Bold,14]]

ModelComparison[noMu,250]

Export[NotebookDirectory[]<>"Plots/ModelComparison.pdf",ModelComparison[noMu,250],ImageResolution->300];

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

ModelCompLegend=LineLegend[{{Black,Dashed},Red,Darker[Green]},{"Neutral Expectation","Moran Model","WF Model"}]

Export[NotebookDirectory[]<>"Plots/ModelCompLegend.pdf",ModelCompLegend,ImageResolution->300];

λ Length

λLenPlot[file_,rep_,t_]:=Block[{data,nlm},data=Table[{λLen/.simPars[file,run,rep],Mean[Hdyn[file,run,rep]][[t]]-Hneut[file,run,rep,t]},{run,1,60}];
nlm=NonlinearModelFit[data,{ Exp[-b (x-1)]-a,0<b},{a,b},x];
R2[rep,t]=nlm["RSquared"];
Show[Plot[nlm[x],{x,1,Max[data[[;;,1]]]},PlotRange->All,PlotStyle->ColorData["CMYKColors"][((rep-1)/6)](*,AxesOrigin->{0,-0.5}*)],ListPlot[data,PlotStyle->ColorData["CMYKColors"][((rep-1)/6)]],PlotRange->All(*,Epilog->Text[Style["(R^2)="<>ToString[nlm["RSquared"]],14,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.1,((-0.2)/rep)}]]*)]
]

λLenPlot2[file_,t_]:=Show[Table[λLenPlot[file,r,t],{r,1,7}],PlotRange->{-0.5,0.1},Frame->True,(*FrameLabel->{"||λ||","(Subscript[H, coev])-(Subscript[H, neut])"},*)LabelStyle->Directive[Black,Bold, 16],FrameTicksStyle->Directive[14],ImageSize->{Automatic,300},Epilog->Table[Text[Style[(*"(R^2)="<>*)ToString[Round[R2[rep,500],0.001]],14,ColorData["CMYKColors"][((rep-1)/6)]],Scaled[{0.1,0.1+((0.5-0.1)/7)(8-rep)}]],{rep,1,7}],PlotRange->{-0.5,0.1}]

λLenPlot2[noMu,500]

Export[NotebookDirectory[]<>"Plots/lambdaLenPlotWF.pdf",λLenPlot2[noMu,500]];

Allele Frequency

(*Doesn't work with mutation*)
NPlot[file_,run_,rep_]:=Show[{ListLinePlot[Select[Table[Transpose[{CoevSel[file,run,rep,deme],(1/(κ/.simPars[file,run,rep]))repData[file,run,rep][[deme,CoevSel[file,run,rep,deme]]]}],{deme,1,100}],#≠{}&],PlotStyle->Gray],ListLinePlot[Select[Table[Transpose[{DirSel[file,run,rep,deme],(1/(κ/.simPars[file,run,rep]))repData[file,run,rep][[deme,DirSel[file,run,rep,deme]]]}],{deme,1,100}],#≠{}&],PlotStyle->Red],ListLinePlot[Select[Table[Transpose[{Fixed[file,run,rep,deme],(1/(κ/.simPars[file,run,rep]))repData[file,run,rep][[deme,Fixed[file,run,rep,deme]]]}],{deme,1,100}],#≠{}&],PlotStyle->Blue]},PlotRange->All,Frame->True,(*FrameLabel->{"time","Allele Frequency"},*)LabelStyle->Directive[Black,Bold, 18],FrameTicksStyle->Directive[14],ImageSize->Large]

(*Works with mutation*)
NPlot2[file_,run_,rep_]:=Show[{ListLinePlot[Select[Table[Transpose[{CoevSel[file,run,rep,deme],(1/(κ/.simPars[file,run,rep]))repData[file,run,rep][[deme,CoevSel[file,run,rep,deme]]]}],{deme,1,100}],#≠{}&],PlotStyle->Gray],ListLinePlot[Flatten[Map[Split[#,#2[[1]]-#1[[1]]<2&]&,Select[Table[Transpose[{DirSel[file,run,rep,deme],(1/(κ/.simPars[file,run,rep]))repData[file,run,rep][[deme,DirSel[file,run,rep,deme]]]}],{deme,1,100}],#≠{}&]],1],PlotStyle->Red],ListLinePlot[Select[Table[Transpose[{Fixed[file,run,rep,deme],(1/(κ/.simPars[file,run,rep]))repData[file,run,rep][[deme,Fixed[file,run,rep,deme]]]}],{deme,1,100}],#≠{}&],PlotStyle->Blue]},PlotRange->All,Frame->True,(*FrameLabel->{"time","Allele Frequency"},LabelStyle->Directive[Black,Bold, 18],*)FrameTicksStyle->Directive[Black,Bold,14],ImageSize->{Automatic,300}]

simPars[noMu,1,2]

NPlot2[noMu,1,2]

Export[NotebookDirectory[]<>"Plots/WF_AF.pdf",NPlot2[noMu,1,2],ImageResolution->300];