About PAR: The distinct evolutionary dynamics of the Pseudoautosomal Region

Clark (1988) showed that, with tight enough linkage to the sex-determining region, a polymorphism would be maintained even with symmetrical overdominant selection with equal strength in both sexes.

Here we show that the same is true even if overdominance is not symmetrical and individuals have fitness:

Genotypes/fitness

basic commands

Off[General::"spell1"];
Off[General::"spell"];
SP[x_]:=Simplify[PowerExpand[x]]

genotypes

female genotypes:

(({{x1, T/T}, {x2, T/t}, {x3, t/t}}))

male genotypes:
(({{y1, T/T}, {y2c
y2t, T/t
t/T}, {y3, t/t}}))

*  y2c refers to T on X, t on Y
*  y2t refers to t on X, T on Y

fitness

fitness in females (natural selection):

(({{WTT, WTt, Wtt}, {1, 1-hfs, 1-s}}))

fitness in males (natural selection):

(({{MTT, MTt, Mtt}, {1-z, 1-hmz, 1}}))

subfit={WTT->1,WTt->1-h s, Wtt->1-s, MTT->1-z, MTt->1-hm z,Mtt->1};

Recursions (mating first, as in Clark)

mating  [ENTER]

*y2c refers to T on X,t on Y
*y2t refers to t on X,T on Y

females:

x1m=x1 y1+x1 (y2c*(1-r)+y2t*r)+1/2 x2 y1+1/2x2 (y2c*(1-r)+y2t*r);

x2m=x1 (y2c*r+y2t*(1-r))+x1 y3+1/2 x2 y1+1/2x2 (y2c*(1-r)+y2t*r)+1/2x2 (y2c*r+y2t*(1-r)) +1/2 x2 y3+x3 y1+x3 (y2c*(1-r)+y2t*r);

x3m=1/2 x2 (y2c*r+y2t*(1-r))+1/2 x2 y3+x3 (y2c*r+y2t*(1-r))+x3 y3;

males:

y1m=x1 y1+x1 (y2c*r+y2t*(1-r))+1/2 x2 y1+1/2x2 (y2c*r+y2t*(1-r));

y2cm=x1 (y2c*(1-r)+y2t*r)+x1 y3+1/2 x2(y2c*(1-r)+y2t*r) +1/2 x2 y3;

y2tm=1/2 x2 y1+1/2 x2 (y2c*r+y2t*(1-r)) +x3 y1+x3 (y2c*r+y2t*(1-r));

y3m=1/2 x2 (y2c*(1-r)+y2t*r)+1/2 x2 y3+x3 (y2c*(1-r)+y2t*r)+x3 y3;

Factor[x1m+x2m+x3m]

Factor[y1m+y2cm+y2tm+y3m]

natural selection  [ENTER]

females:

eqnx1=x1m*WTT/(x1m*WTT+x2m*WTt+x3m*Wtt);

eqnx2=x2m*WTt/(x1m*WTT+x2m*WTt+x3m*Wtt);

eqnx3=x3m*Wtt/(x1m*WTT+x2m*WTt+x3m*Wtt);

Factor[eqnx1+eqnx2+eqnx3]

males:

eqny1=y1m*MTT/(y1m*MTT+(y2cm+y2tm)*MTt+y3m*Mtt);

eqny2c=y2cm*MTt/(y1m*MTT+(y2cm+y2tm)*MTt+y3m*Mtt);

eqny2t=y2tm*MTt/(y1m*MTT+(y2cm+y2tm)*MTt+y3m*Mtt);

eqny3=y3m*Mtt/(y1m*MTT+(y2cm+y2tm)*MTt+y3m*Mtt);

Factor[eqny1+(eqny2c+eqny2t)+eqny3]

Conditions for a protected polymorphism

Here we solve for the conditions that are guaranteed to maintain a polymorphism because the fixation states (T only and t only) are both unstable to the introduction of the alternate allele.

Because we are focused on an autosomal locus and assume 50:50 sex ratios, the frequencies among sons is the same as the frequencies among daughters, so we need only solve for the equilibrium among daughters.

subEqual={x3->1-x1-x2,y3->1-y1-y2c-y2t};

