Code for ‘Effect of assortative mating and sexual selection on polygenic barriers to gene flow’

Throughout the notebook, s=βα is the selection coefficient at any locus, L is the number of  loci, m is the migration rate , c=((γα)^2), where γ is the strength of assortment and α is the effect size per locus. Thus sL corresponds to the strength of viability selection (Subscript[βΔ, 0]) and (cL^2) to the strength of sexual selection (Subscript[((γΔ)^2), 0]) under conditions of maximum divergence between mainland and island ((Subscript[Δ, 0]) = αL) as in the main text.

Definitions

Analytical predictions for allele frequency at a trait locus

SetDirectory["/Users/mnniba014/Desktop/Projects/AM on RI/filesfromnb"];(*Set the directory to the required folder*)

pavpl::usage="pavpl gives the deleterious allele frequency for Model 1 of assortative mating";
pavpl[m_,s_,l_,c_,p_]:=(m/s)Exp[-2 s l(1-p)-(4/3)c (l^2)(((1-p))^2)];

pavan::usage=" pavan gives the deleterious allele frequency Model 2 of assortative mating";
pavan[m_,s_,l_,c_,p_]:=(m/s)Exp[-2s*l(1-p)]*(∏_(i=0)^100) (1/2)(1+(E^(-(4^(-i)) c*(l^2)*(((1-p))^2))))

findfixed::usage="findfixed returns the equilibrium deleterious allele frequency, by iterating to the fixed point for p starting with inip  ";
findfixed[m_,s_,l_,c_,inip_,pav_]:=Module[{p1,px},
p1=inip;Do[px=Min[pav[m,s,l,c,p1],1];If[Abs[p1-px]<(10^(-10)),Break[]];p1=px,{10000}];
px]

Deterministic simulations

