Fst between haploids and diploids in species with discrete ploidal generations

Here we put supplementary information for our manuscript, “Fst between haploids and diploids in species with discrete ploidal generations”.
Note that, you might have to reset the variables as necessary because we reuse same variables in the different sections.

Appendix

Appendix A

Equations

We define the probabilities of identities Qij of genes sampled in deme i and j (QHD=QDH).
When we derive the recursions for Q, we assume that the number of haploid and diploid individuals are fixed as NH=N (Subscript[ρ, H]) and ND=N (Subscript[ρ, D]). We also assume that the probability of identity from diploid population is independent of whether the two alleles are sampled from the same individual or different individuals (no inbreeding).

NextQHH=(((1-μ))^2)*((mHfromH^2)*((1/Hap)+(1-(1/Hap))*QHH)+mHfromH*mHfromD*QHD+mHfromD*mHfromH*QDH+(mHfromD^2)*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD));
NextQHD=(((1-μ))^2)*(mHfromH*mDfromH*((1/Hap)+(1-(1/Hap))*QHH)+mHfromH*mDfromD*QHD+mHfromD*mDfromH*QDH+mHfromD*mDfromD*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD));
NextQDH=(((1-μ))^2)*(mDfromH*mHfromH*((1/Hap)+(1-(1/Hap))*QHH)+mDfromH*mHfromD*QHD+mDfromD*mHfromH*QDH+mDfromD*mHfromD*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD));
NextQDD=(((1-μ))^2)*((mDfromD^2)*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD)+mDfromH*mDfromD*QHD+mDfromD*mDfromH*QDH+(mDfromH^2)*((1/Hap)+(1-(1/Hap))*QHH));

The recursion equations for Q can be described as, (((1-μ))^2)A (Q+c);
(this equation corresponds to Eq. (9.30) in Rousset (2004))

MatrixA=(({{mHfromH*mHfromH, mHfromH*mHfromD, mHfromD*mHfromH, mHfromD*mHfromD}, {mHfromH*mDfromH, mHfromH*mDfromD, mHfromD*mDfromH, mHfromD*mDfromD}, {mDfromH*mHfromH, mDfromH*mHfromD, mDfromD*mHfromH, mDfromD*mHfromD}, {mDfromH*mDfromH, mDfromH*mDfromD, mDfromD*mDfromH, mDfromD*mDfromD}}));
VectC=(({{(1/Hap)*(1-QHH)}, {0}, {0}, {(1/(2*Dip))*(1-QDD)}}));
(({{NextQHH}, {NextQHD}, {NextQDH}, {NextQDD}}))-(((1-μ))^2)*MatrixA.((({{QHH}, {QHD}, {QDH}, {QDD}}))+VectC)//Simplify

Note that, we can write, A = TensorProduct[F, F]

MatrixF=(({{mHfromH, mHfromD}, {mDfromH, mDfromD}}));
TensorProduct[MatrixF,MatrixF]//MatrixForm
MatrixA//MatrixForm

Note that, the matrix A is a probability matrix;

MatrixA.(({{1}, {1}, {1}, {1}}))/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD//Simplify

Furthermore, this can be represented as, (((1-μ))^2)(GQ+Aδ) and (((1-μ))^2)(GQ+(I-G)1);
(this equation corresponds to Eq. (4.4) in Rousset (2004))

MatrixG=(({{(mHfromH^2)*(1-(1/Hap)), mHfromH*mHfromD, mHfromD*mHfromH, (mHfromD^2)*(1-(1/(2*Dip)))}, {mHfromH*mDfromH*(1-(1/Hap)), mHfromH*mDfromD, mHfromD*mDfromH, mHfromD*mDfromD*(1-(1/(2*Dip)))}, {mDfromH*mHfromH*(1-(1/Hap)), mDfromH*mHfromD, mDfromD*mHfromH, mDfromD*mHfromD*(1-(1/(2*Dip)))}, {(mDfromH^2)*(1-(1/Hap)), mDfromH*mDfromD, mDfromD*mDfromH, (mDfromD^2)*(1-(1/(2*Dip)))}}));
Vectδ=(({{(1/Hap)}, {0}, {0}, {(1/(2*Dip))}}));
(({{NextQHH}, {NextQHD}, {NextQDH}, {NextQDD}}))-(((1-μ))^2)*(MatrixG.(({{QHH}, {QHD}, {QDH}, {QDD}}))+MatrixA.Vectδ)//Simplify
(({{NextQHH}, {NextQHD}, {NextQDH}, {NextQDD}}))-(((1-μ))^2)*(MatrixG.(({{QHH}, {QHD}, {QDH}, {QDD}}))+((({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))-MatrixG).(({{1}, {1}, {1}, {1}})))/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD//Simplify

Exact Fst

Recall that these are backwards rates so that mDD+mHD sum to one:

constraints={mHfromH+mHfromD==1,mDfromH+mDfromD==1,ρD+ρH==1};
subconstraints={mHfromH->1-mHfromD,mDfromD->1-mDfromH};

At equilibrium, we have

sol=Simplify[Solve[(NextQHH==QHH&&NextQDD==QDD&&NextQHD==QHD&&NextQDH==QDH),{QHH,QHD,QDH,QDD}],constraints];

When there is not mutation, because of a finite population size, Q must be 1.

sol/.μ->0//Simplify

EqQHH=Simplify[QHH/.sol[[1]],constraints]
EqQHD=Simplify[QHD/.sol[[1]],constraints]
EqQDH=Simplify[QDH/.sol[[1]],constraints]
EqQDD=Simplify[QDD/.sol[[1]],constraints]

Note that, we can show QHD = QDH, as expected:

EqQHD==EqQDH

Using Q, we define the F-like measures for the haploid-diploid population.

FaiHHstrict=Simplify[Limit[((EqQHH-EqQHD)/(1-EqQHD))/.Dip->ρD*Ntot/.Hap->ρH*Ntot,μ->0],constraints]
FaiDDstrict=Simplify[Limit[((EqQDD-EqQHD)/(1-EqQHD))/.Dip->ρD*Ntot/.Hap->ρH*Ntot,μ->0],constraints]

Through a series of collecting factors, as in the following command, and gathering terms

Collect[Numerator[FaiHHstrict],{Ntot},Factor]

we can rewrite FaiHHstrict and FaiDDstrict as:

Eq4a=(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*(2*ρD)-mHfromD*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2*ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH));
%-FaiHHstrict/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH//Factor

Eq4b=(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH));
%-FaiDDstrict/.X->mDfromH+mHfromD/.subconstraints/.ρD->1-ρH//Factor

where X = mDH+mHD, the sum of the movement rates.

Fexact={(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*(2*ρD)-mHfromD*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2*ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH)),(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH))};

If the movement rates are neither near zero nor near one and the population size is large

If N is large (total population size) relative to the migration rates (i.e., none of the migration rates is very rare) these reduce to [dropping terms not involving Ntot] for {haploids, diploids}:

Series[Fexact/.Ntot->(Ntot/ϵ),{ϵ,0,1}]
Factor[Normal[Series[Fexact/.Ntot->(Ntot/ϵ),{ϵ,0,1}]]]/.ϵ->1

FNlarge={(((1-X)*((2-X)*X*(2*ρD)-mHfromD(ρH+2*ρD)))/((2-X)*(X^2) *ρH*(2*ρD) *Ntot)),(( (1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((2-X)*(X^2) *ρH*(2*ρD) *Ntot))};
Simplify[%-%%,constraints]

We can use the above general Fst result when Ntot is large to estimate the degree of sexuality, although this is too complicated to present in the main paper:

FNlarge/.X->mHfromD+mDfromH/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD/.mHfromH->((θHH)/Ntot)/.mDfromD->((θDD)/Ntot)//Simplify;
Simplify[Solve[{%[[1]]==ϕH,%[[2]]==ϕD},{θHH,θDD}],constraints]

In the primarily sexual case (mDD and mHH small):

Factor[Normal[Series[Fexact/.X->mHfromD+mDfromH/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD/.mHfromH->((θHH)/Ntot)/.mDfromD->((θDD)/Ntot)/.Ntot->(Ntot/ϵ),{ϵ,0,0}]]]

%[[1]]-((ρH+(2*ρD))/(ρH+(2*ρD)+4*ρH*(2*ρD)*(θHH+θDD)))//Simplify

Consistent with (B.1).

In the fully symmetric case with ρH=2/3:

Factor[Fexact/.X->mHfromD+mDfromH/.mHfromH->1-m/.mDfromD->1-m/.mHfromD->m/.mDfromH->m/.ρD->1-ρH/.ρH->(2/3)]//Simplify

In agreement with equation (6):

((((1-2*m))^2)/((2/3)*Ntot-(((1-2*m))^2)((2/3)*Ntot-1)));
%%-%//Factor

Extreme values of Fst

Simplify[FaiHHstrict/.subconstraints/.mHfromD->0/.mDfromH->0,constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->0/.mDfromH->0,constraints]

Simplify[FaiHHstrict/.subconstraints/.mHfromD->(1/2)/.mDfromH->(1/2),constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->(1/2)/.mDfromH->(1/2),constraints]

Simplify[FaiHHstrict/.subconstraints/.mHfromD->1/.mDfromH->1,constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->1/.mDfromH->1,constraints]

Simplify[FaiHHstrict/.subconstraints/.mHfromD->0/.mDfromH->1,constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->0/.mDfromH->1,constraints]

Simplify[FaiHHstrict/.subconstraints/.mHfromD->1/.mDfromH->0,constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->1/.mDfromH->0,constraints]

Simplify[FaiHHstrict/.subconstraints/.mHfromD->1-mDfromH,constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->1-mDfromH,constraints]

Appendix B

Equations

FaiHHstrict=(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*(2*ρD)-mHfromD*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2*ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH))/.X->mHfromD+mDfromH//Simplify;
FaiDDstrict=(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH))/.X->mHfromD+mDfromH//Simplify;