Stability matrix when fixed for T:

mat={{D[eqnx1/.subEqual,x1],D[eqnx1/.subEqual,x2],D[eqnx1/.subEqual,y1],D[eqnx1/.subEqual,y2c],D[eqnx1/.subEqual,y2t]},
{D[eqnx2/.subEqual,x1],D[eqnx2/.subEqual,x2],D[eqnx2/.subEqual,y1],D[eqnx2/.subEqual,y2c],D[eqnx2/.subEqual,y2t]},
{D[eqny1/.subEqual,x1],D[eqny1/.subEqual,x2],D[eqny1/.subEqual,y1],D[eqny1/.subEqual,y2c],D[eqny1/.subEqual,y2t]},
{D[eqny2c/.subEqual,x1],D[eqny2c/.subEqual,x2],D[eqny2c/.subEqual,y1],D[eqny2c/.subEqual,y2c],D[eqny2c/.subEqual,y2t]},
{D[eqny2t/.subEqual,x1],D[eqny2t/.subEqual,x2],D[eqny2t/.subEqual,y1],D[eqny2t/.subEqual,y2c],D[eqny2t/.subEqual,y2t]}}/.x1->1/.x2->0/.y1->1/.y2c->0/.y2t->0//Factor

Definining the fitnesses relative to the fitness of the resident (e.g., wTt = WTt/WTT):

mat/.WTt->wTt*WTT/.MTt->mTt*MTT/.Wtt->wtt*WTT/.Mtt->mtt*MTT//Factor

Factor[Det[%-λ*IdentityMatrix[5]]]

This gives the condition (4) given in the Appendix:

char=(1/2) (mTt^2) (1-2 r) wTt+(1/2) (-2 mTt+2 mTt r-wTt) (λ^2)+(λ^3);

A similar analysis can be performed on the opposite equilibrium (with allele t fixed), but the results are symmetrical.

General equations for finding equilibria

The following equations describe a new coordinate system, measuring the frequency of A among the gametes (xf: in females, xm: on X from males, y: on Y from males):

subNew=Solve[{x1+x2/2==xf,y1+(y2c*(1-r)+y2t*r)==xm,y1+(y2c*r+y2t*(1-r))==y,x3==1-x1-x2,y3==1-y1-y2c-y2t},{x1,x2,x3,y1,y2c,y2t,y3}]//Flatten

{x1m,x2m,x3m,y1m,y2cm,y2tm,y3m}/.subNew//Factor