inipop::usage="inipop generates the initial frequency list as frequency of genotype with 0, 1,...l deleterious allele with the wild type frequency 1. ";
inipop[l_]:=Join[{1},ConstantArray[0,l]];
inipop[l_,p_]:=Binomial[l,#](p^(#))(((1-p))^((l-#)))&/@Range[0,l]

W::usage="W gives the fitness of genotypes with 0 to l deleterious alleles with a given selection coefficient s";
W[s_,l_]:= Exp[-s*#]&/@Range[0,l];

pops::usage="pops gives the genoptype frequency after selection";
pops[plist_,s_,l_]:=(plist*W[s,l])/ (plist.W[s,l])

popmig::usage="popmig gives the genptype frequency after migration with a specifed rate m";
popmig[plist_,m_,l_]:=(1-m)*plist +m*Join[ConstantArray[0,l],{1}]

matprobplant::usage="matprobplant gives the probablity of mating with assortment strength c when parent A and parent B has delm and delf no of deleterious alleles";
matingspob::usage="This resturns the mating probablity for every pair of number of deleterious allele combinations";
matprobplant[c_][delm_,delf_]:= Exp[-c (((delm-delf))^2)];
matingspob[c_,l_]:=Array[matprobplant[c][#1-1,#2-1]&,{l+1,l+1}];

recombprob::usage="recombprob calculates the probablity that parents with Y and Z deleterious allele produce an offspring with X deleterious allele for all X";
recombnew[y_,z_,x_,l_]:=Sum[PDF[BinomialDistribution[y+z-2m,0.5],x-m]*PDF[HypergeometricDistribution[z,y,l],m],{m,Max[0,y+z-l],Min[y,z]}];
recombprob[rec_,l_]:=Array[rec[#1,#2,#3,l]&,{l+1,l+1,l+1},0]//N;

fnew[plist_,finaldel_,mat_,rec_]:=Module[{new=mat*rec[[All,All,finaldel+1]]},
plist.new.Transpose[plist]];

norm::usage="norm is the normalizing coefficient in the mating function of the animal model";
norm[c_,plist_,delf_,l_]:=Sum[plist[[i+1]]*Exp[-c(((i-delf))^2)],{i, 0,l}];
normlist[c_,plist_,l_]:=norm[c,plist,#,l]&/@Range[0,l];

matingspobanim::usage="It gives the probability of mating for all combinations of # of deleterious alleles for a given value of c";
matprobanim[delm_,delf_,c_,nlist_]:=Exp[-c(((delm-delf))^2)]/nlist[[delf+1]];
matingspobanim[c_,l_, nlist_]:=Array[matprobanim[#1-1,#2-1, c, nlist]&,{l+1,l+1}];

onegen2::usage="onegen2 gives the genotype frequencies after one round of selection, migraion and mating";
(*Model 1*)
onegen2[l_,m_,s_,c_,mat_,recp_][pop_]:= Module[{popsel,popmigtn, p,pi},
popmigtn=popmig[pop,m,l];
popsel=pops[popmigtn,s,l];
p=fnew[popsel,#,mat,recp]&/@Range[0,l];
pi=p/Total[p];
pi]

(*Model 2*)
onegen2[l_,m_,s_,c_,recp_][pop_]:= Module[{popsel,popmigtn,nlist,mat, p,pi},
popmigtn=popmig[pop,m,l];
popsel=pops[popmigtn,s,l];
nlist=normlist[c,popsel,l];
mat=matingspobanim[c,l,nlist];
p=fnew[popsel,#,mat,recp]&/@Range[0,l];
pi=p/Total[p];
pi]

simall::usage="Performs the dterministic simulation for t generations and returns the deterministic allele frequency";
simall[l_,m_,s_,c_,t_]:=Module[{matall=matingspob[c,l],recall=recombprob[recombnew,l],plist,meandet},
plist=NestList[onegen2[l,m,s,c,matall,recall],inipop[l],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];
simall[l_,m_,s_,c_,t_,mat_,rec_]:=Module[{plist,meandet},
plist=NestList[onegen2[l,m,s,c,mat,rec],inipop[l],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];

simall[inip_,l_,m_,s_,c_,t_]:=Module[{matall=matingspob[c,l],recall=recombprob[recombnew,l],plist,meandet},
plist=NestList[onegen2[l,m,s,c,matall,recall],inipop[l,inip],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];
simall[inip_,l_,m_,s_,c_,t_,mat_,rec_]:=Module[{plist,meandet},
plist=NestList[onegen2[l,m,s,c,mat,rec],inipop[l,inip],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];

simallanim[l_,m_,s_,c_,t_]:=Module[{recall=recombprob[recombnew,l],plist,meandet},
plist=NestList[onegen2[l,m,s,c,recall],inipop[l],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];

simallanim[l_,m_,s_,c_,t_,recall_]:=Module[{plist,meandet},
plist=NestList[onegen2[l,m,s,c,recall],inipop[l],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];

simallanimini[inip_,l_,m_,s_,c_,t_]:=Module[{recall=recombprob[recombnew,l],plist,meandet},
plist=NestList[onegen2[l,m,s,c,recall],inipop[l,inip],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];
simallanimini[inip_,l_,m_,s_,c_,t_,recall_]:=Module[{plist,meandet},
plist=NestList[onegen2[l,m,s,c,recall],inipop[l,inip],t];
Mean[Range[l]*#[[2;;-1]]]&/@plist];

Allele frequency divergence for different strengths of viability and sexual selection

Initializing

l=80;
sLlist={0.2,0.5,1}; (*(Subscript[βΔ, 0]) values*)
cl2list={0,0.5,1,1.5}; (*((Subscript[γΔ, 0])^2) values*)
recall=recombprob[recombnew,l]; 
mbys=Range[0,5,0.2]; (*m/s values*)
mbys2=Range[0,2.5,0.2];

Analytical predictions

pcpl=Table[Table[1-findfixed[#*(sl/l),sl/l,l,cl2/((l^2)),0.001,pavpl]&/@Range[0,5,0.01],{cl2,cl2list}],{sl,sLlist}]; (*Model 1*)
pcan=Table[Table[1-findfixed[#*(sl/l),sl/l,l,cl2/((l^2)),0.001,pavan]&/@Range[0,2.5,0.01],{cl2,cl2list}],{sl,sLlist}];(*Model 2*)

Simulations:
Model1

slcl1=Table[simall[l,#*(0.2/l),0.2/l,i/((l^2)),50000,matingspob[i/(l^2),l],recall]&/@mbys,{i,cl2list}];
slcl2=Table[simall[l,#*(0.5/l),0.5/l,i/((l^2)),50000,matingspob[i/(l^2),l],recall]&/@mbys,{i,cl2list}];
slcl3=Table[simall[l,#*(1/l),1/l,i/((l^2)),50000,matingspob[i/(l^2),l],recall]&/@mbys,{i,cl2list}];

Model 2:

slcl1an=Table[simallanim[l,#*(0.2/l),0.2/l,i/((l^2)),50000,recall]&/@mbys,{i,cl2list}];
slcl2an=Table[simallanim[l,#*(0.5/l),0.5/l,i/((l^2)),50000,recall]&/@mbys,{i,cl2list}];
slcl3an=Table[simallanim[l,#*(1/l),1/l,i/((l^2)),50000,recall]&/@mbys,{i,cl2list}];

Note that these were computed in the cluster

{Show[ListLinePlot[Transpose[{Range[0,5,0.01],#}]&/@pcpl[[1]],PlotStyle->{Gray,Orange,Blue,Red}, PlotLabel->"β(Subscript[(OverBar[z]), m])=0.2", Frame->{True,True,False,False},FrameLabel->{"m/s","Favorable allele frequency"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,5,0.5],Range[0,1,0.2]},FrameStyle->Black,FrameTicksStyle->Black],
ListPlot[Transpose[{mbys,1-#}]&/@slcl1, PlotStyle->{Gray,Orange,Blue,Red}, PlotRange->All]],

Show[ListLinePlot[Transpose[{Range[0,5,0.01],#}]&/@pcpl[[2]],PlotStyle->{Gray,Orange,Blue,Red},PlotLabel->"β(Subscript[(OverBar[z]), m])=0.5",Frame->{True,True,False,False},FrameLabel->{"m/s","Favorable allele frequency"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,5,0.5],Range[0,1,0.2]},FrameStyle->Black,FrameTicksStyle->Black],
ListPlot[Transpose[{mbys,1-#}]&/@slcl2, PlotStyle->{Gray,Orange,Blue,Red}, PlotRange->All]],

Show[ListLinePlot[Transpose[{Range[0,5,0.01],#}]&/@pcpl[[3]],PlotStyle->{Gray,Orange,Blue,Red}, PlotLabel->"β(Subscript[(OverBar[z]), m])=1", Frame->{True,True,False,False},FrameLabel->{"m/s","Favorable allele frequency"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,5,0.5],Range[0,1,0.2]},FrameStyle->Black,FrameTicksStyle->Black,PlotLegends->Placed[PointLegend[Automatic,{0,0.5,1,1.5},LegendLabel->"γ((Subscript[(OverBar[z]), m])^2)"],{0.2,0.3}]],
ListPlot[Transpose[{mbys,1-#}]&/@slcl3, PlotStyle->{Gray,Orange,Blue,Red},PlotRange->All]]}

{Show[ListLinePlot[Transpose[{Range[0,2.5,0.01],#}]&/@pcan[[1]],PlotStyle->{Gray,Orange,Blue,Red},PlotLabel->"β(Subscript[(OverBar[z]), m])=0.2", Frame->{True,True,False,False},FrameLabel->{"m/s","Favorable allele frequency"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.5,0.5],Range[0,1,0.2]},FrameStyle->Black,FrameTicksStyle->Black],
ListPlot[Transpose[{mbys2,1-#}]&/@slcl1an[[All,All,-1]],PlotStyle->{Gray,Orange,Blue,Red}]],
Show[ListLinePlot[Transpose[{Range[0,2.5,0.01],#}]&/@pcan[[2]],PlotStyle->{Gray,Orange,Blue,Red}, PlotLabel->"β(Subscript[(OverBar[z]), m])=0.5",Frame->{True,True,False,False},FrameLabel->{"m/s","Favorable allele frequency"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.5,0.5],Range[0,1,0.2]},FrameStyle->Black,FrameTicksStyle->Black],
ListPlot[Transpose[{mbys2,1-#}]&/@slcl2an[[All,All,-1]],PlotStyle->{Gray,Orange,Blue,Red}]],
Show[ListLinePlot[Transpose[{Range[0,2.5,0.01],#}]&/@pcan[[3]],PlotStyle->{Gray,Orange,Blue,Red},PlotLabel->"β(Subscript[(OverBar[z]), m])=1",PlotLegends->Placed[PointLegend[Automatic,{0,0.5,1,1.5},LegendLabel->"γ((Subscript[(OverBar[z]), m])^2)"],{0.2,0.3}],Frame->{True,True,False,False},FrameLabel->{"m/s","Favorable allele frequency"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.5,0.5],Range[0,1,0.2]},FrameStyle->Black,FrameTicksStyle->Black],
ListPlot[Transpose[{mbys2,1-#}]&/@slcl3an[[All,All,-1]],PlotStyle->{Gray,Orange,Blue,Red}]]}

To find the dependeny on l, change l to required value.

Divergence and critical migration threshold

No assortment

In the absence of assortment, critical divergence  (Subscript[Y, c]) is given by Delcr for any value of (Subscript[βΔ, 0]) and mbyscr gives the critical migration threshold given by (Subscript[m, c])/s =( 1-(Subscript[Y, c]))/ g where g = Exp[-2 (Subscript[βΔ, 0]) g]

Delcr[sl_]:=(2sl-1)/(2sl);

mbyscr[sl_]:=If[Delcr[sl]<0,1,(1-Delcr[sl])/Exp[-2sl Delcr[sl]]]

Model 1

With assortment, critical divergence (Subscript[Y, c]) is Delcrpl for model 1 of assortative mating for a given value of (Subscript[βΔ, 0]) and (Subscript[((γΔ)^2), 0]).  See SM Appendix D for details.

Delcrpl[sl_,cl2_]:=Module[{val=((4cl2-3sl)+Sqrt[9(sl^2)+16 (cl2^2)+24sl cl2-24cl2])/(8cl2)},
If[Element[val,NonNegativeReals],val,0]];

delvalpl=Table[Delcrpl[#,i]&/@Range[0.001,2.1,0.01],{i,{0.5,1,1.5,2}}];

delvalnoc=(val=Delcr[#];If[Element[val,NonNegativeReals],val,0])&/@Range[0.001,2.1,0.01];

delvalplfinal=Join[{delvalnoc},delvalpl];

ListLinePlot[Transpose[{Range[0.001,2.1,0.01],#}]&/@delvalplfinal[[1;;-2]],PlotStyle->{Gray,Orange,Blue,Red,Purple},Frame->{True,True,False,False},FrameLabel->{"(Subscript[βΔ, max])","(Subscript[Y, c])"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.2,0.5],Range[0,1,0.2]},FrameStyle->Black, PlotRange->{{0,2.2},{0,1}}]

mepl gives the gene flow factor corresponding to a divergence level Δ.  mcrbyspl gives the critical migration threshold for given values of (Subscript[βΔ, 0]) and (Subscript[((γΔ)^2), 0])

mepl[Δ_,sl_,cl2_]:=Exp[-2sl Δ- (4/3)cl2 (Δ^2)];

mcrbyspl[sl_,cl2_]:=If[Delcrpl[sl,cl2]==0,1,Max[1,(1-Delcrpl[sl,cl2])/mepl[Delcrpl[sl,cl2],sl,cl2]]];

vals=Table[mcrbyspl[#,cl2]&/@Range[0.001,2.1,0.1],{cl2,{0.5,1,1.5,2}}];

newvals=Join[{mbyscr[#]&/@Range[0.001,2.1,0.1]},vals];

ListLinePlot[Transpose[{Range[0.001,2.1,0.1],#}]&/@newvals[[1;;-2]], PlotStyle->{Gray,Orange,Blue,Red,Purple},Frame->{True,True,False,False},FrameLabel->{"(Subscript[βΔ, max])","(Subscript[m, c])/s"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.5,0.5],Range[0,20,5]},FrameStyle->Black, PlotRange->{{0,2.2},{0,20}}]

Finding minimum strength of sexual selection required to get a tipping point associated with (Subscript[Y, c])>0 and to shift the critical migration threshold beyond 1, (Subscript[m,  c])/s >1 for given value of (Subscript[βΔ, 0]).
getclanl returns the minimum (Subscript[((γΔ)^2), 0]) needed to obtain the critical divergence threshold.

getclanal[sl_]:=If[sl<=0.5,(3/4)(1+Sqrt[1-2*sl]-sl),0]

To calculate minimum (Subscript[((γΔ)^2), 0]) for shifted critical migration thresholds beyond 1  one needs to solve ( 2(Subscript[βΔ, 0])+(8/3)(Subscript[((γΔ)^2), 0]) (Subscript[Δ, c])-1)mepl-(Subscript[Δ, c])]. We simplify this equation , given by solpl1 and solve this using the function clcric, for given values of (Subscript[βΔ, 0]).

FullSimplify[(2sl+(8/3)cl2  del[sl,cl2]-1)mepl[del[sl,cl2],sl,cl2]-del[sl,cl2]]

solpl1[sl_,cl2_]:=-((4 cl2-3 sl+(Sqrt[8 cl2 (-3+2 cl2)+24 cl2 sl+9 (sl^2)]))/(8 cl2))+(1/3) (E^(-(((4 cl2-3 sl+(Sqrt[8 cl2 (-3+2 cl2)+24 cl2 sl+9 (sl^2)])) (4 cl2+9 sl+(Sqrt[8 cl2 (-3+2 cl2)+24 cl2 sl+9 (sl^2)])))/(48 cl2)))) (-3+4 cl2+3 sl+(Sqrt[8 cl2 (-3+2 cl2)+24 cl2 sl+9 (sl^2)]))

clcric[sl_]:=Module[{sol},
sol=x/.FindRoot[solpl1[sl,x],{x,1.6}];
If[sol>=0,sol,0]]

clcricfinal=clcric[#]&/@Range[0,1,0.01];

Show[ListLinePlot[Transpose[{Range[0,1,0.01],clcricfinal}], AxesLabel->{"sL","(cL^2)"}, PlotStyle->Black, PlotLegends->{"((Subscript[m, c])/s)=1"}],ListLinePlot[Transpose[{#,getclanal[#]}]&/@Range[0,1,0.01], PlotStyle->Gray, PlotLegends->{"(Subscript[Δ, c])>0"}]]

Model 2:

For model 2 of assortative mating, critical divergence (Subscript[Y, c]) is given by the function delcricanim forgiven strength of viability selection (Subscript[βΔ, 0]) ,sexual selection (Subscript[((γΔ)^2), 0])  and divergence Δ.

DLnganimal[sl_,cl2_,Δ_]:=-2 sl - Sum[(((2^(1-2 i))  Δ  cl2 )/(1+(E^((4^(-i)) (Δ^2) cl2)))),{i,0,20}]

delfn[sl_,cl2_,Δ_]:=1+ 1/DLnganimal[sl,cl2,Δ]

delcricanim::usage="Returns the critical divergence for the animal model for given strength of viability and sexual selection and iterating until the fixed point starting with an initial guess inidel.";
delcricanim[sl_,cl2_,inidel_]:=Module[{del1,px},
del1=inidel;Do[px=Max[delfn[sl,cl2,del1],0];If[Abs[del1-px]<(10^(-10)),Break[]];del1=px,{10000}];
px]

delcran=Table[delcricanim[#,i,0.5]&/@Range[0.001,2.1,0.01],{i,{0.5,1,1.5,2}}];

delcranfin=Join[{delvalnoc},delcran];

ListLinePlot[Transpose[{Range[0.001,2.1,0.01],#}]&/@delcranfin[[1;;-2]],PlotStyle->{Gray,Orange,Blue,Red,Purple},Frame->{True,True,False,False},FrameLabel->{"(Subscript[βΔ, max])","(Subscript[Y, c])"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.2,0.5],Range[0,1,0.2]},PlotRange->{{0,2.2},{0,1}},FrameStyle->Black,PlotLegends->Placed[PointLegend[{0,0.5,1,1.5},LegendLabel->"((Subscript[γΔ, max])^2)"],{0.8,0.3}]]

mean gives the gene flow factor corresponding to a divergence level Δ for model 2 of assortative mating.  mcrbysan gives the critical migration threshold for given values of (Subscript[βΔ, 0]) and (Subscript[((γΔ)^2), 0]).

mean[sl_,cl2_,Δ_]:=Exp[-2 sl Δ ]*(∏_(i=0)^20) (1/2)(1+(E^(-(4^(-i)) cl2* (Δ^2))))

mcrbysan[sl_,cl2_,inip_]:=If[delcricanim[sl,cl2,inip]==0,1,Max[1,(1-delcricanim[sl,cl2,inip])/mean[sl,cl2,delcricanim[sl,cl2,inip]]]];

mcrsan1=Table[mcrbysan[#,cl2,0.99]&/@Range[0.001,2.1,0.1],{cl2,{0.5,1,1.5,2}}];

newvalsani=Join[{mbyscr[#]&/@Range[0.001,2.1,0.1]},mcrsan1];

ListLinePlot[Transpose[{Range[0.001,2.1,0.1],#}]&/@newvalsani[[1;;-2]],PlotStyle->{Gray,Orange,Blue,Red,Purple},Frame->{True,True,False,False},FrameLabel->{"(Subscript[βΔ, max])","(Subscript[m, c])/s"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,2.5,0.5],Range[0,20,5]},FrameStyle->Black, PlotRange->{{0,2.2},{0,20}},PlotLegends->Placed[PointLegend[{0,0.2,1,1.5},LegendLabel->"((Subscript[γΔ, max])^2)"],{0.2,0.6}]]

getclboth returns the minimum strength of sexual selection  (Subscript[((γΔ)^2), 0]) required to obtain the critical divergence threshold and critical migration thresholds greater than 1 for a given strength of viability selection  (Subscript[βΔ, 0]) .

getclboth[sl_,delfn_,inip_]:=Module[{del,mval,posdel,posm},
del=delfn[sl,#,inip]&/@Range[0,10,0.01];
mval=Table[(1-del[[i]])/mean[sl,Range[0,10,0.01][[i]],del[[i]]],{i,Length[del]}];
posdel=Flatten[Position[del, _?(#> 0&)]];
posm=Flatten[Position[mval, _?(#> 1&)]];
If[Length[posdel]==0,Indeterminate,{sl,del[[posdel[[1]]]],mval[[posdel[[1]]]],mval[[posm[[1]]]],Range[0,10,0.01][[posdel[[1]]]],Range[0,10,0.01][[posm[[1]]]]}]];

clcranimall=getclboth[#,delcricanim,0.5]&/@Range[0.001,1,0.01];

Show[ListLinePlot[DeleteCases[clcranimall,Indeterminate][[All,{1,-2}]] ,PlotStyle->Darker[Gray],Frame->{True,True,False,False},FrameLabel->{"(Subscript[βΔ, max])","((Subscript[γΔ, max])^2)"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameTicks->{Range[0,1,0.1],Range[0,10,1]},FrameStyle->Black,PlotLegends->Placed[{"(Subscript[Y, c])>0"},{0.8,0.45}]],
ListLinePlot[DeleteCases[clcranimall,Indeterminate][[All,{1,-1}]], PlotStyle->Black,LabelStyle->Directive[Black,FontFamily->"Times"],PlotLegends->Placed[{"(Subscript[m, c])/s >1"},{0.8,0.35}]],PlotRange-> All]

Visualizing behaviour of equilibria

f::usage="f refers to the function that satisfies f[Y,c,s,L,m]= 0 at migration- selection- assortment equilibrium as in equation xx in SM";
f[Y_,c_,s_,l_,m_]:=Y*(1-Y)-(m/s)Y mepl[Y,c,s,l]

l=80;s=0.2/80;c=0.5/(80^2);m={0.5,0.7,0.9,1,1.4}*s;fsc=f[Δ,c,s,l,#]&/@m;
Plot[fsc,{Δ,0,1}, PlotLabel->"sL=0.2, (cl^2)=0.5",AxesLabel->{"Δ","f(Δ)"},PlotLegends->{"m/s= 0.5","m/s= 0.7","m/s= 0.9","m/s= 1","m/s= 1.4"}, PlotStyle->{Black, Brown, Gray,Blue, Purple},AxesStyle->Black,Ticks->{Range[0,1,0.1],Range[-0.5,0.1,0.1]}]

s=0.2/80;c=1.3/(80^2);m={0.9,0.93,0.949,0.96,1}*s;
fsc=f[Δ,c,s,l,#]&/@m;
Plot[fsc,{Δ,0,1}, PlotLabel->"sL=0.2, (cl^2)=1.3",AxesLabel->{"Δ","f(Δ)"},PlotLegends->{"m/s= 0.9","m/s= 0.93","m/s= 0.949","m/s= 0.96","m/s= 1","m/s=1.1"}, PlotStyle->{Black, Brown, Blue, Purple,Red},AxesStyle->Black,PlotRange->{{0,1},{-0.05,0.015}},Ticks->{Range[0,1,0.1],Range[-0.05,0.01,0.01]}]

s=0.2/80;c=1.6/(80^2);m={0.8,0.9,1,1.129,1.4}*s;
fsc=f[Δ,c,s,l,#]&/@m;
Plot[fsc,{Δ,0,1}, PlotLabel->"sL=0.2, (cl^2)=1.6",AxesLabel->{"Δ","f(Δ)"},PlotLegends->{"m/s= 0.8","m/s= 0.9","m/s= 1","m/s= 1.129","m/s= 1.4"}, PlotStyle->{Black, Gray,Blue,Brown, Purple},AxesStyle->Black,Ticks->{Range[0,1,0.1],Range[-0.1,0.1,0.05]}]

Effect of assortative mating and viability selection on (Subscript[F, ST]) at an unlinked neutral locus

Let Fst_0 =1/ (1+ 2Nm) =0.02, is the (Subscript[F, ST]) in the absence of any genetic barriers between the mainland and island . This gives 2Nm = 49. Thus, Fst= 1/(1+ 2Nm_e) = 1/ (1+ 2Nmg) = 1/(1+49g)

delhigh::usage="For a given value of (Subscript[βΔ, 0]), delhigh retruns the position from a range of (Subscript[((γΔ)^2), 0]) where divergence crosses 0.95 for model 1 of assortative mating.
delhighan does the same for model 2 of assortative mating";
delhigh[sl_]:=First@Flatten[Position[1-findfixed[1*(sl/80),sl/80,80,#/(80^2),0.001,pavpl]&/@Range[0,2,0.01],_?(#>=0.95&)]];
delhighan[sl_]:=First@Flatten[Position[1-findfixed[1*(sl/80),sl/80,80,#/(80^2),0.001,pavan]&/@Range[0,2,0.01],_?(#>=0.95&)]];

pfuncpl:: usage="returns the divergence level/advatnageous allel frequency for a given value of (Subscript[βΔ, 0]) and (Subscript[((γΔ)^2), 0]) for rate of migration m for model 1 of assortative mating.
pfuncan does the same for model 2 of assortative mating";
pfuncpl[m_,sl_,cl2_,l_]:=1-findfixed[m*(sl/l),sl/l,l,cl2/((l^2)),0.001,pavpl];
pfuncan[m_,sl_,cl2_,l_]:=1-findfixed[m*(sl/l),sl/l,l,cl2/((l^2)),0.001,pavan];

Fstpl::usage=" Fstpl returns Fst = 1/(1+ 2Nmg)= 1/(1+49g) for an unlinked neutral locus when Fst_0 =1/ (1+ 2Nm) =0.02 for model 2 of assortative mating . 
Fstani does the same for model 2 of assortative mating. Refer to the main text for details.";
Fstpl[mbys_,sl_,cl2_,nm_]:= 1/(1+nm mepl[pfuncpl[mbys,sl,cl2,l],sl,cl2]);
Fstani[mbys_,sl_,cl2_,nm_]:=1/(1+nm mean[sl,cl2,pfuncan[mbys,sl,cl2,l]]);

posdel returns the position from a list of (Subscript[((γΔ)^2), 0]) values and delhighall and delhighan retruns both the (Subscript[βΔ, 0]) and (Subscript[((γΔ)^2), 0])  corresponding to when divergence corsses 0.95
delnonzeroan and delnonzeropl returns a list of (Subscript[βΔ, 0]) and (Subscript[((γΔ)^2), 0])  values where for a given (Subscript[βΔ, 0]), the minimum (Subscript[((γΔ)^2), 0]) for which (Subscript[m, c])/s =1 (i.e divergence =0 )

posdel=delhigh[#]&/@Range[0.32,2,0.01];

posdelan=delhighan[#]&/@Range[1.15,2,0.01];

delhighall=Transpose[{Range[0.32,2,0.01],Range[0,2,0.01][[#]]&/@posdel}];

delhighallan=Transpose[{Range[1.15,2,0.01],Range[0,2,0.01][[#]]&/@posdelan}];

delnonzeroan=Join[DeleteCases[clcranimall,Indeterminate][[All,{1,-2}]],Transpose[{Range[1.001,2,0.01],ConstantArray[0,100]}]] ;

delnonzeropl={#,getclanal[#]}&/@Range[0.001,2,0.01];

l=80;

DensityPlot[Fstpl[1,x,y,49],{x,0.0001,2},{y,0,2},PlotLayout->"Row", Epilog->{Thickness[0.005],Line[delnonzeropl],Dashed,Line[delhighall]}, PlotLegends->Automatic,FrameLabel->{"(Subscript[βΔ, 0])","((Subscript[γΔ, 0])^2)"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameStyle->Black,FrameTicksStyle->Black, ColorFunction->"TemperatureMap",Mesh->{10,10}]

DensityPlot[Fstani[1,x,y,49],{x,0.0001,2},{y,0,2}, Epilog->{Thickness[0.005],Line[delnonzeroan[[35;;-1]]],Dashed,Line[delhighallan]}, PlotLegends->Automatic,FrameLabel->{"(Subscript[βΔ, 0])","((Subscript[γΔ, 0])^2)"},LabelStyle->Directive[Black,FontFamily->"Times"],FrameStyle->Black,FrameTicksStyle->Black, ColorFunction->"TemperatureMap",Mesh->{10,10}]