constraints={mHfromH+mHfromD==1,mDfromH+mDfromD==1,ρD+ρH==1};
subconstraints={mHfromH->1-mHfromD,mDfromD->1-mDfromH};

Complete symmetry case

Simplify[FaiHHstrict/.subconstraints/.mHfromD->m/.mDfromH->m/.ρD->1-ρH/.ρH->(2/3),constraints]
Simplify[FaiDDstrict/.subconstraints/.mHfromD->m/.mDfromH->m/.ρD->1-ρH/.ρH->(2/3),constraints]

Eq6=((((1-2*m))^2)/((2/3)*Ntot-(((1-2*m))^2)*((2/3)*Ntot-1)));
Eq6-%%//Simplify

If we assume N large and N m of order one, we get

Simplify[Normal[Series[Eq6/.m->(θ/Ntot)/.Ntot->(Ntot/ϵ),{ϵ,0,0}]]]

which reduces to the classic Fst prediction when expressing the total population size in terms of the effective number of diploids:

Simplify[%/.θ->(3/2)*m*cNtot]

Limit[Eq6,m->0]
Limit[Eq6,m->(1/2)]
Limit[Eq6,m->1]

Rare sexuality in a large population

We consider the case when there is frequent asexual reproduction (a≒1, and then mHfromD and mDfromH is small), and large population size. We assume;
・mHfromD = O(ϵ)
・mDfromH = O(ϵ)
・Ntot = O(1/ϵ)
These assumptions indicates, θHD = mHfromD Ntot = O(1) and θDH = mDfromH Ntot = O(1).

FNlarge[[1]]/.X->mHfromD+mDfromH

