NotebookDirectory[]

f1AA=p1^2; (* frequency of AA genotype before migration in deme 1 *)
f1AB=2 p1*(1-p1); (* frequency of AB genotype before migration in deme 1 *)
f2AA=p2^2; (* frequency of AA genotype before migration in deme 2 *)
f2AB=2 p2*(1-p2); (* frequency of AB genotype before migration in deme 2 *)
f1AAm=f1AA-m*f1AA+m*f2AA; (* frequency of AA genotype after migration in deme 1 *)
f1ABm=f1AB-m*f1AB+m*f2AB; (* frequency of AB genotype after migration in deme 1 *)
f2AAm=f2AA-m*f2AA+m*f1AA; (* frequency of AA genotype after migration in deme 2 *)
f2ABm=f2AB-m*f2AB+m*f1AB; (* frequency of AB genotype after migration in deme 2 *)
wAA1=1+s; (* fitness of AA in deme 1 *)
wAB1=1+0.5 s; (* fitness of AB in deme 1 *)
wBB1=1; (* fitness of BB in deme 1 *)
wAA2=1; (* fitness of AA in deme 2 *)
wAB2=1+0.5 s; (* fitness of AB in deme 2 *)
wBB2=1+s; (* fitness of BB in deme 2 *)
wbar1=FullSimplify[(f1AAm*wAA1)+(f1ABm*wAB1)+((1-(f1AAm+f1ABm))*wBB1)];
(* average fitness in deme 1 *)
wbar2=FullSimplify[(f2AAm*wAA2)+(f2ABm*wAB2)+((1-(f2AAm+f2ABm))*wBB2)];
(* average fitness in deme 2 *)
p1next=FullSimplify[(f1AAm*wAA1/wbar1)+(0.5*f1ABm*wAB1/wbar1)];
p2next=FullSimplify[(f2AAm*wAA2/wbar2)+(0.5*f2ABm*wAB2/wbar2)];
sols=Solve[{p1next==p1,p2next==p2},{p1,p2}]
p1eq=FullSimplify[(p1/.(sols[[3]]))]
Plot3D[p1eq,{m,0.0001,0.5},{s,0.0001,0.5},PlotRange->{0.5,1},AxesLabel->{{"m","s","p1"},FontSize->14}]

p1eq/.{m->0.05,s->0.02}

mvals={0.01,0.02,0.05,0.1}
toPlot=(p1eq-(1-p1eq))/.m->mvals;
Plot[Evaluate[toPlot],{s,0.001,0.5},PlotLegends->mvals,Frame->True,FrameLabel->{"s","Allele freq. diff."},BaseStyle->20]

toPlot

mvs

ToString[mvals]

effectiveS=FullSimplify[Solve[p1eq==p,s]] (* this is Barton 1983's "s*" *)
Plot[((s/.effectiveS)/.m->0.05),{p,0.51,0.8}]

Plot3D[p1eq,{m,0.0001,0.5},{s,0.0001,1},PlotRange->{0,1},AxesLabel->{"m","s","(Subscript[(OverHat[p]), 1])"},BaseStyle->{18,FontFamily->"Arial Black",FontColor->Black}]

ht=2*0.5*0.5;
hs=2*p1*(1-p1);
fst=FullSimplify[((ht-hs)/ht)/.(p1->p1eq)]
Plot3D[fst,{m,0.00001,0.5},{s,0.00001,1},PlotRange->{0,1},AxesLabel->{"m","s","(Subscript[F, ST])"},BaseStyle->{18,FontFamily->"Arial Black",FontColor->Black}]

W[L_,s_]:=Exp[L*s]
W2[L_,s_]:=(1+s)^L
Exp[Range[0,1,0.1]]
1.1^(Range[0,10])

eqFreq[m_,s_]:=((m (0.5` +0.75` s)-0.125` s-0.125` (Sqrt[(s^2)-4.` m (s^2)+4.` (m^2) (((2.` +1.` s))^2)]))/((-0.25`+1.` m) s))
fitApprox[s_,n_]:=(1+0.5*s)^n

(* number of alleles is binomial with parameters 2L and eqFreq *)
avgFitAllIndep[L_,m_,s_]:=Sum[PDF[BinomialDistribution[2*L,eqFreq[m,s]],i]*fitApprox[s,i],{i,0,2*L}]
maxFit[L_,s_]:=(1+s)^L
minFit[L_,s_]:=1/maxFit[L,s]
avgScaled[L_,m_,s_]:=avgFitAllIndep[L,m,s]/maxFit[L,s]
hetScaled[L_,s_]:=((1+0.5 s)^L)/maxFit[L,s]
avgResident[L_,m_,s_]:=Sum[PDF[BinomialDistribution[2*L,eqFreq[m,s]],i]*fitApprox[s,i],{i,0,2*L}]
avgImmigrant[L_,m_,s_]:=Sum[PDF[BinomialDistribution[2*L,1-eqFreq[m,s]],i]*fitApprox[s,i],{i,0,2*L}]

nLOCI=200;
sVal=0.01;
mVal=0.1;
avgFitVals=Table[avgScaled[L,mVal,sVal],{L,1,nLOCI}];
minFitVals=Table[minFit[L,sVal],{L,1,nLOCI}];
ListPlot[{Transpose[{Range[1,nLOCI],avgFitVals}],Transpose[{Range[1,nLOCI],minFitVals}]}]

nLOCI=200;
sVal=0.02;
mVal=0.1;
avgResVals=Table[avgResident[L,mVal,sVal],{L,1,nLOCI}];
avgImmVals=Table[avgImmigrant[L,mVal,sVal],{L,1,nLOCI}];
ListPlot[{Transpose[{Range[1,nLOCI],avgResVals}],Transpose[{Range[1,nLOCI],avgImmVals}]}]

nLOCI=200;
sVal=0.01;
mVal=0.1;
avgResVals=Table[1,{L,1,nLOCI}];
avgImmVals=Table[avgImmigrant[L,mVal,sVal]/avgResident[L,mVal,sVal],{L,1,nLOCI}];
ListPlot[{Transpose[{Range[1,nLOCI],avgResVals}],Transpose[{Range[1,nLOCI],avgImmVals}]},PlotRange->{0,1}]

nLOCI=200;
sVal=0.02;
mVal=0.1;
maxFitVals=Table[1,{L,1,nLOCI}];
avgResVals=Table[avgResident[L,mVal,sVal]/maxFit[L,sVal],{L,1,nLOCI}];
avgImmVals=Table[avgImmigrant[L,mVal,sVal]/maxFit[L,sVal],{L,1,nLOCI}];

ListLinePlot[{Transpose[{Range[1,nLOCI],maxFitVals}],Transpose[{Range[1,nLOCI],avgResVals}],Transpose[{Range[1,nLOCI],avgImmVals}]},PlotRange->{0,1},Frame->True,FrameLabel->{"Number of loci","Fitness relative to max.","m = "<>ToString[mVal]<>", s = "<>ToString[sVal]},BaseStyle->14,PlotLegends->{"max","avg. resident","avg. immigrant"}]

sVal = 0.02;
LogPlot[((1+sVal)^L)/Exp[0.5 sVal 2 L],{L,1,1000}]
LogPlot[((1+sVal)^L)/((1+0.5*sVal)^2*L),{L,1,1000}]
LogPlot[Exp[0.5 sVal 2 L]/((1+0.5*sVal)^2*L),{L,1,1000}]