eqn1=Numerator[eqnx1+eqnx2/2-xf/.subNew//Factor]

eqn2=Numerator[eqny1+eqny2c*(1-r)+eqny2t*r-xm/.subNew//Factor]

eqn3=Numerator[eqny1+eqny2c*r+eqny2t*(1-r)-y/.subNew//Factor]

sub1=Simplify[Flatten[Solve[eqn3==0,xf]]]

sub2=Simplify[Flatten[Solve[(eqn1/.sub1)==0,xm]]]

There are potentially 7 solutions, including the fixation states y = 0 and y =1:

equileqn=Numerator[Factor[eqn2/.sub1/.sub2]]

At r = 0, there is only one possible internal solution with respect to y (recall that we could have different alleles fixed with X and with Y):

equileqn/.r->0/.subfit//Factor

Solve[%==0,y]//Simplify

At r = 1/2, we have the following possible internal solutions:

equileqn/.r->1/2/.subfit//Factor

quadratic=(1-hm z-2 y z+2 hm y z+(y^2) z-2 hm (y^2) z+hm (y^2) (z^2));

The quadratic is always positive at y = 0 and 1, but it's slope is always negative at y=1 (assuming 0<=hm<=1), so the quadratic does not give relevant roots:

quadratic/.y->0//Factor

quadratic/.y->1//Factor

D[quadratic,y]/.y->0//Factor

D[quadratic,y]/.y->1//Factor

The internal roots that exist with r = 1/2 must thus come from the solutions to the following cubic:

cubic=-(s-h s-s y+2 h s y-hm z+h hm s z-y z+2 hm y z-2 s y z+3 h s y z-2 h hm s y z+2 s (y^2) z-4 h s (y^2) z+(hm^2) (z^2)-(hm^2) s (z^2)+3 hm y (z^2)-6 (hm^2) y (z^2)-2 hm s y (z^2)-h hm s y (z^2)+5 (hm^2) s y (z^2)+2 (y^2) (z^2)-10 hm (y^2) (z^2)+12 (hm^2) (y^2) (z^2)-2 h s (y^2) (z^2)+6 hm s (y^2) (z^2)+2 h hm s (y^2) (z^2)-8 (hm^2) s (y^2) (z^2)-2 (y^3) (z^2)+8 hm (y^3) (z^2)-8 (hm^2) (y^3) (z^2)+2 h s (y^3) (z^2)-4 hm s (y^3) (z^2)+4 (hm^2) s (y^3) (z^2))/2;

For weak selection, there is only one internal equilibrium and it's valid if and only if the two corners are unstable:

Normal[Series[cubic/.s->s*eps/.z->z*eps,{eps,0,1}]]

{y,1-y}/.Solve[%==0,y]//Factor

Internal Equilibrium without Sex Differences in Allele Frequency:
	Equal selection in males and females 
	(homozygotes not necessarily equal in fitness)

Next, we ask whether there is an equilibrium which would have no sex differences in allele frequency.  When this is unstable but the polymorphism is protected, we infer that the system, if it reaches an equilibrium, must reach an equilibrium with sex differences in allele frequency.

With selection acting equally on both sexes (allows homozygous fitnesses to differ), there is an equilibrium with T equally frequent on the X and Y:

Factor[{eqn1,eqn2,eqn3}/.Wtt->Mtt/.WTt->1/.WTT->MTT/.MTt->1/.xm->y/.xf->y]

Solve[%==0,y]

Factor[{eqn1,eqn2,eqn3}/.Wtt->Mtt/.WTt->1/.WTT->MTT/.MTt->1/.xm->y/.xf->y/.y->((1-Mtt)/(2-Mtt-MTT))]

The asymmetric equilibrium is, however, still hard to identify (roots of a quartic):

equileqn/.Wtt->Mtt/.WTt->1/.WTT->MTT/.MTt->1//Factor

Internal Equilibrium without Sex Differences in Allele Frequency:
	No selection in females 
	(homozygotes not necessarily equal in fitness)

Next, we ask whether there is an equilibrium which would have no sex differences in allele frequency.  When this is unstable but the polymorphism is protected, we infer that the system, if it reaches an equilibrium, must reach an equilibrium with sex differences in allele frequency.

With selection acting only on males (allows homozygous fitnesses to differ), there is an equilibrium with T equally frequent on the X and Y:

Factor[{eqn1,eqn2,eqn3}/.Wtt->1/.WTt->1/.WTT->1/.MTt->1/.xm->y/.xf->y]

Solve[%==0,y]

Factor[{eqn1,eqn2,eqn3}/.Wtt->1/.WTt->1/.WTT->1/.MTt->1/.xm->y/.xf->y/.y->((1-Mtt)/(2-Mtt-MTT))]

Here are the asymmetrical equilibria with selection in males only:

equileqn/.Wtt->1/.WTt->1/.WTT->1/.MTt->1//Factor

The solution  (-r+2 y-Mtt y-MTT y-2 (y^2)+Mtt (y^2)+MTT (y^2))  = 0 leaves xf and xm undefined (note that the denominator becomes zero):

xf/.sub1/.Wtt->1/.WTt->1/.WTT->1/.MTt->1//Factor

xm/.sub2/.Wtt->1/.WTt->1/.WTT->1/.MTt->1//Factor

Leaving us with five equilibria: the two fixation states, the central one at y->((1-Mtt)/(2-Mtt-MTT))and the other two asymmetrical states:

Solve[(r-Mtt r-2 (r^2)-y+Mtt y+MTT y-Mtt MTT y+2 r y-2 MTT r y+(y^2)-Mtt (y^2)-MTT (y^2)+Mtt MTT (y^2)-2 r (y^2)+Mtt r (y^2)+MTT r (y^2))==0,y]//Simplify

With equally fit homozygotes in males:

roots=Simplify[%/.Mtt->1-s/.MTt->1/.MTT->1-s]

These roots are (not complex) if r < s/4, coinciding with the conditions under which the central equilibrium becomes unstable.

equil[1,1,1,9/10,1,1/2,1/100][[1]]

Stability Analysis:  When do we expect sex differences in allele frequency
	Equal selection in males and females 
	(homozygotes not necessarily equal in fitness)

Finally, we ask when the internal equilibrium without sex differences in allele frequency (i.e., with T equally frequent on the X and Y) is unstable, indicating that the system when it equilibrates cannot be at a state where the allele frequencies are the same on the X and Y (recall that with overdominance in both sexes, the fixation states are unstable).

Recursions for the new variable system:

eqnxf=eqnx1+eqnx2/2/.subNew//Factor

eqnxm=eqny1+eqny2c*(1-r)+eqny2t*r/.subNew//Factor

eqny=eqny1+eqny2c*r+eqny2t*(1-r)/.subNew//Factor

Local stability matrix evaluated at the equilibrium with  T equally frequent on the X and Y:

mat={{D[eqnxf/.subEqual,xf],D[eqnxf/.subEqual,xm],D[eqnxf/.subEqual,y]},{D[eqnxm/.subEqual,xf],D[eqnxm/.subEqual,xm],D[eqnxm/.subEqual,y]},{D[eqny/.subEqual,xf],D[eqny/.subEqual,xm],D[eqny/.subEqual,y]}}/.Wtt->Mtt/.WTt->1/.WTT->MTT/.MTt->1/.xm->y/.xf->y/.y->((-1+Mtt)/(-2+Mtt+MTT))//Factor

char=-Collect[Det[%-λ*IdentityMatrix[3]],λ,Factor]

Notice that at r = 0, the characteristic polynomial is negative at λ = 1 (with overdominance MTT < 1 and Mtt < MTT), implying that the above cubic must have a root greater than one.  Thus, the equilibrium with no sex differences in allele frequency is unstable to perturbations:

Factor[char/.λ->1/.r->0]

Conversely, it is positive at r = 1/2, suggesting that for r large enough, the equilibrium becomes stable:

Factor[char/.λ->1/.r->1/2]

Here we can ask what value of the recombination rate is required for the leading eigenvalue to pass through one, implying a transition from stability to instability of the equilibrium with no sex differences:

Factor[char/.λ->1]

rcrit=Factor[r/.Flatten[Solve[%==0,r]]]

This is the same as Clark's (13) if the two homozygotes are equal in fitness:

Factor[rcrit/.MTT->v1/v2/.Mtt->v1/v2]

But the exact value of rcrit does depend on the relative fitnesses of the homozygotes.

Stability Analysis:  When do we expect sex differences in allele frequency
	No selection in females 
	(homozygotes not necessarily equal in fitness)

Finally, in the absence of selection in females, we ask when the internal equilibrium without sex differences in allele frequency (i.e., with T equally frequent on the X and Y) is unstable, indicating that the system when it equilibrates cannot be at a state where the allele frequencies are the same on the X and Y (recall that with overdominance in males and no selection in females, the fixation states are unstable).

Recursions for the new variable system:

eqnxf=eqnx1+eqnx2/2/.subNew//Factor

eqnxm=eqny1+eqny2c*(1-r)+eqny2t*r/.subNew//Factor

eqny=eqny1+eqny2c*r+eqny2t*(1-r)/.subNew//Factor

Local stability matrix evaluated at the equilibrium with  T equally frequent on the X and Y:

mat={{D[eqnxf/.subEqual,xf],D[eqnxf/.subEqual,xm],D[eqnxf/.subEqual,y]},{D[eqnxm/.subEqual,xf],D[eqnxm/.subEqual,xm],D[eqnxm/.subEqual,y]},{D[eqny/.subEqual,xf],D[eqny/.subEqual,xm],D[eqny/.subEqual,y]}}/.Wtt->1/.WTt->1/.WTT->1/.MTt->1/.xm->y/.xf->y/.y->((1-Mtt)/(2-Mtt-MTT))//Factor

char=-Collect[Det[%-λ*IdentityMatrix[3]],λ,Factor]

Notice that at r = 0, the characteristic polynomial is negative at λ = 1 (with overdominance MTT < 1 and Mtt < MTT), implying that the above cubic must have a root greater than one.  Thus, the equilibrium with no sex differences in allele frequency is unstable to perturbations:

Factor[char/.λ->1/.r->0]

Conversely, it is positive at r = 1/2, suggesting that for r large enough, the equilibrium becomes stable:

Factor[char/.λ->1/.r->1/2]

Here we can ask what value of the recombination rate is required for the leading eigenvalue to pass through one, implying a transition from stability to instability of the equilibrium with no sex differences:

Factor[char/.λ->1]

rcrit=Factor[r/.Flatten[Solve[%==0,r]]]

This is consistent with Clark's (21):

%/.Mtt->v4/v5/.MTT->v4/v5//Factor

assuming no selection in females (v1 = v2):

v1(v5-v4)^2/(2(v2+v1)v5^2+(2v2-6v1)v4 v5)/.v1->v2//Factor

In equation (17), Clark implies that there is a second rcrit that determines when the fixation states are unstable.  This is true for some cases where there are selective differences betwen males and females, but this second rcrit doesn't always apply.

In particular, in this example, with overdominance in males and no selection in females, the fixation states are always unstable and Clark's equation (17) should not apply.  Here we show that the problem is that there can be more than one eigenvalue greater than one.

Returning to the stability of the system when t is fixed:

mat={{D[eqnxf/.subEqual,xf],D[eqnxf/.subEqual,xm],D[eqnxf/.subEqual,y]},{D[eqnxm/.subEqual,xf],D[eqnxm/.subEqual,xm],D[eqnxm/.subEqual,y]},{D[eqny/.subEqual,xf],D[eqny/.subEqual,xm],D[eqny/.subEqual,y]}}/.Wtt->1/.WTt->1/.WTT->1/.MTt->1/.xm->0/.xf->0/.y->0//Factor

char=-Collect[Det[%-λ*IdentityMatrix[3]],λ,Factor]

Factor[char/.λ->1]

rcrit=Factor[r/.Flatten[Solve[%==0,r]]]

This is consistent with Clark's (17):

%/.Mtt->v4/v5/.MTT->v4/v5//Factor

assuming no selection in females (v1 = v2):

(v2 v5^2-2 v1 v4 v5+(2v1-v2)v4^2)/(2v2 v5^2-2v1 v4 v5)/.v1->v2//Factor

But note that in this case (het adv in males only), this critical value only defines where a second eigenvalue crosses one.  There already is an eigenvalue above one:

char/.r->rcrit//Factor

Solve[%==0,λ]

Showing that the  ((1+(Sqrt[1+8 Mtt]))/(4 Mtt)) eigenvalue is already greater than one when the second eigenvalue passes through one: 

	 ((1+(Sqrt[1+8 Mtt]))/(4 Mtt))>1
if

	 1+(Sqrt[1+8 Mtt])>4 Mtt
	
if

	(Sqrt[1+8 Mtt])>4 Mtt - 1
	
if 4 Mtt - 1 is negative, this always holds true.  Even if it is positive, we can continue by squaring:

	1+8 Mtt >(4 Mtt - 1)^2
	
which equals:

	16 (1 - Mtt) Mtt > 0
	
which is always true with overdominance in males.  So there is always an eigenvalue greater than one, even when r passes through rcrit.

Plot[((1+(Sqrt[1+8 Mtt]))/(4 Mtt)),{Mtt,0,1}]

Numerical calculations and figures

Numerical calculations

Clear[equil]
equil[WTT_,WTt_,Wtt_,MTT_,MTt_,Mtt_,r_]:=equil[WTT,WTt,Wtt,MTT,MTt,Mtt,r]=Simplify[NSolve[{-2 Wtt xf+WTt xf+2 Wtt (xf^2)-2 WTt (xf^2)+WTt xm+2 Wtt xf xm-4 WTt xf xm+2 WTT xf xm-2 Wtt (xf^2) xm+4 WTt (xf^2) xm-2 WTT (xf^2) xm==0,-MTt xf+MTt r xf+Mtt xm-Mtt xf xm+MTt xf xm-MTt r y+MTt xf y-MTT xf y-Mtt xm y+MTt xm y+Mtt xf xm y-2 MTt xf xm y+MTT xf xm y==0,-MTt r xf+Mtt y-MTt y+MTt r y-Mtt xf y+2 MTt xf y-MTT xf y-Mtt (y^2)+MTt (y^2)+Mtt xf (y^2)-2 MTt xf (y^2)+MTT xf (y^2)==0},{xm,xf,y}],{xm>0,xf>0,y>0}]

{xm,xf,y}/.equil[1/10,1,1/10,1/10,1,1/10,1/100]

{xm,xf,y}/.equil[1/10,1,1/10,1/10,1,1/10,0.03]

{xm,xf,y}/.equil[1/10,1,1/10,1/10,1,1/10,0.04]

The case with lethal overdominance in both sexes has to be solved explicitly:

Simplify[Solve[{-2 Wtt xf+WTt xf+2 Wtt (xf^2)-2 WTt (xf^2)+WTt xm+2 Wtt xf xm-4 WTt xf xm+2 WTT xf xm-2 Wtt (xf^2) xm+4 WTt (xf^2) xm-2 WTT (xf^2) xm==0,-MTt xf+MTt r xf+Mtt xm-Mtt xf xm+MTt xf xm-MTt r y+MTt xf y-MTT xf y-Mtt xm y+MTt xm y+Mtt xf xm y-2 MTt xf xm y+MTT xf xm y==0,-MTt r xf+Mtt y-MTt y+MTt r y-Mtt xf y+2 MTt xf y-MTT xf y-Mtt (y^2)+MTt (y^2)+Mtt xf (y^2)-2 MTt xf (y^2)+MTT xf (y^2)==0}/.WTt->1/.MTt->1/.WTT->1-s/.Wtt->1-s/.MTT->1-s/.Mtt->1-s/.s->1,{xm,xf,y}]]

equil[0,1,0,0,1,0,r_]:=equil[0,1,0,0,1,0,r]=Simplify[{{xm->0,xf->0,y->0},{xm->(1/2),xf->(1/2),y->(1/2)},{xm->1,xf->1,y->1},{xm->((1+2 r+(Sqrt[1-4 (r^2)]))/(2+4 r)),xf->((1-2 r+(Sqrt[1-4 (r^2)]))/(2-4 r)),y->((1+2 r-(Sqrt[1-4 (r^2)]))/(2+4 r))},{xm->((1+2 r-(Sqrt[1-4 (r^2)]))/(2+4 r)),xf->((-1+2 r+(Sqrt[1-4 (r^2)]))/(-2+4 r)),y->((1+2 r+(Sqrt[1-4 (r^2)]))/(2+4 r))}},{xm>0,xf>0,y>0}]

For r = 0, these equilibria are fine and give complete sex specific alleles, but with higher r, the equilibria become unstable

Series[{((1+2 r+(Sqrt[1-4 (r^2)]))/(2+4 r)),((1-2 r+(Sqrt[1-4 (r^2)]))/(2-4 r)),((1+2 r-(Sqrt[1-4 (r^2)]))/(2+4 r))},{r,0,1}]

Series[{((1+2 r-(Sqrt[1-4 (r^2)]))/(2+4 r)),((-1+2 r+(Sqrt[1-4 (r^2)]))/(-2+4 r)),((1+2 r+(Sqrt[1-4 (r^2)]))/(2+4 r))},{r,0,1}]

Generating Figure 1 (X minus Y)

Figure 1a (selection equal in both males and females):

Clear[tab]
tab[s_]:=tab[s]=Table[{r,Block[{vec=Sort[Select[{xm,xf,y}/.equil[1-s,1,1-s,1-s,1,1-s,r],Im[#[[1]]]==0&&Im[#[[2]]]==0&&Im[#[[3]]]==0&&#[[1]]<(1-10^-10)&&#[[2]]<(1-10^-10)&&#[[3]]<(1-10^-10)&&#[[1]]>(10^-10)&&#[[2]]>(10^-10)&&#[[3]]>(10^-10)&],#1[[1]]>#2[[1]]&][[1]]},(1/3*vec[[1]]+2/3*vec[[2]]-vec[[3]])]},{r,0.0001,0.07,0.0005}]

Plotted is the average frequency on the X, (1/3 xm+2/3 xf), minus on the Y:
[Y-intercept rises from s = 1 at the bottom to s = 0 at the top]

ListPlot[{tab[1/10],tab[1/3],tab[1/2],tab[2/3],tab[9/10],{{0,0},{0,1}}},Joined->True,PlotRange->{{0,0.1},{0,1}},AspectRatio->1,PlotStyle->{RGBColor[0,0,1],RGBColor[0.2,0,0.8],RGBColor[0.4,0,0.6],RGBColor[0.6,0,0.4],RGBColor[0.8,0,0.2],RGBColor[1,0,0]}]

Figure 1b (selection only in males):

Clear[tab2]
tab2[s_]:=tab2[s]=Table[{r,Block[{vec=Sort[Select[{xm,xf,y}/.equil[1,1,1,1-s,1,1-s,r],Im[#[[1]]]==0&&Im[#[[2]]]==0&&Im[#[[3]]]==0&&#[[1]]<(1-10^-10)&&#[[2]]<(1-10^-10)&&#[[3]]<(1-10^-10)&&#[[1]]>(10^-10)&&#[[2]]>(10^-10)&&#[[3]]>(10^-10)&],#1[[1]]>#2[[1]]&][[1]]},(1/3*vec[[1]]+2/3*vec[[2]])-vec[[3]]]},{r,0.0001,0.3,0.001}]

Plotted is the average frequency on the X, (1/3 xm+2/3 xf), minus on the Y:
[X-intercept rises from s = 0 at the left to s = 1 at the right]

ListPlot[{tab2[1/10],tab2[1/3],tab2[1/2],tab2[2/3],tab2[9/10],tab2[1]},Joined->True,PlotRange->{{0,0.5},{0,1}},AspectRatio->1,PlotStyle->{RGBColor[0,0,1],RGBColor[0.2,0,0.8],RGBColor[0.4,0,0.6],RGBColor[0.6,0,0.4],RGBColor[0.8,0,0.2],RGBColor[1,0,0]}]

Generating Figure 2

Here we return to the characteristic polynomial derived in the "Conditions for a protected polymorphism" section.

char=(1/2) (mTt^2) (1-2 r) wTt+(1/2) (-2 mTt+2 mTt r-wTt) (λ^2)+(λ^3);

Assuming additivity for the figure:

Factor[char/.wTt->WTt/WTT/.mTt->MTt/MTT/.subfit/.hm->1/2/.h->1/2/.λ->1]

eqns=s/.Flatten[Solve[%==0,s]]//Simplify

This gives the conditions where fixation on T is unstable (to the right of these curves for different r values):

Show[{Plot[eqns/.r->1/100,{z,0,0.0199},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqns/.r->1/10,{z,0,0.18},PlotRange->{{0,1},{0,1}},PlotStyle->{Dashing[{0.01}]}],
Plot[eqns/.r->1/2,{z,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Orange,Dashing[{0.02}]}]}]

Examining the other fixation state:

matREV={{D[eqnx1/.subEqual,x1],D[eqnx1/.subEqual,x2],D[eqnx1/.subEqual,y1],D[eqnx1/.subEqual,y2c],D[eqnx1/.subEqual,y2t]},
{D[eqnx2/.subEqual,x1],D[eqnx2/.subEqual,x2],D[eqnx2/.subEqual,y1],D[eqnx2/.subEqual,y2c],D[eqnx2/.subEqual,y2t]},
{D[eqny1/.subEqual,x1],D[eqny1/.subEqual,x2],D[eqny1/.subEqual,y1],D[eqny1/.subEqual,y2c],D[eqny1/.subEqual,y2t]},
{D[eqny2c/.subEqual,x1],D[eqny2c/.subEqual,x2],D[eqny2c/.subEqual,y1],D[eqny2c/.subEqual,y2c],D[eqny2c/.subEqual,y2t]},
{D[eqny2t/.subEqual,x1],D[eqny2t/.subEqual,x2],D[eqny2t/.subEqual,y1],D[eqny2t/.subEqual,y2c],D[eqny2t/.subEqual,y2t]}}/.x1->0/.x2->0/.y1->0/.y2c->0/.y2t->0//Factor

matREV/.WTt->wTt*Wtt/.MTt->mTt*Mtt/.Wtt->wtt*Wtt/.Mtt->mtt*Mtt//Factor

Factor[Det[%-λ*IdentityMatrix[5]]]

This gives the same characteristic polynomial as before, but with different definitions for the relative fitnesses:

char=(1/2) (mTt^2) (1-2 r) wTt+(1/2) (-2 mTt+2 mTt r-wTt) (λ^2)+(λ^3);

Factor[char/.wTt->WTt/Wtt/.mTt->MTt/Mtt/.subfit/.hm->1/2/.h->1/2/.λ->1]

eqnsREF=s/.Flatten[Solve[%==0,s]]//Simplify

Show[{Plot[eqns/.r->1/100,{z,0,0.0199},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqns/.r->1/10,{z,0,0.18},PlotStyle->{Dashing[{0.01}]}],
Plot[eqns/.r->1/2,{z,0,1},PlotStyle->{Orange,Dashing[{0.02}]}],
Plot[eqnsREF/.r->1/100,{z,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqnsREF/.r->1/10,{z,0,1},PlotStyle->{Dashing[{0.01}]}],
Plot[eqnsREF/.r->1/2,{z,0,1},PlotStyle->{Orange,Dashing[{0.02}]}]},
AspectRatio->1]

Notice that with linkage, it is substantially easier to maintain polymorphism (black curves define a broader region).

With less fit heterozygotes, it is harder to maintain a polymorphism:

Factor[char/.wTt->WTt/WTT/.mTt->MTt/MTT/.subfit/.hm->2/3/.h->2/3/.λ->1]

eqns=s/.Flatten[Solve[%==0,s]]//Simplify

Factor[char/.wTt->WTt/Wtt/.mTt->MTt/Mtt/.subfit/.hm->2/3/.h->2/3/.λ->1]

eqnsREF=s/.Flatten[Solve[%==0,s]]//Simplify

Show[{Plot[eqns/.r->1/100,{z,0,0.0295},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqns/.r->1/10,{z,0,0.26},PlotStyle->{Dashing[{0.01}]}],
Plot[eqns/.r->1/2,{z,0,1},PlotStyle->{Orange,Dashing[{0.02}]}],
Plot[eqnsREF/.r->1/100,{z,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqnsREF/.r->1/10,{z,0,1},PlotStyle->{Dashing[{0.01}]}],
Plot[eqnsREF/.r->1/2,{z,0,1},PlotStyle->{Orange,Dashing[{0.02}]}]},
AspectRatio->1]

Consider r = 1/2 (orange curves).
In the cross-over (loop region, e.g., s = z = 0.2), both fixation states are stable.  It is only in the upper right with strong selection (e.g., s = z = 0.8) where there is a protected polymorphism.

Conversely, with more fit heterozygotes, it is easier to maintain a polymorphism:

Factor[char/.wTt->WTt/WTT/.mTt->MTt/MTT/.subfit/.hm->1/3/.h->1/3/.λ->1]

eqns=s/.Flatten[Solve[%==0,s]]//Simplify

Factor[char/.wTt->WTt/Wtt/.mTt->MTt/Mtt/.subfit/.hm->1/3/.h->1/3/.λ->1]

eqnsREF=s/.Flatten[Solve[%==0,s]]//Simplify

Show[{Plot[eqns/.r->1/100,{z,0,0.0295},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqns/.r->1/10,{z,0,0.15},PlotStyle->{Dashing[{0.01}]}],
Plot[eqns/.r->1/2,{z,0,1},PlotStyle->{Orange,Dashing[{0.02}]}],
Plot[eqnsREF/.r->1/100,{z,0,1},PlotRange->{{0,1},{0,1}},PlotStyle->{Black,Thickness[0.008]}],
Plot[eqnsREF/.r->1/10,{z,0,1},PlotStyle->{Dashing[{0.01}]}],
Plot[eqnsREF/.r->1/2,{z,0,1},PlotStyle->{Orange,Dashing[{0.02}]}]},
AspectRatio->1]