FaiHH◆RS=Simplify[Normal[Series[FaiHHstrict/.subconstraints/.ρD->1-ρH/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

FaiDD◆RS=Simplify[Normal[Series[FaiDDstrict/.subconstraints/.ρD->1-ρH/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

Eq9a=((1+(2*ρD)*(θHD+θDH)-(ρH*θHD-(2*ρD)*θDH))/(1+2*(θHD+θDH)+2*ρH*(2*ρD)*(((θHD+θDH))^2)-(ρH-(2*ρD))*θDH));
Eq9b=((1+ρH*(θHD+θDH)+(ρH*θHD-(2*ρD)*θDH))/(1+2*(θHD+θDH)+2*ρH*(2*ρD)*(((θHD+θDH))^2)-(ρH-(2*ρD))*θDH));
FaiHH◆RS-Eq9a/.ρD->1-ρH//Simplify
FaiDD◆RS-Eq9b/.ρD->1-ρH//Simplify

FaiHH◆RS/.ρH->(2/3)/.θHD->m*Ntot/.θDH->m*Ntot/.Ntot->(3/2)*Nlocal//Factor
FaiDD◆RS/.ρH->(2/3)/.θHD->m*Ntot/.θDH->m*Ntot/.Ntot->(3/2)*Nlocal//Factor

sol=Simplify[Solve[(Eq9a/.ρD->1-ρH)==φH&&(Eq9b/.ρD->1-ρH)==φD,{θHD,θDH}],0<ρH<1]

Eq10a=(((2*ρD)*(1-φH)*(φH+2*φD)-ρH*φH*(1-φD))/(2*ρH*(2*ρD)*(((φH+φD))^2)));
Eq10b=((ρH*(1-φD)*(2*φH+φD)-(2ρD)*φD*(1-φH) )/(2*ρH*(2*ρD)*(((φH+φD))^2)));
(θHD/.sol[[1]])-Eq10a/.ρD->1-ρH//Simplify
(θDH/.sol[[1]])-Eq10b/.ρD->1-ρH//Simplify

This matches equation 10.

Common sexuality in a large population

We consider the case when the species has frequent sexual reproduction (a≒0, and then mHH and mDD is small), large population size.
We assume mHH = O(ϵ), mDD = O(ϵ), N = O(1/ϵ), mHH N = O(1), and mDD N = O(1).

FaiHH◆RA=Simplify[Normal[Series[FaiHHstrict/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD/.mHfromH->((θHH)/Ntot)/.mDfromD->((θDD)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

FaiDD◆RA=Simplify[Normal[Series[FaiDDstrict/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD/.mHfromH->((θHH)/Ntot)/.mDfromD->((θDD)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

These are equal to one another.

EqB1=((ρH+(2*ρD))/(ρH+(2*ρD)+4* ρH*(2*ρD)*(θHH+θDD)));
FaiHH◆RA-EqB1/.ρD->1-ρH//Simplify
FaiDD◆RA-EqB1/.ρD->1-ρH//Simplify

sol=Simplify[Solve[((ρH+(2*ρD))/(ρH+(2*ρD)+4* ρH*(2*ρD)*sum))==φ,{sum}],0<ρH<1]

(sum/.sol[[1]])-(((1-φ)*(ρH+(2*ρD)))/(4*φ*ρH*(2*ρD)))/.ρD->1-ρH//Simplify

which matches equation (B.2).

Rare sexuality in a very large population

We consider the case when the specie represents frequent asexual reproduction (a≒1, and then mDtoH and mHtoD is small), large population size.
We assume mDtoH = O(ϵ), mHtoD = O(ϵ), N = O(1/(ϵ^2)), mDtoH N = O(1/ϵ), and mHtoD N = O(1/ϵ).

Series[FaiHHstrict/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH/.mHfromD->cmHfromD*ϵ/.mDfromH->cmDfromH*ϵ/.Ntot->cNtot*(1/(ϵ^2)),{ϵ,0,1}]//FullSimplify
FaiHH◆RSlargemN=Normal[%]/.cmDfromH->(mDfromH/ϵ)/.cmHfromD->(mHfromD/ϵ)/.cNtot->Ntot*(ϵ^2)/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)//FullSimplify

Series[FaiDDstrict/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH/.mHfromD->cmHfromD*ϵ/.mDfromH->cmDfromH*ϵ/.Ntot->cNtot*(1/(ϵ^2)),{ϵ,0,1}]//FullSimplify
FaiDD◆RSlargemN=Normal[%]/.cmDfromH->(mDfromH/ϵ)/.cmHfromD->(mHfromD/ϵ)/.cNtot->Ntot*(ϵ^2)/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)//FullSimplify

FaiHH◆RSlargemN-(((2*θDH+θHD)*(2*ρD)-θHD*ρH)/(4*(((θDH+θHD))^2)*ρH*ρD))/.ρD->1-ρH//Simplify
FaiDD◆RSlargemN-(((2*θHD+θDH)*ρH-θDH*(2*ρD))/(4*(((θDH+θHD))^2)*ρH*ρD))/.ρD->1-ρH//Simplify

sol=Simplify[Solve[(FaiHH◆RSlargemN/.ρD->1-ρH)==φH&&(FaiDD◆RSlargemN/.ρD->1-ρH)==φD,{θHD,θDH}],0<ρH<1]

θHD+θDH/.sol[[1]]//Simplify

We note, this approximation can be derived from the approximation of large N limit.

Simplify[Normal[Series[FNlarge/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]
%-{FaiHH◆RSlargemN,FaiDD◆RSlargemN}

Appendix C

Definition of Fst in Whitlock and Barton (1997)

The generalized Wright’s Fst coefficient can be defined as;
Corr[ij]=((Fij - Fave)/(1-Fave))
Because Fij could be approximated as;
Fij(t)=1-((rijλ)^t)
Fave = ((H^2)/(((H+2D))^2))FHH+((H(2D))/(((H+2D))^2))FHD+(((2D)H)/(((H+2D))^2))FDH+((((2D))^2)/(((H+2D))^2))FHH
using leading eigenvalue λ and right eigenvector u for matrix G, the Fst becomes
 Corrijave[ij]=((Fij - Fave)/(1-Fave))=(((((H^2)/(((H+2D))^2))rHH+((H(2D))/(((H+2D))^2))rHD+(((2D)H)/(((H+2D))^2))rDH+((((2D))^2)/(((H+2D))^2))rHH-rij)(λ^t))/((((H^2)/(((H+2D))^2))rHH+((H(2D))/(((H+2D))^2))rHD+(((2D)H)/(((H+2D))^2))rDH+((((2D))^2)/(((H+2D))^2))rHH)(λ^t)))=1-(rij/(((H^2)/(((H+2D))^2))rHH+((H(2D))/(((H+2D))^2))rHD+(((2D)H)/(((H+2D))^2))rDH+((((2D))^2)/(((H+2D))^2))rHH))
Hence, we will need the leading right eigenvector for the matrix G.

MatrixG=(({{(mHfromH^2)*(1-(1/Hap)), mHfromH*mHfromD, mHfromD*mHfromH, (mHfromD^2)*(1-(1/(2*Dip)))}, {mHfromH*mDfromH*(1-(1/Hap)), mHfromH*mDfromD, mHfromD*mDfromH, mHfromD*mDfromD*(1-(1/(2*Dip)))}, {mDfromH*mHfromH*(1-(1/Hap)), mDfromH*mHfromD, mDfromD*mHfromH, mDfromD*mHfromD*(1-(1/(2*Dip)))}, {(mDfromH^2)*(1-(1/Hap)), mDfromH*mDfromD, mDfromD*mDfromH, (mDfromD^2)*(1-(1/(2*Dip)))}}));

We  first consider the right leading eigenvalue of matrix G.

MatrixG/.subconstraints/.Hap->ρH*Ntot/.Dip->(1-ρH)*Ntot/.Ntot->cNtot*(1/ϵ);
Limit[%,ϵ->0];
%//MatrixForm
Eigenvalues[%%]
Eigenvectors[%%%]

CharEq=Det[MatrixG-λ*IdentityMatrix[4]]/.subconstraints/.Hap->ρH*Ntot/.Dip->(1-ρH)*Ntot//Simplify;

CharEqlargeN=CharEq/.Ntot->cNtot*(1/ϵ)/.λ->λ0+λ1*ϵ+λ2*(ϵ^2)+λ3*(ϵ^3)+λ4*(ϵ^4)+λ5*(ϵ^5)+λ6*(ϵ^6)+λ7*(ϵ^7)+λ8*(ϵ^8)+λ9*(ϵ^9)+λ10*(ϵ^10)//Simplify;
Limit[%,ϵ->0]
Solve[%==0,λ0]

Series[CharEqlargeN,{ϵ,0,1}];
%/.λ0->1//Simplify

solλ1=Simplify[Solve[(Normal[%])==0,λ1],0<cNtot]

Order0λ=1;
Order1λ=λ1/.solλ1[[1]];

λGapprox=Order0λ+Order1λ*ϵ/.cNtot->Ntot*ϵ//Simplify

We  next derive the right leading eigenvector of matrix G.

RightVapprox=(({{rHHorder0+rHHorder1*ϵ+rHHorder2*(ϵ^2)}, {rHDorder0+rHDorder1*ϵ+rHDorder2*(ϵ^2)}, {rDHorder0+rDHorder1*ϵ+rDHorder2*(ϵ^2)}, {rDDorder0+rDDorder1*ϵ+rDDorder2*(ϵ^2)}}));
condRightV=(MatrixG.RightVapprox)-(λ*RightVapprox)/.λ->λGapprox+λ2*(ϵ^2)/.subconstraints/.Hap->ρH*Ntot/.Dip->(1-ρH)*Ntot/.Ntot->cNtot*(1/ϵ)//Simplify;

We consider order 0, we have

Series[condRightV,{ϵ,0,0}];
Normal[%]//Simplify

sol0=Solve[%==(({{0}, {0}, {0}, {0}}))&&rHHorder0+rHDorder0+rDHorder0+rDDorder0==const,{rHHorder0,rHDorder0,rDHorder0,rDDorder0}]

First order term

Series[condRightV,{ϵ,0,1}];
Normal[%]/.sol0[[1]]//Simplify

sol1=Solve[%==(({{0}, {0}, {0}, {0}}))&&rHHorder1+rHDorder1+rDHorder1+rDDorder1==0,{rHHorder1,rHDorder1,rDHorder1,rDDorder1}]

rightVapprox=(({{rHHorder0+rHHorder1*ϵ}, {rHDorder0+rHDorder1*ϵ}, {rDHorder0+rDHorder1*ϵ}, {rDDorder0+rDDorder1*ϵ}}))/.sol0[[1]]/.sol1[[1]]/.cNtot->Ntot*ϵ//Simplify

Total[rightVapprox]//Simplify
rightVapprox[[2]]==rightVapprox[[3]]

Using eigenvectors, we define Fst

CorrHHave=1-(rHH/(((Hap^2)/(((Hap+2*Dip))^2))*rHH+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*rHD+((((2*Dip))^2)/(((Hap+2*Dip))^2))*rDD))/.rHH->rightVapprox[[1,1]]/.rHD->rightVapprox[[2,1]]/.rDD->rightVapprox[[4,1]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.(OverTilde[Ntot])->Ntot*ϵ//FullSimplify
CorrHDave=1-(rHD/(((Hap^2)/(((Hap+2*Dip))^2))*rHH+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*rHD+((((2*Dip))^2)/(((Hap+2*Dip))^2))*rDD))/.rHH->rightVapprox[[1,1]]/.rHD->rightVapprox[[2,1]]/.rDD->rightVapprox[[4,1]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.(OverTilde[Ntot])->Ntot*ϵ//FullSimplify
CorrDDave=1-(rDD/(((Hap^2)/(((Hap+2*Dip))^2))*rHH+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*rHD+((((2*Dip))^2)/(((Hap+2*Dip))^2))*rDD))/.rHH->rightVapprox[[1,1]]/.rHD->rightVapprox[[2,1]]/.rDD->rightVapprox[[4,1]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.(OverTilde[Ntot])->Ntot*ϵ//FullSimplify

For large values of the population size, we do not obtain the same measure of Fst as Rousset’s method (by contrast to the Whitlock & Barton method using HD as a comparator, rather than the average, see below):

Factor[Normal[Series[{CorrHHave,CorrDDave}/.Ntot->(Ntot/ϵ),{ϵ,0,1}]]/.ϵ->1]

Factor[((%)/FNlarge)/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH]//FullSimplify

For example, in the symmetric case with large N and small (but not very small) movement rates, this measure yields Fst values that are half as large:

Normal[Series[%/.mHfromD->m/.mDfromH->m/.ρH->(2/3),{m,0,1}]]//Factor

Alternatively, we can focus on the case where sexuality is nearly complete (asexuality rare, with mHfromH and mDfromD of order ϵ):

Simplify[Normal[Series[CorrHHave/.mHfromD->1-mDfromD/.mDfromH->1-mHfromH/.mHfromH->((θHH)/Ntot)/.mDfromD->((θDD)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]],constraints]

which can be written as:

EqB1WB=((4*ρH(2*ρD) )/(-(((ρH-(2*ρD)))^2)+8*ρH*(2*ρD)*(ρH+(2*ρD))*(θHH+θDD)));
%%-%/.ρD->1-ρH//Factor

This differs substantially from Eq. (B.1) based on Rousset’s method (becoming negative as sexuality approaches one rather than approaching one as does Eq. (B.1)), likely because the coalescent events are poorly described by the leading eigenvalue only when the populations are nearly fully sexual:

EqB1=((ρH+(2*ρD))/(ρH+(2*ρD)+4* ρH*(2*ρD)*(θHH+θDD)));

Fst measured relative to HD by Whitlock and Barton (1997)

Using the eigenvectors obtained above, we define Fst

CorrHHhd=1-(rHH/rHD)/.rHH->rightVapprox[[1,1]]/.rHD->rightVapprox[[2,1]]/.rDD->rightVapprox[[4,1]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.cNtot->Ntot*ϵ//FullSimplify
CorrHDhd=1-(rHD/rHD)/.rHH->rightVapprox[[1,1]]/.rHD->rightVapprox[[2,1]]/.rDD->rightVapprox[[4,1]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.cNtot->Ntot*ϵ//FullSimplify
CorrDDhd=1-(rDD/rHD)/.rHH->rightVapprox[[1,1]]/.rHD->rightVapprox[[2,1]]/.rDD->rightVapprox[[4,1]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.cNtot->Ntot*ϵ//FullSimplify

For large values of the population size, we obtain the same measure of Fst as Rousset’s method:

Factor[Normal[Series[{CorrHHhd,CorrDDhd}/.Ntot->(Ntot/ϵ),{ϵ,0,1}]]/.ϵ->1]

Factor[%-FNlarge/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH]

Alternatively, we can focus on the case where sexuality is nearly complete (asexuality rare, with mHfromH and mDfromD of order ϵ):

Simplify[Normal[Series[CorrHHhd/.mHfromD->1-mDfromD/.mDfromH->1-mHfromH/.mHfromH->((θHH)/Ntot)/.mDfromD->((θDD)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]],constraints]

altEqB1WB=2*((ρH+(2*ρD))/(ρH+(2*ρD)+8*ρH*(2*ρD)*(θDD+θHH)));
%%-%/.ρD->1-ρH//Factor

EqB1WB

This differs substantially from Eq. (B.1) based on Rousset’s method (being nearly twice as large when asexuality is very rare and θDD+θHH near zero) , likely because the coalescent events are poorly described by the leading eigenvalue only when the populations are nearly fully sexual, as assumed by the Whitlock and Barton method.

The following graph illustrates the discrepancies that arise with the Whitlock & Barton method when asexuality is very rare, where the asymptotic properties focusing on the leading eigenvalue fail to capture the dynamics when there is a nearly strict alternation of generations.

Plot[Evaluate[{EqB1,EqB1WB,altEqB1WB}/.ρD->1-ρH/.(θHH+θDD)->sum/.ρH->(1/2)],{sum,0,2},PlotRange->{All,{-10,2}},PlotStyle->{{Black},{Blue}, {Red}}]

Complete symmetry case

MatrixG/.subconstraints/.mHfromD->m/.mDfromH->m/.Hap->ρH*Ntot/.Dip->(1-ρH)*Ntot/.ρH->(2/3)//Simplify;
%//MatrixForm
FullSimplify[Eigenvalues[%%],0<Ntot]
RightV=FullSimplify[Eigenvectors[%%%],0<Ntot]

RightVG=RightV[[4]]
RightVG[[2]]==RightVG[[3]]

CorrHHsym=1-(rHH/(((Hap^2)/(((Hap+2*Dip))^2))*rHH+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*rHD+((((2*Dip))^2)/(((Hap+2*Dip))^2))*rDD))/.rHH->RightVG[[1]]/.rHD->RightVG[[2]]/.rDD->RightVG[[4]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.ρH->(2/3)//FullSimplify
CorrHDsym=1-(rHD/(((Hap^2)/(((Hap+2*Dip))^2))*rHH+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*rHD+((((2*Dip))^2)/(((Hap+2*Dip))^2))*rDD))/.rHH->RightVG[[1]]/.rHD->RightVG[[2]]/.rDD->RightVG[[4]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.ρH->(2/3)//FullSimplify
CorrDDsym=1-(rDD/(((Hap^2)/(((Hap+2*Dip))^2))*rHH+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*rHD+((((2*Dip))^2)/(((Hap+2*Dip))^2))*rDD))/.rHH->RightVG[[1]]/.rHD->RightVG[[2]]/.rDD->RightVG[[4]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)/.ρH->(2/3)//FullSimplify

EqC4a=(1/3)*((Sqrt[9-4*m*(1-m)*(9-m*(1-m)*(((4*Ntot-3))^2))])-2*m*(1-m)*(4Ntot-3));
EqC4a-CorrHHsym//Simplify

CorrHHsym==CorrDDsym

Limit[CorrHHsym,m->0]
Limit[CorrHHsym,m->1]

Limit[CorrHHsym,m->(1/2)]
Simplify[%,1<Ntot]

Series[CorrHHsym/.m->cm*ϵ/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]//Simplify
Normal[%]/.cm->(m/ϵ)/.cNtot->Ntot*ϵ/.Ntot->3*Ndeem//Simplify

Series[CorrHHsym/.m->cm*ϵ/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]//Simplify
Normal[%]/.cm->(m/ϵ)/.cNtot->Ntot*ϵ/.Ntot->(3/2)*Nlocal//Simplify

This approximation is only half that obtained by Rousset’s method ((1/(1+4 θlocal)), where θlocal=m Nlocal) when θlocal is large:

Limit[((-4 θlocal+(Sqrt[1+16 ((θlocal)^2)]))/(1/(1+4 θlocal))),θlocal->Infinity]

Series[CorrHHsym/.m->cm*(ϵ^2)/.Ntot->cNtot*(1/ϵ),{ϵ,0,1}]//Simplify
Normal[%]/.cm->(m/(ϵ^2))/.cNtot->Ntot*ϵ/.Ntot->(3/2)*Nlocal//Simplify

Simplify[Series[CorrHHsym/.m->cm*ϵ/.Ntot->cNtot*(1/(ϵ^2)),{ϵ,0,1}],0<cm&&0<cNtot]
Simplify[Normal[%]/.cm->(m/ϵ)/.cNtot->Ntot*(ϵ^2)/.Ntot->(3/2)*Nlocal,{ϵ>0}]

Alternatively, we can focus on the case where sexuality is nearly complete (asexuality rare, with mHH and mDD of order ϵ):

FullSimplify[Normal[Series[CorrHHsym/.m->1-((θii)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]/.cNtot->Ntot*ϵ,{0<θii,Ntot>0}]

This approaches 1, as it should, when there is very little asexual reproduction (θii approaches 0):

Limit[%,θii->0]

But note that this is not generally true when there is very rare sexuality and the case is not fully symmetric.

Fst measured relative to HD for the complete symmetry case

If we instead compare sampling within a ploidy population to between ploidy populations (instead of the average) we get:

CorrHHsym2=1-(rHH/rHD)/.rHH->RightVG[[1]]/.rHD->RightVG[[2]]/.rDD->RightVG[[4]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)//FullSimplify
CorrDDsym2=1-(rDD/rHD)/.rHH->RightVG[[1]]/.rHD->RightVG[[2]]/.rDD->RightVG[[4]]/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)//FullSimplify

CorrHHsym2==CorrDDsym2

This is no longer EqC4a:

EqC4a=(1/3)*((Sqrt[9-4*m*(1-m)*(9-m*(1-m)*(((4*Ntot-3))^2))])-2*m*(1-m)*(4Ntot-3));
FullSimplify[EqC4a-CorrHHsym2,{m>0,Ntot>0}]

By multiplying top and bottom of the fraction in CorrHH2 by 3+6 (-1+m) m-(Sqrt[9+4 (-1+m) m (9+(-1+m) m (((3-4 Ntot))^2))]), we can simplify it to:

1+((8 (-1+m) m Ntot(3+6 (-1+m) m-(Sqrt[9+4 (-1+m) m (9+(-1+m) m (((3-4 Ntot))^2))])))/((((3+6 (-1+m) m))^2)-(9+4 (-1+m) m (9+(-1+m) m (((3-4 Ntot))^2)))));
FullSimplify[Factor[%],{m>0,Ntot>0}]

EqC4b=((3+2*m*(1-m)*(4*Ntot-9)-(Sqrt[9-4*m*(1-m)*(9-m*(1-m)*(((4*Ntot-3))^2))]))/(4*m*(1-m)*(2*Ntot-3)));
FullSimplify[EqC4b-CorrHHsym2,{m>0,Ntot>0}]

The alternative definition (dashed) is nearly twice the Fst value obtained by comparing to the average (solid), especially when the population size is large, but this alternative measure better matches Rousset’s Fst (black):

Eq3=((((1-2*m))^2)/((2/3)*Ntot-(((1-2*m))^2)*((2/3)*Ntot-1)));

GraphicsGrid[
{{Show[
Plot[Eq3/.Ntot->30,{m,0,1},PlotStyle->Black,AxesLabel->{"m","ϕ"}],
Plot[EqC4a/.Ntot->30,{m,0,1}],
Plot[EqC4b/.Ntot->30,{m,0,1},PlotStyle->Dashed]
],
Show[
Plot[Eq3/.Ntot->3000,{m,0,1},PlotStyle->Black,AxesLabel->{"m","ϕ"}],
Plot[EqC4a/.Ntot->3000,{m,0,1}],
Plot[EqC4b/.Ntot->3000,{m,0,1},PlotStyle->Dashed]
]}}]

which reduces to the classic Fst prediction when expressing the total population size in terms of the effective number of diploids:

Limit[CorrHHsym2,m->0]
Limit[CorrHHsym2,m->1]

Limit[CorrHHsym2,m->(1/2)]
Simplify[%,1<Ntot]

Series[CorrHHsym2/.m->cm*ϵ/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]//Simplify
Normal[%]/.cm->(m/ϵ)/.cNtot->Ntot*ϵ/.Ntot->3*Ndeem//Simplify

Series[CorrHHsym2/.m->cm*ϵ/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]//Simplify
Normal[%]/.cm->(m/ϵ)/.cNtot->Ntot*ϵ/.Ntot->(3/2)*Nlocal//Simplify

Multiplying the top and bottom of the first by 3-(Sqrt[9+64 (cm^2) (cNtot^2)]) and simplifying:

(1-((8 cm cNtot (3-(Sqrt[9+64 (cm^2) (cNtot^2)])))/(9-(9+64 (cm^2) (cNtot^2)))))//Simplify

Multiplying the top and bottom of the second by 1-(Sqrt[1+16 (m^2) (Nlocal^2)]) and simplifying:

1-((4 m Nlocal(1-(Sqrt[1+16 (m^2) (Nlocal^2)])))/(1-(1+16 (m^2) (Nlocal^2)))) //Simplify

For large θlocal=m Nlocal, this matches the Rousset result ((1/(1+4 θlocal))):

Limit[((1-((-1+(Sqrt[1+16 ((θlocal)^2)]))/(4 θlocal)))/(1/(1+4 θlocal))),θlocal->Infinity]

Very rate sexuality and large population for the complete symmetry case

Eq6=((((1-2*m))^2)/((2/3)*Ntot-(((1-2*m))^2)*((2/3)*Ntot-1)));
Series[%/.m->cm*(ϵ^2)/.Ntot->cNtot*(1/ϵ),{ϵ,0,1}]//Simplify
Normal[%]/.cm->(m/(ϵ^2))/.cNtot->Ntot*ϵ/.Ntot->3*Ndeem//Simplify

Series[CorrHHsym/.m->cm*(ϵ^2)/.Ntot->cNtot*(1/ϵ),{ϵ,0,1}]//Simplify
Normal[%]/.cm->(m/(ϵ^2))/.cNtot->Ntot*ϵ/.Ntot->3*Ndeem//Simplify

Series[CorrHHsym2/.m->cm*(ϵ^2)/.Ntot->cNtot*(1/ϵ),{ϵ,0,1}]//Simplify
Normal[%]/.cm->(m/(ϵ^2))/.cNtot->Ntot*ϵ/.Ntot->3*Ndeem//Simplify

Simplify[Series[CorrHHsym2/.m->cm*ϵ/.Ntot->cNtot*(1/(ϵ^2)),{ϵ,0,1}],0<cm&&0<cNtot]
Simplify[Normal[%]/.cm->(m/ϵ)/.cNtot->Ntot*(ϵ^2)/.Ntot->(3/2)*Nlocal,{ϵ>0}]

FullSimplify[Normal[Series[(1/3)*((Sqrt[9-4*m*(1-m)*(9-m*(1-m)*(((4*Ntot-3))^2))])-2*m*(1-m)*(4Ntot-3))/.m->(θ/Ntot)/.Ntot->Ntot/ϵ,{ϵ,0,0}]],{θ>0,Ntot>0}]

Defining θlocal = m Nlocal = (2/3) m Ntot, we can compare the three approximations, finding that the alternative definition using HD sampling (dashed) better matches Rousset (black) than does using the average (blue):

Show[
Plot[{-4 θlocal+(Sqrt[1+16 ((θlocal)^2)])},{θlocal,0,10},PlotStyle->Blue],
Plot[1-((-1+(Sqrt[1+16 ((θlocal)^2)]))/(4 θlocal)),{θlocal,0,10},PlotStyle->{Blue,Dashed}],
Plot[(1/(1+4 θlocal)),{θlocal,0,10},PlotStyle->Black]
]

unless migration is very rare (θlocal<<1), in which case it is the average and Rousset’s method that match (at ~1-4θlocal):

Show[
Plot[{-4 θlocal+(Sqrt[1+16 ((θlocal)^2)])},{θlocal,0,1},PlotRange->{Automatic,{0,1}},PlotStyle->{Blue}],
Plot[1-((4 θlocal)/(1+(Sqrt[1+16 ((θlocal)^2)]))),{θlocal,0,1},PlotStyle->{Blue,Dashed}],
Plot[(1/(1+4 θlocal)),{θlocal,0,1},PlotStyle->Black]
]

Alternatively, we can focus on the case where sexuality is nearly complete (asexuality rare, with mHH and mDD of order ϵ):

FullSimplify[Normal[Series[CorrHH2/.m->1-((θii)/Ntot)/.Ntot->(OverTilde[Ntot])*(1/ϵ),{ϵ,0,0}]]/.(OverTilde[Ntot])->ϵ Ntot,{0<θii,Ntot>0}]

This approaches 1, as it should, when there is very little asexual reproduction (θii approaches 0):

Limit[%,θii->0]

But note that this is not generally true when there is very rare sexuality and the case is not fully symmetric.

Appendix D

Variance effective size in Bessho and Otto (2022)

varianceNe=((ρH*ρD)/((cD^2)*ρH+2*(cH^2)*ρD))*Ntot;

For the global regulation model, we have

solH=(aH*wH+aD*wD-2*wD+(Sqrt[4*(1-aH)*(1-aD)*((f*wH)/2)*wD+(((aH*wH-aD*wD))^2)]))/(2*(aH*wH+((f*wH)/2)*(1-aH)-wD));
NeG=varianceNe/.cH->(((1-aH)*(f/2)*wH*((ρH)^2))/((1-aH)*(f/2)*wH*((ρH)^2)+(1-aD)*wD*((ρD)^2)))/.cD->(((1-aD)*wD*((ρD)^2))/((1-aH)*(f/2)*wH*((ρH)^2)+(1-aD)*wD*((ρD)^2)))/.ρD->1-ρH/.ρH->solH//Simplify

For the local regulation model, we have

NeL=varianceNe/.cH->((1+((aH*wH*ρH)/((1-aD)*wD*ρD)))/(2+((aH*wH*ρH)/((1-aD)*wD*ρD))+(2/f)*((aD*wD*ρD)/((1-aH)*wH*ρH))))/.cD->((1+(2/f)*((aD*wD*ρD)/((1-aH)*wH*ρH)))/(2+((aH*wH*ρH)/((1-aD)*wD*ρD))+(2/f)*((aD*wD*ρD)/((1-aH)*wH*ρH))))/.ρD->1-ρH;

Equations (from above)

Migration rates are backwards rates so that mHfromH+mHfromD sum to one:

constraints={mHfromH+mHfromD==1,mDfromH+mDfromD==1,ρD+ρH==1};

subconstraints={mHfromH->1-mHfromD,mDfromD->1-mDfromH};

FNlarge={(((1-X)*((2-X)*X*(2*ρD)-mHfromD(ρH+2*ρD)))/((2-X)*(X^2) *ρH*(2*ρD) *Ntot)),(( (1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((2-X)*(X^2) *ρH*(2*ρD) *Ntot))};

NextQHH=(((1-μ))^2)*((mHfromH^2)*((1/Hap)+(1-(1/Hap))*QHH)+mHfromH*mHfromD*QHD+mHfromD*mHfromH*QDH+(mHfromD^2)*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD));
NextQHD=(((1-μ))^2)*(mHfromH*mDfromH*((1/Hap)+(1-(1/Hap))*QHH)+mHfromH*mDfromD*QHD+mHfromD*mDfromH*QDH+mHfromD*mDfromD*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD));
NextQDH=(((1-μ))^2)*(mDfromH*mHfromH*((1/Hap)+(1-(1/Hap))*QHH)+mDfromH*mHfromD*QHD+mDfromD*mHfromH*QDH+mDfromD*mHfromD*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD));
NextQDD=(((1-μ))^2)*((mDfromD^2)*((1/(2*Dip))+(1-(1/(2*Dip)))*QDD)+mDfromH*mDfromD*QHD+mDfromD*mDfromH*QDH+(mDfromH^2)*((1/Hap)+(1-(1/Hap))*QHH));

The recursion equations for Q can be described as, (((1-μ))^2)A (Q+c);
(this equation corresponds to Eq. (9.30) in Rousset (2004))

MatrixA=(({{mHfromH*mHfromH, mHfromH*mHfromD, mHfromD*mHfromH, mHfromD*mHfromD}, {mHfromH*mDfromH, mHfromH*mDfromD, mHfromD*mDfromH, mHfromD*mDfromD}, {mDfromH*mHfromH, mDfromH*mHfromD, mDfromD*mHfromH, mDfromD*mHfromD}, {mDfromH*mDfromH, mDfromH*mDfromD, mDfromD*mDfromH, mDfromD*mDfromD}}));
VectC=(({{(1/Hap)*(1-QHH)}, {0}, {0}, {(1/(2*Dip))*(1-QDD)}}));
(({{NextQHH}, {NextQHD}, {NextQDH}, {NextQDD}}))-(((1-μ))^2)*MatrixA.((({{QHH}, {QHD}, {QDH}, {QDD}}))+VectC)//Simplify

Note that, we can write, A = TensorProduct[F, F]

MatrixF=(({{mHfromH, mHfromD}, {mDfromH, mDfromD}}));
TensorProduct[MatrixF,MatrixF]//MatrixForm
MatrixA//MatrixForm

Note that, the matrix A is a probaility matrix;

MatrixA.(({{1}, {1}, {1}, {1}}))/.mHfromD->1-mHfromH/.mDfromH->1-mDfromD//Simplify

Furthermore, this can be represented as, (((1-μ))^2)(GQ+Aδ) and (((1-μ))^2)(GQ+(I-G)1);
(this equation corresponds to Eq. (4.4) in Rousset (2004))

MatrixG=(({{(mHfromH^2)*(1-(1/Hap)), mHfromH*mHfromD, mHfromD*mHfromH, (mHfromD^2)*(1-(1/(2*Dip)))}, {mHfromH*mDfromH*(1-(1/Hap)), mHfromH*mDfromD, mHfromD*mDfromH, mHfromD*mDfromD*(1-(1/(2*Dip)))}, {mDfromH*mHfromH*(1-(1/Hap)), mDfromH*mHfromD, mDfromD*mHfromH, mDfromD*mHfromD*(1-(1/(2*Dip)))}, {(mDfromH^2)*(1-(1/Hap)), mDfromH*mDfromD, mDfromD*mDfromH, (mDfromD^2)*(1-(1/(2*Dip)))}}));
Vectδ=(({{(1/Hap)}, {0}, {0}, {(1/(2*Dip))}}));
(({{NextQHH}, {NextQHD}, {NextQDH}, {NextQDD}}))-(((1-μ))^2)*(MatrixG.(({{QHH}, {QHD}, {QDH}, {QDD}}))+MatrixA.Vectδ)//Simplify
(({{NextQHH}, {NextQHD}, {NextQDH}, {NextQDD}}))-(((1-μ))^2)*(MatrixG.(({{QHH}, {QHD}, {QDH}, {QDD}}))+((({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))-MatrixG).(({{1}, {1}, {1}, {1}})))/.subconstraints//Simplify

Asymptotic inbreeding effective size by Rousset (2004)

At equilibrium, we have
Q = (((1-μ))^2)(((I -(((1-μ))^2)G ))^(-1))(I-G)1.

We describe (I-G)1 by the linear combination by eigen vector of matrix G (v1, v2, v3, v4),
(I-G)1 = a1 v1 + a2 v2 +a3 v3 + a4 v4.

Using this, we have
Q=(((1-μ))^2)((1/(1-λ1))a1 v1+(1/(1-λ2))a2 v2+(1/(1-λ3))a3v3+(1/(1-λ4))a4v4),
where λ1, λ2, λ3, λ4 are eigen value of matrix G.

Furthermore, because of (((1-μ))^2)(λ^0)+(((1-μ))^4)(λ^1)+(((1-μ))^6)(λ^2)+...=((((1-μ))^2)/(1-(((1-μ))^2)λ)), we have
Q=((((1-μ))^2)(Subscript[λ, 1]^0)a1 v1+(((1-μ))^4)(Subscript[λ, 1]^1)a1 v1+...)+((((1-μ))^2)(Subscript[λ, 2]^0)a2 v2+(((1-μ))^4)(Subscript[λ, 2]^1)a2 v2+...)+((((1-μ))^2)(Subscript[λ, 3]^0)a3 v3+(((1-μ))^4)(Subscript[λ, 3]^1)a3 v3+...)+((((1-μ))^2)(Subscript[λ, 4]^0)a4 v4+(((1-μ))^4)(Subscript[λ, 4]^1)a4 v4+...).
It becomes
Q=(((1-μ))^2)((Subscript[λ, 1]^0)a1 v1+(Subscript[λ, 2]^0)a2 v2+(Subscript[λ, 3]^0)a3 v3+(Subscript[λ, 4]^0)a4 v4)+(((1-μ))^4)((Subscript[λ, 1]^1)a1 v1+(Subscript[λ, 2]^1)a2 v2+(Subscript[λ, 3]^1)a3 v3+(Subscript[λ, 4]^1)a4 v4)+....

Here the probability of IBD, Q, can be described as
Q=(((1-μ))^2)C(1)+(((1-μ))^4)C(2)+...,
where C indicates that the probability of coalecence.

Comparing these equations, we have the probabilty of coalecence using the eigenvalue of matrix G.

Using the leading eigenvalue of the matrix G, we have the inbreeding effective population size,
C(t)=((Subscript[λ, 1]^(t-1))a1 v1+(Subscript[λ, 2]^(t-1))a2 v2+(Subscript[λ, 3]^(t-1))a3 v3+(Subscript[λ, 4]^(t-1))a4 v4)

The inbreeding effective size comparing haploid WF population is defined as,
(1/Ne)=Limit[((C(t+1))/(1-(C(1)+C(2)+...+C(t)))),t-> infinity]

Using the probability of coalecence, we have
Limit[((C(t+1))/(1-(C(1)+C(2)+...+C(t)))),t-> infinity]=1-λ1,
where λ1 is the leading eigenvalue of matrix G.

Hence, we have the inbreeding effective population size,
Ne=(1/(1-λ1))

If we define the effective population size comparing with diploid WF model, it becomes
Ne=(1/2)(1/(1-λ1))

We consider the perturbation, G - A;

MatrixG-MatrixA//MatrixForm//Simplify

Because this perturbation is very small when we assume a large population size, we can derive the change in leading eigenvalue, using Δλ=((LeftEigenVector*(G-A)*RightEigenvector)/(LeftEigenVector*v)).
Because the matrix A is probability matrix, the leading eigenvalue and right eigenvector is one and unit vector.

LeftEigenVectorA=(({{VHH}, {VHD}, {VDH}, {VDD}}));
RightEigenVectorA=(({{1}, {1}, {1}, {1}}));
DeltaApproximationλG=((Transpose[LeftEigenVectorA].(MatrixG-MatrixA).RightEigenVectorA)/(Transpose[LeftEigenVectorA].RightEigenVectorA))//Simplify

The left eigenvector of matrix A (VHH, VHD, VDH, VDD) can be described as (VH*VH, VH*VD, VH*VD, VD*VD), where the VH and VD are the left eigenvector of the simpler 2x2 matrix F, normalized by their sum (VH=(mDfromH/(mHfromD+mDfromH)), VD=(mHfromD/(mHfromD+mDfromH))).

Eigenvalues[Simplify[Transpose[MatrixA]/.subconstraints]]
Eigenvectors[Simplify[Transpose[MatrixA]/.subconstraints]]
Eigenvalues[Simplify[Transpose[MatrixF]/.subconstraints]]
Eigenvectors[Simplify[Transpose[MatrixF]/.subconstraints]]

{((%[[1,1]])/(%[[1,1]]+%[[1,2]])),((%[[1,2]])/(%[[1,1]]+%[[1,2]]))}//Simplify

subV={VH->(mDfromH/(mHfromD+mDfromH)),VD->(mHfromD/(mHfromD+mDfromH))};

The leading eigenvalue of the matrix G is approximated as 1/Ne=λG=1-((VH^2)/( Hap))-(( (VD^2))/(2*Dip)).
This means that two gene lineages are in the same ploidy with probaility V and they can coalesce with probability 1/number of chromosomes (Eq.9.32, in Rousset 2004).

Part[DeltaApproximationλG[[1]],1]/.VHH->VH*VH/.VHD->VH*VD/.VDH->VD*VH/.VDD->VD*VD/.subconstraints/.subV//Simplify
%-(-((VH^2)/( Hap))-(( (VD^2))/(2*Dip)))/.subV//Simplify

Then, the leading eigenvalue of the matrix G is approximated as 1-((VH^2)/( Hap))-(( (VD^2))/(2*Dip)) and we have approximated inbreeding effective size comparing with classical diploid WF model, Ne=(1/2)(1/(1-λ1));

ApproximationλG=1-((VH^2)/( Hap))-(( (VD^2))/(2*Dip));
InbreedingNe=(1/2)*(1/(1-ApproximationλG))/.subV/.mHD->(((1-aH)*(f/2)*wH*ρH)/((1-aH)*(f/2)*wH*ρH+aD*wD*ρD))/.mDH->(((1-aD)*wD*ρD)/(aH*wH*ρH+(1-aD)*wD*ρD))/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.ρD->1-ρH//Simplify;

NeG-(InbreedingNe/.ρH->solH)//Simplify

NeL-InbreedingNe//Simplify

Effective size by Whitlock and Barton (1997)

Using the formulae Eq.(8) and (9) in Whitlock and Barton (1997), we have

Hap2=(Hap/2);
Nave=((Hap2+Dip)/2);
φH=((Hap2*mHfromH+Dip*mDfromH)/Hap2); 
(*Numerator represents number of individuals in the next generation that descend from haploids, with mHD being the fraction of the next generation of diploids that descend from haploids.*)
φD=((Hap2*mHfromD+Dip*mDfromD)/Dip);
Eq8inWB=(1/(4*(Nave^2)))*((((Hap2*((φH)^2))/2)*(1-CorrHHave)+((Dip*((φD)^2))/2)*(1-CorrDDave))+((Hap2^2)*((φH)^2)*CorrHHave+2*Hap2*Dip*φH*φD*CorrHDave+(Dip^2)*((φD)^2)*CorrDDave))//Simplify;

NeWB=(1/(2*Eq8inWB))/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)//Simplify

Series[(NeWB/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]//Simplify
NeWBlargeN=Normal[%*Ntot/.mDfromH->(((1-aH)*(f/2)*wH*ρH)/((1-aH)*(f/2)*wH*ρH+aD*wD*ρD))/.mHfromD->(((1-aD)*wD*ρD)/(aH*wH*ρH+(1-aD)*wD*ρD))/.ρD->1-ρH]//Simplify

NeG-(NeWBlargeN/.ρH->solH)//Simplify

NeL-NeWBlargeN//Simplify

Note that even if we measure Fst relative to a sample from haploids and diploids (HD) rather than the average sample across the full population, we must still rephrase the HD-based Fst values in terms of the average-based Fst values in order to use equations (8) and (9) of Barton and Whitlock:

((Corrij-fave)/(1-fave))/.Corrij->((Fij-FHD)/(1-FHD))/.fave->((Hap^2)/(((Hap+2*Dip))^2))*((FHH-FHD)/(1-FHD))+((Hap*(2*Dip))/(((Hap+2*Dip))^2))*((FHD-FHD)/(1-FHD))+((Hap*(2*Dip))/(((Hap+2*Dip))^2))*((FDH-FHD)/(1-FHD))+((((2*Dip))^2)/(((Hap+2*Dip))^2))*((FDD-FHD)/(1-FHD))//Simplify;
((Fij-Fave)/(1-Fave))/.Fave->((Hap^2)/(((Hap+2*Dip))^2))*FHH+((Hap*(2*Dip))/(((Hap+2*Dip))^2))*FHD+((Hap*(2*Dip))/(((Hap+2*Dip))^2))*FDH+((((2*Dip))^2)/(((Hap+2*Dip))^2))*FDD//Simplify;
%-%%//Simplify

altEq8inWB=(1/(4*(Nave^2)))*((((Hap2*((φH)^2))/2)*(1-(((CorrHHhd-fave)/(1-fave))))+((Dip*((φD)^2))/2)*(1-(((CorrDDhd-fave)/(1-fave)))))+((Hap2^2)*((φH)^2)*(((CorrHHhd-fave)/(1-fave)))+2*Hap2*Dip*φH*φD*(((CorrHDhd-fave)/(1-fave)))+(Dip^2)*((φD)^2)*(((CorrDDhd-fave)/(1-fave)))))/.fave->((Hap^2)/(((Hap+2*Dip))^2))*CorrHHhd+((2*Hap*(2*Dip))/(((Hap+2*Dip))^2))*CorrHDhd+((((2*Dip))^2)/(((Hap+2*Dip))^2))*CorrDDhd//Simplify;

Eq8inWB=(1/(4*(Nave^2)))*((((Hap2*((φH)^2))/2)*(1-CorrHHave)+((Dip*((φD)^2))/2)*(1-CorrDDave))+((Hap2^2)*((φH)^2)*CorrHHave+2*Hap2*Dip*φH*φD*CorrHDave+(Dip^2)*((φD)^2)*CorrDDave))//Simplify;

(recall that CorrHD2=0).

The above is the same equation as derived above, once put back in the terms used by Whitlock and Barton, and so gives the same Ne:

altEq8inWB-Eq8inWB/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*(1-ρH)//Factor

Figures

Figure 1-4: Fst with different m

Equations

subconstraints={mHfromH->1-mHfromD,mDfromD->1-mDfromH};

FaiHHstrict=(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*(2*ρD)-mHfromD*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2*ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH))/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH//Simplify;

FaiDDstrict=(((((1-X))^3)+Ntot*(1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((((1-X))^3)-Ntot*(1-X)*((((1-X))^2)*(ρH+2ρD)-mDfromD*ρH-mHfromH*(2*ρD))+(Ntot^2)*(X^2)*(2-X)*(2*ρD)*ρH))/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH//Simplify;

FaiAverage=((ρH)/(ρH+2*ρD))*FaiHHstrict+((2*ρD)/(ρH+2*ρD))*FaiDDstrict/.ρD->1-ρH//Simplify

FaiHHstrict/.mHfromD->1/.mDfromH->1//Simplify
FaiHHstrict/.mHfromD->1/.mDfromH->0//Simplify
FaiHHstrict/.mHfromD->0/.mDfromH->1//Simplify
FaiHHstrict/.mHfromD->0/.mDfromH->0//Simplify

FaiDDstrict/.mHfromD->1/.mDfromH->1//Simplify
FaiDDstrict/.mHfromD->1/.mDfromH->0//Simplify
FaiDDstrict/.mHfromD->0/.mDfromH->1//Simplify
FaiDDstrict/.mHfromD->0/.mDfromH->0//Simplify

Figure 1: Equal movement, Equal chromosome

FaiHHstrict/.mHfromD->m/.mDfromH->m//Simplify
FaiDDstrict/.mHfromD->m/.mDfromH->m//Simplify

FaiHHstrict/.mHfromD->m/.mDfromH->m/.ρH->(2/3)//Simplify
FaiDDstrict/.mHfromD->m/.mDfromH->m/.ρH->(2/3)//Simplify
%-%%//Simplify

Eq3=((((1-2*m))^2)/((2/3)*Ntot-(((1-2*m))^2)*((2/3)*Ntot-1)));
Eq3-FaiHHstrict/.mHfromD->m/.mDfromH->m/.ρH->(2/3)//Simplify

Limit[Eq3,m->0]
Limit[Eq3,m->(1/2)]
Limit[Eq3,m->1]

Plot[Eq3/.Ntot->10,{m,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Thickness[0.005],Blue},AspectRatio->0.75];
Plot[Eq3/.Ntot->100,{m,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Thickness[0.005],Black},AspectRatio->0.75];
Plot[Eq3/.Ntot->1000,{m,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Thickness[0.005],Red},AspectRatio->0.75];
Show[%,%%,%%%];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Figure 2: Equal movement

FaiHHstrict/.mHfromD->m/.mDfromH->m//Simplify
FaiDDstrict/.mHfromD->m/.mDfromH->m//Simplify

FaiHHstrict/.mHfromD->m/.mDfromH->m/.m->0//Simplify
FaiHHstrict/.mHfromD->m/.mDfromH->m/.m->1//Simplify
FaiDDstrict/.mHfromD->m/.mDfromH->m/.m->0//Simplify
FaiDDstrict/.mHfromD->m/.mDfromH->m/.m->1//Simplify

FaiHHstrict/.mHfromD->m/.mDfromH->m/.m->(1/2)//Simplify
FaiDDstrict/.mHfromD->m/.mDfromH->m/.m->(1/2)//Simplify

g1=Plot[FaiHHstrict/.mHfromD->m/.mDfromH->m/.ρH->0.9/.Ntot->100,{m,0,1},PlotRange->{{0,1},{-0.1,0.3}},PlotStyle->{Thickness[0.005],Red},AspectRatio->0.75];
g2=Plot[FaiHHstrict/.mHfromD->m/.mDfromH->m/.ρH->(2/3)/.Ntot->100,{m,0,1},PlotRange->{{0,1},{-0.1,0.3}},PlotStyle->{Thickness[0.005],Black,Dashed},AspectRatio->0.75];
g3=Plot[FaiHHstrict/.mHfromD->m/.mDfromH->m/.ρH->0.1/.Ntot->100,{m,0,1},PlotRange->{{0,1},{-0.1,0.3}},PlotStyle->{Thickness[0.005],Blue},AspectRatio->0.75];
Show[g1,g2,g3,PlotRange->{{0,1},{-0.05,0.2}}];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

g1=Plot[FaiDDstrict/.mHfromD->m/.mDfromH->m/.ρH->0.9/.Ntot->100,{m,0,1},PlotRange->{{0,1},{-0.15,0.5}},PlotStyle->{Thickness[0.005],Red},AspectRatio->0.75];
g2=Plot[FaiDDstrict/.mHfromD->m/.mDfromH->m/.ρH->(2/3)/.Ntot->100,{m,0,1},PlotRange->{{0,1},{-0.15,0.5}},PlotStyle->{Thickness[0.005],Black,Dashed},AspectRatio->0.75];
g3=Plot[FaiDDstrict/.mHfromD->m/.mDfromH->m/.ρH->0.1/.Ntot->100,{m,0,1},PlotRange->{{0,1},{-0.15,0.5}},PlotStyle->{Thickness[0.005],Blue},AspectRatio->0.75];
Show[g1,g2,g3,PlotRange->{{0,1},{-0.15,0.2}}];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Figure 3: Asymmetry movement, Equal chromosome number

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

FQHH=FaiHHstrict/.ρH->(2/3)/.Ntot->100;
FQDD=FaiDDstrict/.ρH->(2/3)/.Ntot->100;

OmegaQHH=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
OmegaQDD=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQHH[[repeat,1]]=x1;
OmegaQHH[[repeat,2]]=x2;
OmegaQHH[[repeat,3]]=FQHH/.mHfromD->x1/.mDfromH->x2;

OmegaQDD[[repeat,1]]=x1;
OmegaQDD[[repeat,2]]=x2;
OmegaQDD[[repeat,3]]=FQDD/.mHfromD->x1/.mDfromH->x2;
];
];

ListPlot3D[OmegaQHH]
ListPlot3D[OmegaQHH,PlotRange->All]

ListPlot3D[OmegaQDD]
ListPlot3D[OmegaQDD,PlotRange->All]

ListContourPlot[OmegaQHH,Contours->{{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

ListContourPlot[OmegaQDD,Contours->{{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->(2/3)/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.005,0.05}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->(2/3)/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.005,0.05}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiDDstrict/.ρH->(2/3)/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.005,0.05}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiDDstrict/.ρH->(2/3)/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.005,0.05}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Figure 4: Asymmetry movement

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

FQHH=FaiHHstrict/.ρH->0.1/.Ntot->100;
FQDD=FaiDDstrict/.ρH->0.1/.Ntot->100;

OmegaQHH=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
OmegaQDD=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQHH[[repeat,1]]=x1;
OmegaQHH[[repeat,2]]=x2;
OmegaQHH[[repeat,3]]=FQHH/.mHfromD->x1/.mDfromH->x2;

OmegaQDD[[repeat,1]]=x1;
OmegaQDD[[repeat,2]]=x2;
OmegaQDD[[repeat,3]]=FQDD/.mHfromD->x1/.mDfromH->x2;
];
];

ListContourPlot[OmegaQHH,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05},{0.5}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

ListContourPlot[OmegaQDD,Contours->{{-0.1},{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05},{0.5}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->0.1/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->0.1/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

FQHH=FaiHHstrict/.ρH->0.5/.Ntot->100;
FQDD=FaiDDstrict/.ρH->0.5/.Ntot->100;

OmegaQHH=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
OmegaQDD=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQHH[[repeat,1]]=x1;
OmegaQHH[[repeat,2]]=x2;
OmegaQHH[[repeat,3]]=FQHH/.mHfromD->x1/.mDfromH->x2;

OmegaQDD[[repeat,1]]=x1;
OmegaQDD[[repeat,2]]=x2;
OmegaQDD[[repeat,3]]=FQDD/.mHfromD->x1/.mDfromH->x2;
];
];

ListContourPlot[OmegaQHH,Contours->{{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

ListContourPlot[OmegaQDD,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->0.5/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->0.5/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

FQHH=FaiHHstrict/.ρH->0.9/.Ntot->100;
FQDD=FaiDDstrict/.ρH->0.9/.Ntot->100;

OmegaQHH=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
OmegaQDD=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];
repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQHH[[repeat,1]]=x1;
OmegaQHH[[repeat,2]]=x2;
OmegaQHH[[repeat,3]]=FQHH/.mHfromD->x1/.mDfromH->x2;

OmegaQDD[[repeat,1]]=x1;
OmegaQDD[[repeat,2]]=x2;
OmegaQDD[[repeat,3]]=FQDD/.mHfromD->x1/.mDfromH->x2;
];
];

ListContourPlot[OmegaQHH,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

ListContourPlot[OmegaQDD,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->0.9/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiHHstrict/.ρH->0.1/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Supplementary: Average Fst, Equal chromosome number

FaiAverage/.ρH->(2/3)//FullSimplify
Eq8=((((1-X))^2)/((2/3)*Ntot-(((1-X))^2)*((2/3)*Ntot-1)))/.X->mHfromD+mDfromH;
%-%%//Simplify

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

OmegaQave=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];

repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQave[[repeat,1]]=x1;
OmegaQave[[repeat,2]]=x2;
OmegaQave[[repeat,3]]=FaiAverage/.ρH->(2/3)/.Ntot->100/.mHfromD->x1/.mDfromH->x2;

];
];

ListContourPlot[OmegaQave,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05},{0.5}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->(2/3)/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->(2/3)/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Figure 5: Average Fst

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

OmegaQave=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];

repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQave[[repeat,1]]=x1;
OmegaQave[[repeat,2]]=x2;
OmegaQave[[repeat,3]]=FaiAverage/.ρH->0.1/.Ntot->100/.mHfromD->x1/.mDfromH->x2;

];
];

ListContourPlot[OmegaQave,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05},{0.5}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->0.1/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->0.1/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

OmegaQave=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];

repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQave[[repeat,1]]=x1;
OmegaQave[[repeat,2]]=x2;
OmegaQave[[repeat,3]]=FaiAverage/.ρH->0.5/.Ntot->100/.mHfromD->x1/.mDfromH->x2;

];
];

ListContourPlot[OmegaQave,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05},{0.5}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->0.5/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->0.5/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

MaxV=0.99;
MinV=0.01;
RepeatMax=101;

OmegaQave=Table[10,{i,(RepeatMax-1)*(RepeatMax-1)},{j,3}];

repeat=0;
For[repeat1=1,repeat1≤RepeatMax-1,repeat1=repeat1+1,
For[repeat2=1,repeat2≤RepeatMax-1,repeat2=repeat2+1,
repeat=repeat+1;
x1=MinV+(((MaxV-MinV))/RepeatMax)*(repeat1-1);
x2=MinV+(((MaxV-MinV))/RepeatMax)*(repeat2-1);

OmegaQave[[repeat,1]]=x1;
OmegaQave[[repeat,2]]=x2;
OmegaQave[[repeat,3]]=FaiAverage/.ρH->0.9/.Ntot->100/.mHfromD->x1/.mDfromH->x2;

];
];

ListContourPlot[OmegaQave,Contours->{{-0.01},{-0.001},{0,Darker[Green]},{0.001},{0.01},{0.05},{0.5}},PlotRange->All,ContourShading->None];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->0.9/.Ntot->100/.mDfromH->Table[i,{i,0.05,0.95,0.05}]],{mHfromD,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Plot[Evaluate[FaiAverage/.ρH->0.9/.Ntot->100/.mHfromD->Table[i,{i,0.05,0.95,0.05}]],{mDfromH,0,1},PlotStyle->Table[RGBColor[i,0,1-i],{i,0.05,0.95,0.05}],
PlotRange->{{0,1},{-0.05,0.1}}
];

Show[%,AxesOrigin->{0,0},AxesStyle->Directive[Black,AbsoluteThickness[0.5]],Frame->True,FrameStyle->AbsoluteThickness[0.5],FrameTicksStyle->Directive[AbsoluteThickness[0.5],FontFamily->"Times New Roman",FontSize->18],FrameLabel->{Style["",20,FontFamily->"Times New Roman"],Style["",20,FontFamily->"Times New Roman"]}]

Supplementary: Coalescence process considering FIS

Equations

NextQHH=(((1-μ))^2)*((mHfromH^2)*((1/Hap)+(1-(1/Hap))*QHH)+mHfromH*mHfromD*QHD+mHfromD*mHfromH*QDH+(mHfromD^2)*((1/(2*Dip))+(1/(2*Dip))*QDDw+(1-(2/(2*Dip)))*QDDb));
NextQHD=(((1-μ))^2)*(mHfromH*mDfromH*((1/Hap)+(1-(1/Hap))*QHH)+mHfromD*mDfromH*QDH+mHfromH*mDfromD*QHD+mHfromD*mDfromD*((1/(2*Dip))+(1/(2*Dip))*QDDw+(1-(2/(2*Dip)))*QDDb));
NextQDH=(((1-μ))^2)*(mDfromH*mHfromH*((1/Hap)+(1-(1/Hap))*QHH)+mDfromD*mHfromH*QDH+mDfromH*mHfromD*QHD+mDfromD*mHfromD*((1/(2*Dip))+(1/(2*Dip))*QDDw+(1-(2/(2*Dip)))*QDDb));
NextQDDb=(((1-μ))^2)*((mDfromD^2)*((1/(2*Dip))+(1/(2*Dip))*QDDw+(1-(2/(2*Dip)))*QDDb)+mDfromH*mDfromD*QHD+mDfromD*mDfromH*QDH+(mDfromH^2)*((1/Hap)+(1-(1/Hap))*QHH));
NextQDDw=(((1-μ))^2)*(mDfromD*QDDw+mDfromH*((1/Hap)+(1-(1/Hap))*QHH));

Recall that these are backwards rates so that mDD+mHD sum to one:

constraints={mHfromH+mHfromD==1,mDfromH+mDfromD==1,ρD+ρH==1};
subconstraints={mHfromH->1-mHfromD,mDfromD->1-mDfromH};

Exact FST and FIS

At equilibrium, we have

sol=Simplify[Solve[(NextQHH==QHH&&NextQDDb==QDDb&&NextQDDw==QDDw&&NextQHD==QHD&&NextQDH==QDH),{QHH,QHD,QDH,QDDb,QDDw}],constraints];

When there is not mutation, because of a finite population size, Q must be 1.

sol/.subconstraints/.μ->0//Simplify

EqQHH=Simplify[QHH/.sol[[1]],constraints];
EqQHD=Simplify[QHD/.sol[[1]],constraints];
EqQDH=Simplify[QDH/.sol[[1]],constraints];
EqQDDb=Simplify[QDDb/.sol[[1]],constraints];
EqQDDw=Simplify[QDDw/.sol[[1]],constraints];

Note that, we can show QHD = QDH,

EqQHD==EqQDH

Using Q, we define the F-like measures for the haploid-diploid population.

FaiHstrict=Simplify[Limit[((EqQHH-EqQHD)/(1-EqQHD)),μ->0],constraints]

FaiDbstrict=Simplify[Limit[((EqQDDb-EqQHD)/(1-EqQHD)),μ->0],constraints]

FaiDwstrict=Simplify[Limit[((EqQDDw-EqQHD)/(1-EqQHD)),μ->0],constraints]

FISstrict=Simplify[Limit[((EqQDDw-EqQDDb)/(1-EqQDDb))/.subconstraints,μ->0],constraints]

we can rewrite FaiHstrict as:

numerator=(((1-X))^3)+Hap*mHfromD*mDfromH*(1-X)+2*Dip*mDfromH*((((1-X))^2)-mHfromH)*(1-X)-(((1-X))^2)*(1+mDfromH*(1-2*X))//Simplify;
denominator=(((1-X))^3)+Hap*((((1-X))^2)-mDfromD)*(X+mDfromH*(1-2*X))+ 2*Dip*mDfromH*((((1-X))^2)-mHfromH)*(1-X)-2*Hap*Dip*mDfromH*(2-X)*(X^2)-(((1-X))^2)*(1+mDfromH*(1-2*X))//Simplify;
(numerator/denominator)//FullSimplify;
%-FaiHstrict/.X->mHfromD+mDfromH/.subconstraints/.ρD->1-ρH//Factor

While this doesn’t equal FaiHHstrict obtained without considering between and within diploid sampling, it does when the population size is large (here assuming intermediate migration rates):

FaiHstrict/.Hap->Ntot*ρH/.Dip->Ntot*ρD//Simplify

If N is large (total population size) relative to the migration rates (i.e., none of the migration rates is very rare) these reduce to the same FNlarge obtained above without considering sampling within and between diploids:

Factor[Normal[Series[{FaiHstrict,FaiDbstrict}/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.Ntot->(Ntot/ϵ),{ϵ,0,1}]]]/.ϵ->1

From above (ignoring between/within), we derived the following in this case:

FNlarge={(((1-X)*((2-X)*X*(2*ρD)-mHfromD(ρH+2*ρD)))/((2-X)*(X^2) *ρH*(2*ρD) *Ntot)),(( (1-X)*((2-X)*X*ρH-mDfromH*(ρH+2*ρD)))/((2-X)*(X^2) *ρH*(2*ρD) *Ntot))};

As expected this equals the same large N approximation as found without considering within and between sampling:

Factor[%%-FNlarge/.X->mDfromH+mHfromD/.ρD->1-ρH/.subconstraints]

Factor[Normal[Series[FaiDwstrict/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.Ntot->(Ntot/ϵ),{ϵ,0,1}]]]/.ϵ->1
Factor[%-((mHfromD*(-(2-X)*X*ρH+mDfromH*(ρH+2*ρD)))/(mDfromH*(2-X)*(X^2)*ρH*(2*ρD)*Ntot))/.X->mDfromH+mHfromD/.ρD->1-ρH/.subconstraints]

Factor[Normal[Series[FISstrict/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.Ntot->(Ntot/ϵ),{ϵ,0,1}]]]/.ϵ->1
Factor[%-((mDfromD*(-(2-X)*X*ρH+mDfromH*(ρH+2*ρD)))/(mDfromH*(2-X)*X*ρH*(2*ρD)*Ntot))/.X->mDfromH+mHfromD/.ρD->1-ρH/.subconstraints]

Rare sexuality in a large population

We consider the case when there is frequent asexual reproduction (a≒1, and then mHfromD and mDfromH is small), and large population size. We assume;
・mHfromD = O(ϵ)
・mDfromH = O(ϵ)
・Ntot = O(1/ϵ)
These assumptions indicates, θHD = mHfromD Ntot = O(1) and θDH = mDfromH Ntot = O(1).

FaiH◆RS=Simplify[Normal[Series[FaiHstrict/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

FaiDb◆RS=Simplify[Normal[Series[FaiDbstrict/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

FaiDw◆RS=Simplify[Normal[Series[FaiDwstrict/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

FIS◆RS=Simplify[Normal[Series[FISstrict/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]

FaiH◆RS

FaiDave◆RS=Simplify[Normal[Series[(1/(2*Dip-1))*FaiDwstrict+((2*Dip-2)/(2*Dip-1))*FaiDbstrict/.subconstraints/.Hap->Ntot*ρH/.Dip->Ntot*ρD/.mHfromD->((θHD)/Ntot)/.mDfromH->((θDH)/Ntot)/.Ntot->cNtot*(1/ϵ),{ϵ,0,0}]]]