Off[General::"spell"]

One Patch (or mean field)

Equations

x1prime=x1*(x1*w1+x2*w2D+x3*w1+x4*w2H)-r*w2H*(x1*x4-x2*x3);
x2prime=x2*(x1*w2D+x2*w3+x3*w2H+x4*w3)+r*w2H*(x1*x4-x2*x3);
x3prime=x3*(x1*w1+x2*w2H+x3*w1+x4*w2R)+r*w2H*(x1*x4-x2*x3);
x4prime=x4*(x1*w2H+x2*w3+x3*w2R+x4*w3)-r*w2H*(x1*x4-x2*x3);

tot=Factor[x1prime+x2prime+x3prime+x4prime]

polyOLD=(w2D - w3)/(2*w2D - w1 - w3)

totstart=Factor[tot/.x3->0/.x4->0/.x2->1-x1/.x1->polyOLD]

sub={w1->(1-s)*a+(1-a),
w2D->(1-h11*s)*a+(1-k11*t)*(1-a),
w2H->(1-h12*s)*a+(1-k12*t)*(1-a),
w2R->(1-h22*s)*a+(1-k22*t)*(1-a),
w3->a+(1-t)*(1-a)};

Polymorphic equilibrium without modifier

x1prime/tot/.x3->0/.x4->0
x2prime/tot/.x3->0/.x4->0

Factor[x1-(x1prime/tot/.x3->0/.x4->0)/.x2->1-x1]

Solve[%==0,x1]

polyOLD=(w2D - w3)/(2*w2D - w1 - w3)

For this to be stable the following must be positive:

Simplify[{w2D-w1,w2D-w3}/.sub]

poly=Factor[polyOLD/.sub]

Simplify[{w2D-w1,w2D-w3}/.sub/.h11->1/2/.k11->1/2]

Simplify[{w1,w2D,w3}/.sub]

Simplify[{w2D-w1,w2D-w3}/.sub]

{a/.Solve[((w2D-w1)/.sub)==0,a],a/.
	Solve[((w2D-w3)/.sub)==0,a]}

Factor[D[poly,h11]]

Stability to invasion of modifier (general)

Getting the  matrix

a1=Factor[Normal[Series[x3prime/tot/.x3->e3*small/.x4->e4*small,{small,0,1}]]]/.small->1

a2=Factor[Normal[Series[x4prime/tot/.x3->e3*small/.x4->e4*small,{small,0,1}]]]/.small->1

ma=Factor[Coefficient[Expand[a1],e3]]

mb=Factor[Coefficient[Expand[a1],e4]]

mc=Factor[Coefficient[Expand[a2],e3]]

md=Factor[Coefficient[Expand[a2],e4]]

matrix1={{ma,mb},{mc,md}}

Matrix and eigenvalues

matrix1={{(w1*x1 + w2H*x2 - r*w2H*x2)/(w1*x1^2 + 2*w2D*x1*x2 + w3*x2^2),
(r*w2H*x1)/(w1*x1^2 + 2*w2D*x1*x2 + w3*x2^2)},
{(r*w2H*x2)/(w1*x1^2 + 2*w2D*x1*x2 + w3*x2^2),
(w2H*x1 - r*w2H*x1 + w3*x2)/(w1*x1^2 + 2*w2D*x1*x2 + w3*x2^2)}};

The difference between the eigenvalue and one is a measure of the spread rate of the modifier:

eigen1=Simplify[PowerExpand[Simplify[PowerExpand[(Eigenvalues[matrix1]-1)]]]]

matrix2=Simplify[matrix1/.x1->polyOLD/.x2->1-polyOLD];

eigen2=Simplify[PowerExpand[Simplify[PowerExpand[(Eigenvalues[matrix2]-1)]]]]

Factor[pa^2*w1+2*pa*(1-pa)*w2D+(1-pa)^2*w3/.pa->polyOLD]

Rewriting the eigenvalue in simpler terms:

notsqrtterm=(1/(2 ((w2D^2)-w1 w3)))(w1 w2D-2 (w2D^2)-w1 w2H+r w1 w2H+2 w2D w2H-2 r w2D w2H+w2D w3-w2H w3+r w2H w3)//Factor
sqrtterm=((((1/(2 ((w2D^2)-w1 w3)))))^2)((w1^2) (((w2D+(-1+r) w2H))^2)+(((w2D-w2H))^2) (w3^2)+2 r (w2D-w2H) w2H (w3^2)+(r^2) (w2H^2) (((-2 w2D+w3))^2)-2 w1 ((r^2) (w2H^2) (2 w2D-w3)+(((w2D-w2H))^2) w3+2 r (w2D-w2H) w2H w3))//Factor

Simplifying the notsqrtterm (replacing the mean fitness with just “wbar”):

notsqrtterm*totstart/wbar//Factor

Simplifying the sqrtpart:

Factor[sqrtpart/.r->0]

Factor[Coefficient[sqrtpart,r]]

Factor[Coefficient[sqrtpart,r^2]]

We can thus rewrite the leading eigenvalue (with the + sqrt term) as:

sm=((w2H-r w2H-w2D)/(2 wbar))+Sqrt[(((((w2D-w2H))^2) (((w1-w3))^2))/(4 (((w1-2 w2D+w3))^2) (wbar^2)))+r*(((w2D-w2H) w2H (((w1-w3))^2))/(2 (((w1-2 w2D+w3))^2) (wbar^2)))+r^2*((w2H^2)/(4 (wbar^2)))]

Or using the allele frequencies, pa->(w2D-w1)/(2 w2D-w3-w1) and pA->(w2D-w3)/(2 w2D-w3-w1):

sm2=((w2H-w2D-r w2H)/(2 wbar))+1/(2*wbar)*Sqrt[(((w2H-w2D-w2H*r))^2) (((pA-pa))^2)+4*pa*pA*r^2*(w2H^2)]

Weak modifier approximation

For a weak modifier (w2H ~ w2D), the eigenvalue is nearly:

Simplify[PowerExpand[Simplify[PowerExpand[Normal[Series[sm2,{w2H,w2D,1}]]/.pa->1-pA]]]]

as given by equation (3).

The next order term involves the recombination rate, showing that the eigenvalue declines as r rises (approximately as an inverse function of r):

Simplify[PowerExpand[Simplify[PowerExpand[Normal[Series[sm2,{w2H,w2D,2}]]-%/.pa->1-pA]]]]

Note that technically we have taken the Taylor series assuming that recombination isn’t too small relative to the difference between w2D and w2H (i.e., r is not so small that the next order term dominates).

Tight linkage approximation

Simplify[Simplify[((matrix2/.r->0)-IdentityMatrix[2])*totstart/Wbar]/.
	(w2D-w3)->pA*(2 w2D-w3-w1)/.
	(w1-w2D)->-pa*(2 w2D-w3-w1)]

The first describes the eigenvalue of the  rare mA haplotype and the second for the rare ma haplotype.  Invasion of the modifier occurs if these are greater than one, i.e., if the following are positive:

Simplify[(Eigenvalues[matrix2/.r->0]-1)*totstart/Wbar]

{%[[1]]*pA/polyOLD,%[[2]]*(1-pA)/(1-polyOLD)}//Factor

As given by equation (2). The modifier will invade in both cases if w2H>w2D.  The modifier allele linked to the rare allele will have a higher fitness when rare (because of its increased expression in heterozygotes).

Next exploring how the eigenvalues change with a small amount of recombination:

Simplify[D[((Eigenvalues[matrix2]-1)*totstart/Wbar)[[2]],r]];

Simplify[PowerExpand[Simplify[PowerExpand[%/.r->0]]]];
%(1-p)/(1-polyOLD)//Factor

Simplify[D[((Eigenvalues[matrix2]-1)*totstart/Wbar)[[1]],r]];

Simplify[PowerExpand[Simplify[PowerExpand[%/.r->0]]]];
%(p)/(polyOLD)//Factor

If linkage is tight, then increasing recombination will always reduce the selection for the modifier (reducing both eigenvalues).

Behaviour of selection with r (always goes down)

More generally, we can show that recombination causes the eigenvalue to decrease:

Simplify[D[sm2,r]*2*wbar/w2H]

The part on the right is always less than one in magnitude.  To prove this, 
first square the right hand side (if its less than one, it will remain so):

a1=Simplify[(((4 pa pA r w2H-(((pa-pA))^2) ((1-r) w2H-w2D))^2)/(4 pa pA (r^2) (w2H^2)+(((pa-pA))^2)((((1-r) w2H-w2D))^2)))]

This is equal to

a2=Simplify[((( r w2H-(((pa-pA))^2) (w2H-w2D))^2)/(r^2 w2H^2-2 (((pa-pA))^2) r w2H (w2H-w2D)+(((pa-pA))^2) (w2H-w2D)^2))]

Factor[a1/a2/.pa->1-pA]

This takes the form of (a-bc)^2 / (a^2-2abc+bc^2).  But since b=(((pa-pA))^2) is less 
than (or at most equal to) one, the denominator 

	a^2-2abc+bc^2  > a^2-2abc+b^2c^2 = (((a-b c))^2)

Hence the fraction is always less than (or at most equal to) one and so the derivative
with respect to r will always be negative (or at most zero).

lambda=1 analysis

charpoly=Factor[Det[matrix1-lambda*IdentityMatrix[2]]/.
	x1->polyOLD/.x2->1-polyOLD]

This is an upwards-facing quadratic:

Coefficient[charpoly,lambda^2]//Factor

The value at lambda=1 is (((1-pA) pA (w2H-w2D) (w2H(1-2 r)- w2D))/(Wbar^2)):

Factor[Factor[charpoly*totstart/Wbar*totstart/Wbar/.lambda->1]*pA(1-pA)/(polyOLD (1-polyOLD))]

as given in (A3a).

The slope of the characteristic polynomial at lambda=1 is:

Factor[D[charpoly*totstart/Wbar,lambda]/.lambda->1]

as given in (A3b).

For a modifier that decreases the fitness of the heterozygote, the slope at lambda = 1 is positive, so we must have that the characteristic polynomial is negative at lambda = 1 for the leading eigenvalue to be greater than one.  This requires that:
	(w2H (1-2 r) - w2D) > 0
which is never true for a modifier that decreases the fitness of the heterozygote.  Thus, such modifiers cannot invade.

For a modifier that increases the fitness of the heterozygote, the char poly will be negative at lambda=1 if:
	(w2H (1-2 r) - w2D) < 0
This is always true for r large enough, but is not true for r small.  But for r small, the slope of the charpoly is negative at lambda = 1, ensuring that both eigenvalues are greater than one.  Specifically, for r<(w2H-w2D)/(2*w2H) the char poly becomes positive at lambda =1, but from r=0 to r<(w2H-w2D)/(w2H) [a higher value] the slope is negative at lambda=1.  Therefore, the largest root is either imaginary or greater than one.

The Perron-Froebenius theorem guarantees the latter, demonstrating that the leading eigenvalue will be greater than one for modifiers that increase heterozygous fitness.

For graphing

General definitions

sub={w1->(1-s)*a+(1-a),
w2D->(1-h11*s)*a+(1-k11*t)*(1-a),
w2H->(1-h12*s)*a+(1-k12*t)*(1-a),
w2R->(1-h22*s)*a+(1-k22*t)*(1-a),
w3->a+(1-t)*(1-a)};

Clear[eigen2]
eigen2[h11_,h12_,s_,k11_,k12_,t_,a_,r_]=
	{(-(w1*w2D) + 2*w2D^2 + w1*w2H - r*w1*w2H - 2*w2D*w2H +

     2*r*w2D*w2H - w2D*w3 + w2H*w3 - r*w2H*w3 +

     (w1^2*w2D^2 - 2*w1^2*w2D*w2H + 2*r*w1^2*w2D*w2H +

        w1^2*w2H^2 - 2*r*w1^2*w2H^2 + r^2*w1^2*w2H^2 -

        4*r^2*w1*w2D*w2H^2 + 4*r^2*w2D^2*w2H^2 - 2*w1*w2D^2*w3 +

        4*w1*w2D*w2H*w3 - 4*r*w1*w2D*w2H*w3 - 2*w1*w2H^2*w3 +

        4*r*w1*w2H^2*w3 + 2*r^2*w1*w2H^2*w3 -

        4*r^2*w2D*w2H^2*w3 + w2D^2*w3^2 - 2*w2D*w2H*w3^2 +

        2*r*w2D*w2H*w3^2 + w2H^2*w3^2 - 2*r*w2H^2*w3^2 +

        r^2*w2H^2*w3^2)^(1/2))/(2*(-w2D^2 + w1*w3)),

  (-(w1*w2D) + 2*w2D^2 + w1*w2H - r*w1*w2H - 2*w2D*w2H +

     2*r*w2D*w2H - w2D*w3 + w2H*w3 - r*w2H*w3 -

     (w1^2*w2D^2 - 2*w1^2*w2D*w2H + 2*r*w1^2*w2D*w2H +

        w1^2*w2H^2 - 2*r*w1^2*w2H^2 + r^2*w1^2*w2H^2 -

        4*r^2*w1*w2D*w2H^2 + 4*r^2*w2D^2*w2H^2 - 2*w1*w2D^2*w3 +

        4*w1*w2D*w2H*w3 - 4*r*w1*w2D*w2H*w3 - 2*w1*w2H^2*w3 +

        4*r*w1*w2H^2*w3 + 2*r^2*w1*w2H^2*w3 -

        4*r^2*w2D*w2H^2*w3 + w2D^2*w3^2 - 2*w2D*w2H*w3^2 +

        2*r*w2D*w2H*w3^2 + w2H^2*w3^2 - 2*r*w2H^2*w3^2 +

        r^2*w2H^2*w3^2)^(1/2))/(2*(-w2D^2 + w1*w3))}/.sub//Evaluate;

Clear[eigen3]
eigen3[h11_,h12_,s_,k11_,k12_,t_,a_,r_]:=
	eigen3[h11,h12,s,k11,k12,t,a,r]=
	If[
	a*s - a*h11*s - k11*t + a*k11*t>0 &&
	-(a*h11*s) + t - a*t - k11*t + a*k11*t>0,
	Max[eigen2[h11,h12,s,k11,k12,t,a,r][[1]],
		eigen2[h11,h12,s,k11,k12,t,a,r][[2]]],
	0]

Clear[test]
test[h11_,h12_,s_,k11_,k12_,t_,a_,r_]:={w1,w2D,w3}/.sub//Evaluate

Clear[sol]
sol[h11_,h12_,s_,k11_,k12_,t_,a_,r_]=
	(2*(w2H-w2D)*(w1 - w2D)*(-w2D + w3))/(Wbar*(w1 - 2*w2D + w3)^2)/.
	Wbar->(-w2D^2 + w1*w3)/(w1 - 2*w2D + w3)/.sub;

Clear[solr0]
solr0[h11_,h12_,s_,k11_,k12_,t_,a_,r_]=
	Max[((w1 - w2D)*(-w2D + w2H))/(Wbar*(w1 - 2*w2D + w3)),
	((-w2D + w2H)*(-w2D + w3))/(Wbar*(w1 - 2*w2D + w3))]/.
	Wbar->(-w2D^2 + w1*w3)/(w1 - 2*w2D + w3)/.sub;

Plots with r

h11=0.5, h12=0.1, s=0.5
k11=0.1, k12=0.1, t=0.1
alpha =0.1, r=0.5

Plot[eigen3[0.5,0.1,0.5,0.1,0.1,0.1,0.1,r],
	{r,0,0.5},PlotRange->{{0,0.5},{0,0.02}}]

h11=0.5, h12=0.1, s=0.5
k11=0.1, k12=0.1, t=0.4
alpha =0.5

test[0.5,0.1,0.5,0.1,0.1,0.4,0.5,r]

h11=0.5, h12=0.4, s=0.5
k11=0.1, k12=0.1, t=0.1
alpha =0.1

test[0.5,0.4,0.5,0.1,0.1,0.1,0.1,r]

Show[
ListPlot[Table[{r,
	eigen3[0.5,0.1,0.5,0.1,0.1,0.4,0.5,r]},
	{r,0,0.5,0.01}],Joined->True,PlotRange->All],
Plot[sol[0.5,0.1,0.5,0.1,0.1,0.4,0.5,r],
	{r,0,0.5},PlotStyle->Dashing[{0.02,0.04}]],
ListPlot[Table[{r,
	solr0[0.5,0.1,0.5,0.1,0.1,0.4,0.5,r]},
	{r,0,0,0.01}],PlotStyle->PointSize[0.03]],
PlotRange->{{-0.01,0.5},{0,0.1}},AxesOrigin->{0,0}]

h11=0.5, h12=0.6, s=0.5
k11=0.1, k12=0.1, t=0.1
alpha =0.1

Show[
ListPlot[Table[{r,
	eigen3[0.5,0.4,0.5,0.1,0.1,0.1,0.1,r]},
	{r,0,0.024,0.001}],Joined->True,PlotRange->All],
ListPlot[Table[{r,
	eigen3[0.5,0.4,0.5,0.1,0.1,0.1,0.1,r]},
	{r,0.024,0.5,0.003}],Joined->True,PlotRange->All],
Plot[sol[0.5,0.4,0.5,0.1,0.1,0.1,0.1,r],
	{r,0,0.5},PlotStyle->Dashing[{0.02,0.04}]],
ListPlot[Table[{r,
	solr0[0.5,0.4,0.5,0.1,0.1,0.1,0.1,r]},
	{r,0,0,0.01}],PlotStyle->PointSize[0.03]],
PlotRange->{{-0.01,0.5},{0,0.005}},AxesOrigin->{0,0}]

Simulations

Getting the equations

x1prime=x1*(x1*w1+x2*w2D+x3*w1+x4*w2H)-r*w2H*(x1*x4-x2*x3);
x2prime=x2*(x1*w2D+x2*w3+x3*w2H+x4*w3)+r*w2H*(x1*x4-x2*x3);
x3prime=x3*(x1*w1+x2*w2H+x3*w1+x4*w2R)+r*w2H*(x1*x4-x2*x3);
x4prime=x4*(x1*w2H+x2*w3+x3*w2R+x4*w3)-r*w2H*(x1*x4-x2*x3);

tot=Factor[x1prime+x2prime+x3prime+x4prime]

polyOLD=(w2D - w3)/(2*w2D - w1 - w3)

totstart=Factor[tot/.x3->0/.x4->0/.x2->1-x1/.x1->polyOLD]

Simplify[x4prime/tot/.sub/.sub2]

The equations

sub={w1->(1-s)*a+(1-a),
w2D->(1-h11*s)*a+(1-k11*t)*(1-a),
w2H->(1-h12*s)*a+(1-k12*t)*(1-a),
w2R->(1-h22*s)*a+(1-k22*t)*(1-a),
w3->a+(1-t)*(1-a)};

sub2={x1->x1[g-1,a,s,h11,h12,h22,t,k11,k12,k22],
		x2->x2[g-1,a,s,h11,h12,h22,t,k11,k12,k22],
		x3->x3[g-1,a,s,h11,h12,h22,t,k11,k12,k22],
		x4->x4[g-1,a,s,h11,h12,h22,t,k11,k12,k22]};

Clear[poly]
poly[a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]=(w2D - w3)/(2*w2D - w1 - w3)/.sub;

Clear[tot]
tot[a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]=((-(w2D^2)+w1 w3)/(w1-2 w2D+w3))/.sub;

hetchange[a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]=w2H-w2D/.sub;

hetchange2[a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]=w2R-w2H/.sub;

Clear[x1,x2,x3,x4]
	x1[g_,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x1[g,a,s,h11,h12,h22,t,k11,k12,k22,r]=-(x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+a s) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a h11 s+k11 t-a k11 t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a s) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a h12 s+k12 t-a k12 t) x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+a h12 s+k12 t-a k12 t) (x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]))/((-1+a s) ((x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h11 s+k11 t-a k11 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a s) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h12 s+k12 t-a k12 t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(-1+a s) ((x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h12 s+k12 t-a k12 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-2 (1+(-1+a) t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h22 s+k22 t-a k22 t) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2));
x2[g_,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x2[g,a,s,h11,h12,h22,t,k11,k12,k22,r]=
		-(x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+a h11 s+k11 t-a k11 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(1+(-1+a) t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a h12 s+k12 t-a k12 t) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(1+(-1+a) t) x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+a h12 s+k12 t-a k12 t) (x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]))/((-1+a s) ((x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h11 s+k11 t-a k11 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a s) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h12 s+k12 t-a k12 t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(-1+a s) ((x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h12 s+k12 t-a k12 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-2 (1+(-1+a) t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h22 s+k22 t-a k22 t) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2));
x3[g_,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x3[g,a,s,h11,h12,h22,t,k11,k12,k22,r]=-(x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+a s) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a h12 s+k12 t-a k12 t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a s) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a h22 s+k22 t-a k22 t) x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+a h12 s+k12 t-a k12 t) (x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]))/((-1+a s) ((x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h11 s+k11 t-a k11 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a s) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h12 s+k12 t-a k12 t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(-1+a s) ((x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h12 s+k12 t-a k12 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-2 (1+(-1+a) t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h22 s+k22 t-a k22 t) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2));
		x4[g_,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x4[g,a,s,h11,h12,h22,t,k11,k12,k22,r]=-(x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+a h12 s+k12 t-a k12 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(1+(-1+a) t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+a h22 s+k22 t-a k22 t) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(1+(-1+a) t) x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+a h12 s+k12 t-a k12 t) (x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]))/((-1+a s) ((x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h11 s+k11 t-a k11 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a s) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h12 s+k12 t-a k12 t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+(-1+a s) ((x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2)+2 (-1+a h12 s+k12 t-a k12 t) x1[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-2 (1+(-1+a) t) x2[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]+2 (-1+a h22 s+k22 t-a k22 t) x3[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r]-(1+(-1+a) t) ((x4[-1+g,a,s,h11,h12,h22,t,k11,k12,k22,r])^2));

One-locus start

x1[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.1;
x2[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.9;
x3[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
x4[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;

ListPlot[Table[x1[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,250}],Joined->True]

polyOLD/.w1->0.95/.w3->0.91/.w2D->0.966

Factor[polyOLD/.sub]

Factor[1-polyOLD/.sub]

Let’s start at the equilibrium and confirm it remains there.

Reenter the previous section “The equations” and then Enter the following:

x1[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=((a h11 s-t+a t+k11 t-a k11 t)/(-a s+2 a h11 s-t+a t+2 k11 t-2 a k11 t));
x2[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=((-a s+a h11 s+k11 t-a k11 t)/(-a s+2 a h11 s-t+a t+2 k11 t-2 a k11 t));
x3[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
x4[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;

ListPlot[Table[x1[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,20}],Joined->True,PlotRange->{Automatic,{0,1}}]

Two-locus start

Now introducing the modifier at a frequency of 1%.

Reenter the section “The equations” and then enter the following:

x1[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=((a h11 s-t+a t+k11 t-a k11 t)/(-a s+2 a h11 s-t+a t+2 k11 t-2 a k11 t))*0.99;
x2[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=((-a s+a h11 s+k11 t-a k11 t)/(-a s+2 a h11 s-t+a t+2 k11 t-2 a k11 t))*0.99;
x3[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=((a h11 s-t+a t+k11 t-a k11 t)/(-a s+2 a h11 s-t+a t+2 k11 t-2 a k11 t))*0.01;
x4[0,a_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=((-a s+a h11 s+k11 t-a k11 t)/(-a s+2 a h11 s-t+a t+2 k11 t-2 a k11 t))*0.01;

poly[0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]

Generating Figure 3:

Show[
	ListPlot[Table[{gg,x3[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.01,0.01}]],
	ListPlot[Table[{gg,x3[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.001]+x4[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.001]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.02,0.01}]],
	ListPlot[Table[{gg,x3[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+x4[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.04,0.01}]],
	ListPlot[Table[{gg,x3[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+x4[gg,0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->2*poly[0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]*(1-poly[0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0])*(hetchange[0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0])/tot[0.1,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]/.p0->0.01,{t,0,6000},PlotStyle->Thin]
]

One Patch - with a cost

Equations

Let F be the fitness of Mm individuals and F^2 the fitness of mm individuals due to direct selection.
We will be testing to make sure that if the fitness of the modifier is 1/lambda, then the 
eigenvalue will equal one.

x1prime=x1*(x1*w1+x2*w2D+x3*w1*F+x4*w2H*F)-r*w2H*F*(x1*x4-x2*x3);
x2prime=x2*(x1*w2D+x2*w3+x3*w2H*F+x4*w3*F)+r*w2H*F*(x1*x4-x2*x3);
x3prime=x3*F*(x1*w1+x2*w2H+x3*w1*F+x4*w2R*F)+r*w2H*F*(x1*x4-x2*x3);
x4prime=x4*F*(x1*w2H+x2*w3+x3*F*w2R+x4*F*w3)-r*w2H*F*(x1*x4-x2*x3);

polyOLD=(w2D - w3)/(2*w2D - w1 - w3)

tot=Factor[x1prime+x2prime+x3prime+x4prime]

totstart=Factor[tot/.x3->0/.x4->0/.x2->1-x1/.x1->polyOLD]

Stability to invasion of modifier

General

a1=Normal[Normal[Series[x3prime/tot/.x3->e3*small/.x4->e4*small,{small,0,1}]]]/.small->1

a2=Normal[Normal[Series[x4prime/tot/.x3->e3*small/.x4->e4*small,{small,0,1}]]]/.small->1

ma=Factor[Coefficient[Expand[a1],e3]]

mb=Factor[Coefficient[Expand[a1],e4]]

mc=Factor[Coefficient[Expand[a2],e3]]

md=Factor[Coefficient[Expand[a2],e4]]

matrix1={{ma,mb},{mc,md}}

matrix2=Simplify[matrix1/.x1->polyOLD/.x2->1-polyOLD];

Recall from above that the eigenvalue without direct selection is 1+sm:

sm

If we set direct selection, F, to 1/(1+sm), we do get that the leading eigenvalue is 1:

FullSimplify[PowerExpand[Simplify[PowerExpand[Simplify[PowerExpand[(Eigenvalues[matrix2/.F->1/(1+sm)/.wbar->totstart])]]]]]]

Note that complicated term in the above is less than one (the term in front of the squareroot is negative on the numerator but positive in the denominator).

Two Patches - NEW

Equations

x1p=x1*(x1*w1+x2*w2D+x3*w1+x4*w2H)-r*w2H*(x1*x4-x2*x3);
x2p=x2*(x1*w2D+x2*w3+x3*w2H+x4*w3)+r*w2H*(x1*x4-x2*x3);
x3p=x3*(x1*w1+x2*w2H+x3*w1+x4*w2R)+r*w2H*(x1*x4-x2*x3);
x4p=x4*(x1*w2H+x2*w3+x3*w2R+x4*w3)-r*w2H*(x1*x4-x2*x3);

y1p=y1*(y1*v1+y2*v2D+y3*v1+y4*v2H)-r*v2H*(y1*y4-y2*y3);
y2p=y2*(y1*v2D+y2*v3+y3*v2H+y4*v3)+r*v2H*(y1*y4-y2*y3);
y3p=y3*(y1*v1+y2*v2H+y3*v1+y4*v2R)+r*v2H*(y1*y4-y2*y3);
y4p=y4*(y1*v2H+y2*v3+y3*v2R+y4*v3)-r*v2H*(y1*y4-y2*y3);

totx=Factor[x1p+x2p+x3p+x4p]

toty=Factor[y1p+y2p+y3p+y4p]

x1prime=(1-m1)*(x1p/totx)+m1*(y1p/toty);
x2prime=(1-m1)*(x2p/totx)+m1*(y2p/toty);
x3prime=(1-m1)*(x3p/totx)+m1*(y3p/toty);
x4prime=(1-m1)*(x4p/totx)+m1*(y4p/toty);

y1prime=(1-m2)*(y1p/toty)+m2*(x1p/totx);
y2prime=(1-m2)*(y2p/toty)+m2*(x2p/totx);
y3prime=(1-m2)*(y3p/toty)+m2*(x3p/totx);
y4prime=(1-m2)*(y4p/toty)+m2*(x4p/totx);

sub={w1->(1-s),
w2D->(1-h11*s),
w2H->(1-h12*s),
w2R->(1-h22*s),
w3->1,
v1->1,
v2D->(1-k11*t),
v2H->(1-k12*t),
v2R->(1-k22*t),
v3->1-t};

sub2={w1->(1-s),
w2D->(1-h11*s),
w2H->(1-h12*s),
w2R->(1-h22*s),
w3->1,
v1->1,
v2D->(1-k*t),
v2H->(1-k*t),
v2R->(1-k*t),
v3->1-t}/.s->1/.k->1/2/.h11->1/2;

Polymorphic equilibrium assuming weak
migration

The two equations that must equal zero at the equilibrium are:

eq1=Factor[x1-x1prime/.x3->0/.x4->0/.y3->0/.y4->0/.x2->1-x1/.y1->1-y2]

eq2=Factor[y2-y2prime/.x3->0/.x4->0/.y3->0/.y4->0/.x2->1-x1/.y1->1-y2]

Factor[Normal[Series[eq1/.
	m1->m1*small/.m2->m2*small/.x1->x1*small/.y2->y2*small,
	{small,0,1}]]]/.w3->1

Factor[Normal[Series[eq2/.
	m1->m1*small/.m2->m2*small/.x1->x1*small/.y2->y2*small,
	{small,0,1}]]]/.v1->1

The solutions to these equations give (5):

migx1=m1/(1-w2D);
migy2=m2/(1-v2D);

migx1alt=m1/(h11*s);
migy2alt=m2/(k*t);

Stability to invasion of modifier for r=0 and weak migration

matrix

sub4={w3->1,v1->1}

b1=Factor[Normal[Series[
x3prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small/.r->0,{small,0,1}]]]/.small->1;

row1={Factor[Coefficient[Expand[b1],ex3]],
Factor[Coefficient[Expand[b1],ex4]],
Factor[Coefficient[Expand[b1],ey3]],
Factor[Coefficient[Expand[b1],ey4]]}

b2=Factor[Normal[Series[
x4prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small/.r->0,{small,0,1}]]]/.small->1;

row2={Factor[Coefficient[Expand[b2],ex3]],
Factor[Coefficient[Expand[b2],ex4]],
Factor[Coefficient[Expand[b2],ey3]],
Factor[Coefficient[Expand[b2],ey4]]}

b3=Factor[Normal[Series[
y3prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small/.r->0,{small,0,1}]]]/.small->1;

row3={Factor[Coefficient[Expand[b3],ex3]],
Factor[Coefficient[Expand[b3],ex4]],
Factor[Coefficient[Expand[b3],ey3]],
Factor[Coefficient[Expand[b3],ey4]]}

b4=Factor[Normal[Series[
y4prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small/.r->0,{small,0,1}]]]/.small->1;

row4={Factor[Coefficient[Expand[b4],ex3]],
Factor[Coefficient[Expand[b4],ex4]],
Factor[Coefficient[Expand[b4],ey3]],
Factor[Coefficient[Expand[b4],ey4]]}

matrix3={row1,row2,row3,row4}

submat1 describes the invasion of x3/y3.
submat2 describes the invasion of x4/y4.

submat1={{((1 - m1)*(w1*x1 + w2H*x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2),  
   (m1*(y1 + v2H*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
  {(m2*(w1*x1 + w2H*x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
   ((1 - m2)*(y1 + v2H*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}};

submat2={{((1 - m1)*(w2H*x1 + x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
   (m1*(v2H*y1 + v3*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
  {(m2*(w2H*x1 + x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2),
   ((1 - m2)*(v2H*y1 + v3*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}};

Eigenvalue of submatrix1

Factor[Det[lambda*IdentityMatrix[2]-submat1]]

Here we calculate the leading eigenvalue as 1+sm, solving for sm:

Factor[Normal[Series[%/.lambda->1+sm*small/.y1->1-y2/.x2->1-x1/.y2->y2*small/.x1->x1*small/.m1->m1*small/.m2->m2*small,
	{small,0,1}]]]/.small->1

Simplify[Flatten[Solve[%==0,sm]]]

Factor[sm/.%/.y2->migy2/.x1->migx1]

Simplify[%*y2/migy2]

This gives equation (6b)

Eigenvalue of submatrix2

Factor[Det[lambda*IdentityMatrix[2]-submat2]]

Here we calculate the leading eigenvalue as 1+sm, solving for sm:

Factor[Normal[Series[%/.lambda->1+sm*small/.y1->1-y2/.x2->1-x1/.y2->y2*small/.x1->x1*small/.m1->m1*small/.m2->m2*small,
	{small,0,1}]]]/.small->1

Simplify[Flatten[Solve[%==0,sm]]]

Factor[sm/.%/.y2->migy2/.x1->migx1]

Simplify[%*x1/migx1]

This gives equation (6a)

Stability to invasion of modifier for general r

matrix

migx1=m1/(1-w2D);
migy2=m2/(1-v2D);

sub4={w3->1,v1->1}

b1=Normal[Series[
x3prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small,{small,0,1}]]/.small->1;

row1={Factor[Coefficient[Expand[b1],ex3]],
Factor[Coefficient[Expand[b1],ex4]],
Factor[Coefficient[Expand[b1],ey3]],
Factor[Coefficient[Expand[b1],ey4]]};

b2=Normal[Series[
x4prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small,{small,0,1}]]/.small->1;

row2={Factor[Coefficient[Expand[b2],ex3]],
Factor[Coefficient[Expand[b2],ex4]],
Factor[Coefficient[Expand[b2],ey3]],
Factor[Coefficient[Expand[b2],ey4]]};

b3=Normal[Series[
y3prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small,{small,0,1}]]/.small->1;

row3={Factor[Coefficient[Expand[b3],ex3]],
Factor[Coefficient[Expand[b3],ex4]],
Factor[Coefficient[Expand[b3],ey3]],
Factor[Coefficient[Expand[b3],ey4]]};

b4=Normal[Series[
y4prime/.sub4/.
x3->ex3*small/.x4->ex4*small/.
	y3->ey3*small/.y4->ey4*small,{small,0,1}]]/.small->1;

row4={Factor[Coefficient[Expand[b4],ex3]],
Factor[Coefficient[Expand[b4],ex4]],
Factor[Coefficient[Expand[b4],ey3]],
Factor[Coefficient[Expand[b4],ey4]]};

matrix3={row1,row2,row3,row4}

wbx=Factor[totx/.sub4/.x3->0/.x4->0/.y3->0/.y4->0]

wby=Factor[toty/.sub4/.x3->0/.x4->0/.y3->0/.y4->0]

Evaluated at lambda=1 assuming weak migration

We need the following to be negative for invasion (the characteristic polynomial at lambda = 1):

Det[IdentityMatrix[4]-matrix3]/.y1->1-y2/.x2->1-x1;

Factor[Normal[Series[%/.y2->y2*small/.x1->x1*small/.m1->m1*small/.m2->m2*small,
	{small,0,2}]]]/.small->1;

fac=Factor[%/.y1->1-y2/.x2->1-x1/.y2->migy2/.x1->migx1]

If the modifier increases fitness in one patch but has no effect in the other, then the following will be negative and invasion will occur:

Factor[fac/.w2H->w2D]

Factor[fac/.v2H->v2D]

The sum of these plus a cross-product term equals fac:

Factor[fac-%-%%]

For a weak modifier, we can write the value of the characteristic polynomial at lambda = 1 as:

-(((m1 m2 r v2D (w2H-w2D) (1-w2D(1-2 r)))/(1-w2D))+((m1 m2 r w2D(v2H-v2D) (1-v2D(1-2 r)) )/(1-v2D)));

Factor[Normal[Series[(%-fac/.v2H->v2D+epsV*small/.w2H->w2D+epsW*small),{small,0,1}]]/.small->1]

Ψ

Thus if equation (7) is positive, the value of the char poly at one is negative, and invasion is guaranteed.

Ψ=(w2H-w2D) v2D(1-v2D) (1-w2D(1-2 r))+(v2H-v2D) w2D(1-w2D)(1-v2D(1-2 r));

-(((m1 m2 r v2D (w2H-w2D) (1-w2D(1-2 r)))/(1-w2D))+((m1 m2 r w2D(v2H-v2D) (1-v2D(1-2 r)) )/(1-v2D)))/Ψ//Factor

Leading eigenvalue assuming weak migration and a weak modifier r>0

migx1=m1/(1-w2D);
migy2=m2/(1-v2D);

matrix4=lam*IdentityMatrix[4]-matrix3;

cp=Factor[Det[matrix4]];

To perform the stability analysis while keeping the exact relationships between the variables and parameters, we solve the equilibrium equations for m1 and m2 (these equations are linear in these parameters, making this requirement at the equilibrium simple to write):

subm=Flatten[Simplify[Solve[{eq1==0,eq2==0},{m1,m2}]]/.sub4]

Assuming weak migration and a weak modifier, the characteristic polynomial is approximately:

Factor[Normal[Series[cp/.y1->1-y2/.x2->1-x1/.subm/.lam->1+olam*small/.x1->x1*small/.y2->y2*small/.v2H->v2D+epsV*small/.w2H->w2D+epsW*small,
	{small,0,1}]]/.small->1]

We need the next order:

Factor[Normal[Series[cp/.y1->1-y2/.x2->1-x1/.subm/.lam->1+olam*small/.x1->x1*small/.y2->y2*small/.v2H->v2D+epsV*small/.w2H->w2D+epsW*small,
	{small,0,2}]]/.small->1]

Simplify[olam/.Solve[%==0,olam]]

The above indicates that olam = 0 (the parenthetical factor does not involve the modifier’s effects, epsW and epsV, and so is irrrelevant), and the eigenvalue then needs to be written allowing even smaller order terms, lam->1+olam2*(small^2):

Factor[Normal[Series[cp/.y1->1-y2/.x2->1-x1/.subm/.lam->1+olam2*(small^2)/.x1->x1*small/.y2->y2*small/.v2H->v2D+epsV*small/.w2H->w2D+epsW*small,
	{small,0,3}]]/.small->1]

This gives the selection coefficient on the modifier:

Collect[olam2/.Flatten[Solve[%==0,olam2]],{epsV,epsW},Simplify]

which can be written as equation (8), recognizing that x1 (1-w2D) is m1 at the equilibrium and y2 (1-v2D) is m2:

sol2=((x1 y2 Ψ)/(x1 (1-w2D) w2D(1-v2D(1-r))+y2 (1-v2D) v2D(1-w2D(1-r))));
%%-sol2/.epsW->(w2H-w2D)/.epsV->(v2H-v2D)//Factor

As r gets increasingly large (relative to 1-v2D and 1-w2D), the above approaches:

Factor[Limit[sol2,r->Infinity]]

which can be written as equation (9):

(((w2H-w2D) 2 x1 m2)/(m1+m2))+(((v2H-v2D) 2 y2 m1)/(m1+m2))/.m1->x1 (1-w2D) /.m2->y2 (1-v2D);
%-%%//Factor

Behavior with r

Factor[D[sol2,r]]

For modifiers that increase heterozygous fitness, the eigenvalue now increases with r if m1 and m2 are similar,
but it can decrease with r if migration rates are very different.

Collect[D[sol2,r]/.v2H->v2D+epsV/.w2H->w2D+epsW/.x1->migx1/.y2->migy2/.m1->m/.m2->m,{epsW,epsV},Factor]

Here is an example where the rates can go down with r:

Collect[D[sol2,r]/.v2H->v2D+epsV/.w2H->w2D+epsW/.x1->migx1/.y2->migy2/.m1->m/.m2->4m/.v2D->1/2/.w2D->3/4/.epsW->0,{epsW,epsV},Factor]

For graphing

General Definitions

For the numerical evaluation of the stability
matrix (which doesn't assume migration is small),
the precise allele frequency needs to be calculated:

Clear[allele]
allele[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	allele[h11,h12,s,k11,k12,t,m1,m2,r]=
If[Evaluate[x1/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]<0.999999 &&
   Evaluate[y2/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]<0.999999 &&
   Evaluate[x1/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]>0 &&
   Evaluate[y2/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]>0,
   FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}],
   FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.501},{y2,0.501}]]

Clear[eig2]
eig2[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	eig2[h11,h12,s,k11,k12,t,m1,m2,r]=
	({{((-1 + m1)*(-(w1*x1) - w2H*x2 + r*w2H*x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((1 - m1)*r*w2H*x1)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m1*(y1 + v2H*y2 - r*v2H*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   (m1*r*v2H*y1)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
 
  {((1 - m1)*r*w2H*x2)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((-1 + m1)*(-(w2H*x1) + r*w2H*x1 - x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m1*r*v2H*y2)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   (m1*(v2H*y1 - r*v2H*y1 + v3*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
 
  {(m2*(w1*x1 + w2H*x2 - r*w2H*x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m2*r*w2H*x1)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((-1 + m2)*(-y1 - v2H*y2 + r*v2H*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   ((1 - m2)*r*v2H*y1)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
 
  {(m2*r*w2H*x2)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m2*(w2H*x1 - r*w2H*x1 + x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((1 - m2)*r*v2H*y2)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   ((-1 + m2)*(-(v2H*y1) + r*v2H*y1 - v3*y2))/
   (y1^2 + 2*v2D*y1*y2 + v3*y2^2)}}/.
   y1->1-y2/.x2->1-x1/.
   allele[h11,h12,s,k11,k12,t,m1,m2,r]/.{w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t})

Clear[eig3]
eig3[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	eig3[h11,h12,s,k11,k12,t,m1,m2,r]=
	-1+Max[Re[Eigenvalues[eig2[h11,h12,s,k11,k12,t,m1,m2,r]]]]

allele[0.5,0.5,0.5,0.1,0.1,0.1,0.001,0.001,0.1]

eig3[0.5,0.1,0.5,0.1,0.1,0.1,0.001,0.001,0.1]

eig3[0.5,0.5,0.5,0.1,0.1,0.1,0.001,0.001,0.1]

Clear[solutr0]
solutr0[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	Max[{(m2*(v2D - v2H))/(-1 + v2D),
		(m1*(w2D - w2H))/(-1 + w2D)}/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t}]

Clear[solutrL]
solutrL[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	(w2H-w2D)*2*m1/(h11*s)*m2/(m1+m2)+	(v2H-v2D)*2*m2/(k11*t)*m1/(m1+m2)/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t}

Clear[solut]
solut[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	((m1/(1-w2D)*m2/(1-v2D)*((v2H-v2D)*(1 - v2D + 2*r*v2D)*(1-w2D)*w2D + 
 
      (w2H-w2D)*(1-v2D)*v2D*(1 - w2D + 2*r*w2D)))/
 
  (m1*(1-v2D+r*v2D)*w2D+m2*v2D*(1-w2D+r*w2D))/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t});

Clear[solut2]
solut2[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	(x1*y2*((v2H-v2D) (1-v2H+2 r v2H) (1-w2H) w2H+(1-v2H) v2H (w2H-w2D) (1-w2H+2 r w2H))/(v2H(1-v2H) (1-w2H+r w2H) y2+(1-v2H+r v2H) w2H(1- w2H) x1))/.
		{x1->(m1/(1-w2D)),y2->(m2/(1-v2D))}/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t};

solut[0.5,0.5,0.5,0.1,0.1,0.1,0.001,0.001,0.1]

solut[0.5,0.1,0.5,0.1,0.1,0.1,0.001,0.001,0.1]

solut2[0.5,0.1,0.5,0.1,0.1,0.1,0.001,0.001,0.1]

For comparison, the small-scale patchiness result was:

Clear[eigen3]
eigen3[h11_,h12_,s_,k11_,k12_,t_,a_,r_]:=
	eigen3[h11,h12,s,k11,k12,t,a,r]=
	If[
	a*s - a*h11*s - k11*t + a*k11*t>0 &&
	-(a*h11*s) + t - a*t - k11*t + a*k11*t>0,
	Max[eigen2[h11,h12,s,k11,k12,t,a,r][[1]],
		eigen2[h11,h12,s,k11,k12,t,a,r][[2]]],
	0]

Plots with r

Figure 5:
h11=0.5, h12=0.4, s=0.5
k11=0.1, k12=0.1, t=0.1
m1=0.005 m2 =0.001  (red)
vs
m1=0.001 m2 =0.005

Show[
ListPlot[Table[{r,
	eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.005,0.001,r]},
	{r,0,0.5,0.01}],Joined->True,PlotStyle->Red],
Plot[solut[0.5,0.4,0.5,0.1,0.1,0.1,0.005,0.001,r],
	{r,0,0.5},PlotStyle->{Red,Dashing[{0.02,0.04}]}],
ListPlot[Table[{r,
	solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.005,0.001,r]},
	{r,0,0,0.01}],PlotStyle->{Red,PointSize[0.03]}],
ListPlot[Table[{r,
	solutrL[0.5,0.4,0.5,0.1,0.1,0.1,0.005,0.001,r]},
	{r,0.5,0.5,0.01}],PlotStyle->{Red,PointSize[0.03]}],
ListPlot[Table[{r,
	eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.001,0.005,r]},
	{r,0,0.5,0.01}],Joined->True],
Plot[solut[0.5,0.4,0.5,0.1,0.1,0.1,0.001,0.005,r],
	{r,0,0.5},PlotStyle->Dashing[{0.02,0.04}]],
ListPlot[Table[{r,
	solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.001,0.005,r]},
	{r,0,0,0.01}],PlotStyle->PointSize[0.03]],
ListPlot[Table[{r,
	solutrL[0.5,0.4,0.5,0.1,0.1,0.1,0.001,0.005,r]},
	{r,0.5,0.5,0.01}],PlotStyle->PointSize[0.03]],
PlotRange->{{-0.02,0.52},{0,0.0016}},
AxesOrigin->{0,0}]

Plots with m

Figure 4a:
h11=0.5, h12=0.4, s=0.5
k11=0.1, k12=0.1, t=0.5
m1=m2, r =0.2

Show[
ListLogLinearPlot[Join[
	Table[{m,
	eig3[0.5,0.4,0.5,0.1,0.1,0.5,m,m,0.2]},
	{m,0.001,0.01,0.001}],Table[{m,
	eig3[0.5,0.4,0.5,0.1,0.1,0.5,m,m,0.2]},
	{m,0.01,0.5,0.01}]],Joined->True],
LogLinearPlot[solut[0.5,0.4,0.5,0.1,0.1,0.5,m,m,0.2],
	{m,0.001,0.5},PlotStyle->Dashing[{0.02,0.04}]],
ListLogLinearPlot[Table[{0.5,
	eigen3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.2]},
	{m,0,0,0.01}],PlotStyle->PointSize[0.03]],
PlotRange->{{Log[0.001],Log[0.5]},{0,0.02}},
AxesOrigin->{Log[0.001],0}]

Figure 4b:
h11=0.5, h12=0.4, s=0.5
k11=0.1, k12=0.1, t=0.1
m1=m2, r =0.2

Show[
ListLogLinearPlot[Join[
	Table[{m,
	eig3[0.5,0.4,0.5,0.1,0.1,0.1,m,m,0.2]},
	{m,0.001,0.01,0.001}],Table[{m,
	eig3[0.5,0.4,0.5,0.1,0.1,0.1,m,m,0.2]},
	{m,0.01,0.5,0.01}]],Joined->True],
LogLinearPlot[solut[0.5,0.4,0.5,0.1,0.1,0.1,m,m,0.2],
	{m,0.001,0.5},PlotStyle->Dashing[{0.02,0.04}]],
ListLogLinearPlot[Table[{0.5,
	eigen3[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.2]},
	{m,0,0,0.01}],PlotStyle->PointSize[0.03]],
PlotRange->{{Log[0.001],Log[0.5]},{0,0.02}},
AxesOrigin->{Log[0.001],0}]

Role of m2/(m1+m2)

In the neutral case, what proportion of the final population will come originally from patch1?  Here we label the initial individuals from patch 1 as “1” and from patch 2 as “0” to track which ones come from patch 1 originally.

RSolve[{pA1[t+1]==(1-m1)*pA1[t]+m1*pA2[t],pA2[t+1]==(1-m2)*pA2[t]+m2*pA1[t],pA1[0]==1,pA2[0]==0},{pA1[t],pA2[t]},t]

pA1 and pA2 are both approaching m2/(m1 + m2).  Thus if m2 is relatively high, then a very
large fraction of the ultimate population will come from patch 1 (and vice versa).

Clear[pA1,pa1,pA2,pa2,m1,m2]
pA1[t_,m1_,m2_]:=pA1[t,m1,m2]=(1-m1)*pA1[t-1,m1,m2]+m1*pA2[t-1,m1,m2]
pa1[t_,m1_,m2_]:=pa1[t,m1,m2]=(1-m1)*pa1[t-1,m1,m2]+m1*pa2[t-1,m1,m2]
pA2[t_,m1_,m2_]:=pA2[t,m1,m2]=(1-m2)*pA2[t-1,m1,m2]+m2*pA1[t-1,m1,m2]
pa2[t_,m1_,m2_]:=pa2[t,m1,m2]=(1-m2)*pa2[t-1,m1,m2]+m2*pa1[t-1,m1,m2]

pA1[0,m1_,m2_]=1;
pa1[0,m1_,m2_]=0;
pA2[0,m1_,m2_]=0;
pa2[0,m1_,m2_]=1;

The rate of approach is:

Simplify[pA1[3,m1,m2]-m2/(m1+m2)]

Simplify[pA1[4,m1,m2]-m2/(m1+m2)]

Simplify[pA1[5,m1,m2]-m2/(m1+m2)]

Simplify[pA2[5,m1,m2]-m2/(m1+m2)]

Simulations (for exploration)

Getting the equations

x1p=x1*(x1*w1+x2*w2D+x3*w1+x4*w2H)-r*w2H*(x1*x4-x2*x3);
x2p=x2*(x1*w2D+x2*w3+x3*w2H+x4*w3)+r*w2H*(x1*x4-x2*x3);
x3p=x3*(x1*w1+x2*w2H+x3*w1+x4*w2R)+r*w2H*(x1*x4-x2*x3);
x4p=x4*(x1*w2H+x2*w3+x3*w2R+x4*w3)-r*w2H*(x1*x4-x2*x3);

y1p=y1*(y1*v1+y2*v2D+y3*v1+y4*v2H)-r*v2H*(y1*y4-y2*y3);
y2p=y2*(y1*v2D+y2*v3+y3*v2H+y4*v3)+r*v2H*(y1*y4-y2*y3);
y3p=y3*(y1*v1+y2*v2H+y3*v1+y4*v2R)+r*v2H*(y1*y4-y2*y3);
y4p=y4*(y1*v2H+y2*v3+y3*v2R+y4*v3)-r*v2H*(y1*y4-y2*y3);

totx=Factor[x1p+x2p+x3p+x4p]

toty=Factor[y1p+y2p+y3p+y4p]

x1prime=(1-m1)*(x1p/totx)+m1*(y1p/toty);
x2prime=(1-m1)*(x2p/totx)+m1*(y2p/toty);
x3prime=(1-m1)*(x3p/totx)+m1*(y3p/toty);
x4prime=(1-m1)*(x4p/totx)+m1*(y4p/toty);

y1prime=(1-m2)*(y1p/toty)+m2*(x1p/totx);
y2prime=(1-m2)*(y2p/toty)+m2*(x2p/totx);
y3prime=(1-m2)*(y3p/toty)+m2*(x3p/totx);
y4prime=(1-m2)*(y4p/toty)+m2*(x4p/totx);

sub={w1->(1-s),
w2D->(1-h11*s),
w2H->(1-h12*s),
w2R->(1-h22*s),
w3->1,
v1->1,
v2D->(1-k11*t),
v2H->(1-k12*t),
v2R->(1-k22*t),
v3->1-t};

sub2={x1->x1[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],x2->x2[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],
		x3->x3[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],x4->x4[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],
 	 y1->y1[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],
		y2->y2[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],y3->y3[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],y4->y4[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22]};

Simplify[{x1prime,x2prime,x3prime,x4prime,y1prime,y2prime,y3prime,y4prime}/.sub/.sub2]

The equations [MUST RE-ENTER FOR EACH NEW STARTING POSITION]

sub={w1->(1-s),
w2D->(1-h11*s),
w2H->(1-h12*s),
w2R->(1-h22*s),
w3->1,
v1->1,
v2D->(1-k11*t),
v2H->(1-k12*t),
v2R->(1-k22*t),
v3->1-t};

sub2={x1->x1[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],x2->x2[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],
		x3->x3[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],x4->x4[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],
 	 y1->y1[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],
		y2->y2[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],y3->y3[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22],y4->y4[g-1,m1,m2,s,h11,h12,h22,t,k11,k12,k22]};

Clear[x1,x2,x3,x4,y1,y2,y3,y4]
x1[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x1[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=((1-m1) (x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h11 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h12 s) x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+h12 s) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+(m1 (y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k11 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k12 t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+k12 t) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2);

x2[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x2[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=((1-m1) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((1-h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h12 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+h12 s) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+(m1 (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((1-k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k12 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+k12 t) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2);

x3[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=
  x3[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=((1-m1) (x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h22 s) x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+h12 s) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+(m1 (y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((-1+r) (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-k22 t y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+r (1-k12 t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2); 

x4[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=
  x4[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=((1-m1) (x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((-1+r) (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-h22 s x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((r-h12 r s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+(m1 (y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((1-k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+k12 t) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2); 

y1[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=
  y1[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=(m2 (x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h11 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h12 s) x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+h12 s) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+((1-m2) (y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k11 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k12 t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+k12 t) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2); 

y2[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=
  y2[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=(m2 (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((1-h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h12 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+h12 s) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+((1-m2) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((1-k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k12 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+k12 t) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2); 

y3[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=
  y3[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=(m2 (x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-h22 s) x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+r (-1+h12 s) (x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+((1-m2) (y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((-1+r) (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-k22 t y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] (y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+r (1-k12 t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2); 

y4[g_,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=
  y4[g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]=(m2 (x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((-1+r) (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-h22 s x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])+x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((r-h12 r s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(-(-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h11 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h12 s) x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2-2 (-1+h12 s) x1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+2 x2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+h22 s) x3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((x4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2))+((1-m2) (y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] ((1-k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+(1-k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])-r (-1+k12 t) (y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])))/(((y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k11 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2+2 y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k12 t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]+((y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r])^2)-2 (-1+k12 t) y1[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+t) y2[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-2 (-1+k22 t) y3[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r] y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]-(-1+t) y4[-1+g,m1,m2,s,h11,h12,h22,t,k11,k12,k22,r]^2);

General Definitions

For the numerical evaluation of the stability
matrix (which doesn't assume migration is small),
the precise allele frequency needs to be calculated:

Clear[allele]
allele[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	allele[h11,h12,s,k11,k12,t,m1,m2,r]=
If[Evaluate[x1/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]<0.999999 &&
   Evaluate[y2/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]<0.999999 &&
   Evaluate[x1/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]>0 &&
   Evaluate[y2/.FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}]]>0,
   FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.5*m1/(h11*s)},{y2,0.5*m2/(k11*t)}],
   FindRoot[{
x1==((1 - m1)*x1*(1 - h11*s - s*x1 + h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  (m1*(1 - y2)*(1 - k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2),
y2==(m2*(1 - x1)*(1 - h11*s*x1))/
 
   (1 - 2*h11*s*x1 - s*x1^2 + 2*h11*s*x1^2) + 
 
  ((1 - m2)*y2*(1 - k11*t - t*y2 + k11*t*y2))/
 
   (1 - 2*k11*t*y2 - t*y2^2 + 2*k11*t*y2^2)},
   {x1,0.501},{y2,0.501}]]

Clear[eig2]
eig2[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	eig2[h11,h12,s,k11,k12,t,m1,m2,r]=
	({{((-1 + m1)*(-(w1*x1) - w2H*x2 + r*w2H*x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((1 - m1)*r*w2H*x1)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m1*(y1 + v2H*y2 - r*v2H*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   (m1*r*v2H*y1)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
 
  {((1 - m1)*r*w2H*x2)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((-1 + m1)*(-(w2H*x1) + r*w2H*x1 - x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m1*r*v2H*y2)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   (m1*(v2H*y1 - r*v2H*y1 + v3*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
 
  {(m2*(w1*x1 + w2H*x2 - r*w2H*x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m2*r*w2H*x1)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((-1 + m2)*(-y1 - v2H*y2 + r*v2H*y2))/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   ((1 - m2)*r*v2H*y1)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2)}, 
 
  {(m2*r*w2H*x2)/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   (m2*(w2H*x1 - r*w2H*x1 + x2))/(w1*x1^2 + 2*w2D*x1*x2 + x2^2), 
 
   ((1 - m2)*r*v2H*y2)/(y1^2 + 2*v2D*y1*y2 + v3*y2^2), 
 
   ((-1 + m2)*(-(v2H*y1) + r*v2H*y1 - v3*y2))/
   (y1^2 + 2*v2D*y1*y2 + v3*y2^2)}}/.
   y1->1-y2/.x2->1-x1/.
   allele[h11,h12,s,k11,k12,t,m1,m2,r]/.{w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t})

Clear[eig3]
eig3[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	eig3[h11,h12,s,k11,k12,t,m1,m2,r]=
	-1+Max[Re[Eigenvalues[eig2[h11,h12,s,k11,k12,t,m1,m2,r]]]]

Clear[solutr0]
solutr0[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	Max[{(m2*(v2D - v2H))/(-1 + v2D),
		(m1*(w2D - w2H))/(-1 + w2D)}/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t}]

Clear[solutrL]
solutrL[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	(w2H-w2D)*2*m1/(h11*s)*m2/(m1+m2)+	(v2H-v2D)*2*m2/(k11*t)*m1/(m1+m2)/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t}

Clear[solut]
solut[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	((m1/(1-w2D)*m2/(1-v2D)*((v2H-v2D)*(1 - v2D + 2*r*v2D)*(1-w2D)*w2D + 
 
      (w2H-w2D)*(1-v2D)*v2D*(1 - w2D + 2*r*w2D)))/
 
  (m1*(1-v2D+r*v2D)*w2D+m2*v2D*(1-w2D+r*w2D))/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t});

Clear[solut2]
solut2[h11_,h12_,s_,k11_,k12_,t_,m1_,m2_,r_]:=
	(x1*y2*((v2H-v2D) (1-v2H+2 r v2H) (1-w2H) w2H+(1-v2H) v2H (w2H-w2D) (1-w2H+2 r w2H))/(v2H(1-v2H) (1-w2H+r w2H) y2+(1-v2H+r v2H) w2H(1- w2H) x1))/.
		{x1->(m1/(1-w2D)),y2->(m2/(1-v2D))}/.
   {w1->(1-s),
	w2D->(1-h11*s),
	w2H->(1-h12*s),
	w2R->(1-h22*s),
	w3->1,
	v1->1,
	v2D->(1-k11*t),
	v2H->(1-k12*t),
	v2R->(1-k22*t),
	v3->1-t};

One-locus start

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=m1/(h11*s);
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=1-m1/(h11*s);
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=1-m2/(k11*t);
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=m2/(k11*t);
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;

Show[
	ListPlot[Table[x1[gg,0.001,0.001,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,200}],Joined->True],
	ListPlot[Table[y2[gg,0.001,0.001,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,200}],Joined->True],
PlotRange->{Automatic,{0,0.1}},
AxesOrigin->{0,0}
]

Show[
	ListPlot[Table[x1[gg,0.0001,0.0001,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,200}],Joined->True],
	ListPlot[Table[y2[gg,0.0001,0.0001,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,200}],Joined->True],
PlotRange->{Automatic,{0,0.01}},
AxesOrigin->{0,0}
]

Confirming no change occurs if started at the equilibrium:

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0;

Show[
	ListPlot[Table[x1[gg,0.001,0.001,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,200}],Joined->True],
	ListPlot[Table[y2[gg,0.001,0.001,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.1],{gg,1,200}],Joined->True],
PlotRange->{Automatic,{0,0.1}},
AxesOrigin->{0,0}
]

Two-locus start - Strong selection in patch 2 (low mig)

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];

Looking at r=0

For[gg=1,gg<7001,gg++,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True],		ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True]
]

Because r=0 and the mA haplotype does not have a selective advantage,
the ma haplotype (4)  invades to replace the Ma haplotype (2), but the
most frequent A haplotype remains MA (1) not mA (3).

Show[
	ListPlot[Table[{gg,x1[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True],		ListPlot[Table[{gg,y1[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True]
]

solutr0[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]

eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True,PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,7200},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,7200},PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,70}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.04,0.01}]]
]

Looking at r=0.5

For[gg=1,gg<7001,gg++,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],		ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]/.p0->0.01,{t,0,6000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]/.p0->0.01,{t,0,6000},PlotStyle->Dashing[{0.04,0.01}]]
]

solut[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]

solut2[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]

eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]

Looking at r=0.01

For[gg=1,gg<7001,gg++,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,6000,100}],Joined->True],		ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.06,0.06}],Thickness[0.01]}],
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.01,0.01}],Thickness[0.01]}],
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,7000,100}],Joined->True,PlotStyle->Thickness[0.01]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]/.p0->0.01,{t,0,7000}]
]

Show[
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.06,0.06}],Thickness[0.01]}],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.01,0.01}],Thickness[0.01]}],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,7000,100}],Joined->True,PlotStyle->Thickness[0.01]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]/.p0->0.01,{t,0,7000}]
]

Show[
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.04,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]/.p0->0.01,{t,0,7000}]
]

solut[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]

solut2[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]

eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.01]

eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.0]

eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.01,0.01,0.5]

Two-locus start - Weak selection in patch 2 (low mig)

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];

Looking at r=0

For[gg=1,gg<7001,gg++,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True],		ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True]
]

Because r=0 and the mA haplotype does not have a selective advantage,
the ma haplotype (4)  invades to replace the Ma haplotype (2), but the
most frequent A haplotype remains MA (1) not mA (3).

Show[
	ListPlot[Table[{gg,x1[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True],		ListPlot[Table[{gg,y1[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True]
]

solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]

eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,7200},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,7200},PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,70}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.04,0.01}]]
]

Looking at r=0.5

For[gg=1,gg<7001,gg++,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],		ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,6000,100}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]/.p0->0.01,{t,0,6000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]/.p0->0.01,{t,0,6000},PlotStyle->Dashing[{0.04,0.01}]]
]

solut[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]

solut2[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]

eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]

Looking at r=0.01

For[gg=1,gg<7001,gg++,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,6000,100}],Joined->True],		ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,6000,100}],Joined->True,PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.04,0.01}]]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.06,0.06}],Thickness[0.01]}],
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.01,0.01}],Thickness[0.01]}],
	ListPlot[Table[{gg,x3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+x4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,7000,100}],Joined->True,PlotStyle->Thickness[0.01]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]/.p0->0.01,{t,0,7000}]
]

Show[
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.06,0.06}],Thickness[0.01]}],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,7000,100}],Joined->True,PlotStyle->{Dashing[{0.01,0.01}],Thickness[0.01]}],
	ListPlot[Table[{gg,y3[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+y4[gg,0.01,0.01,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,7000,100}],Joined->True,PlotStyle->Thickness[0.01]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]/.p0->0.01,{t,0,7000}]
]

Show[
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.04,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]/.p0->0.01,{t,0,7000},PlotStyle->Dashing[{0.01,0.01}]]
]

solut[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]

solut2[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]

eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.01]

eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.0]

eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.01,0.01,0.5]

Two-locus start - Strong selection in patch 2 (high mig)

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];

Looking at r=0

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.03,0.06}],PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.005,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.0]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.05,0.05}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.0]/.p0->0.01,{t,0,1200},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All],
	PlotRange->{{0,1000},{0,1}}
]

{x1[0,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0],x2[0,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0],
	y1[0,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0],y2[0,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0]}

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0]]

{x3[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0],x4[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0],
	y3[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0],y4[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.3,0.1,0.1,0.1,0.0]}

{x1[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],x2[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],
	y1[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],y2[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]}

{x3[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],x4[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],
	y3[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],y4[0,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]}

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]]

{x1[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],x2[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],
	y1[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],y2[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]}

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]]

{x3[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],x4[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],
	y3[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0],y4[1000,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]}

This result is different than the weak migration case and more similar to the Dempster
model (small-scale patch model).  The bump is missing in this r=0 case because the two
haplotypes are nearly equal in frequency.  If we reduce t to 0.25, then the initial allele frequency
of A is lower and we see the bump again (next section) since the ma haplotype rises more slowly.

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,70}],Joined->True],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.05,0.05}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.04,0.01}]]
]

Looking at r=0.5

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.03,0.06}],PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.005,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.5]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.05,0.05}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.5]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All]
]

Weak migration doesn't (and shouldn't) give a good estimate for selection in this case,
but the exact eigenvalue solution still does.  Note that high migration erases the difference
between the dynamics in patch 1 and patch 2.

Looking at r=0.01

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.03,0.06}],PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.005,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.01]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.05,0.05}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.01]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->{Dashing[{0.06,0.06}],Thickness[0.01]}],
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.01]},{gg,1,1000,10}],Joined->True,PlotStyle->{Dashing[{0.01,0.01}],Thickness[0.01]}],
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.5,0.1,0.1,0.1,0.5]},{gg,1,1000,10}],Joined->True,PlotStyle->Thickness[0.01]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.0]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.01]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.5,0.5,0.5,0.5]/.p0->0.01,{t,0,1000}]
]

Two-locus start - Medium selection in patch 2 (high mig)

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];

Looking at r=0

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.03,0.06}],PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.005,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.05,0.05}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01*x1/.allele[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0],{t,0,1000},PlotStyle->Dashing[{0.01,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01,{t,0,1200},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All],
	PlotRange->{{0,1000},{0,1}}
]

This shows that the selection coefficient obtained from the eigenvalue does 
adequately describe the increase in frequency of the modifier, but only
when it's realized that selection is more strongly favoring the mA haplotype,
which is starting at a much lower frequency than 0.01.

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,70}],Joined->True,PlotRange->All],
	ListPlot[Table[{gg,x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,70}],Joined->True,PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.05,0.05}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01*x1/.allele[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0],{t,0,72},PlotStyle->Dashing[{0.05,0.05}]]
]

Looking at r=0.5

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.5]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.5]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.5]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.03,0.06}],PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.5]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.5]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.005,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.5]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.05,0.05}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.5]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All]
]

Weak migration doesn't (and shouldn't) give a good estimate for selection in this case,
but the exact eigenvalue solution still does.  Note that high migration erases the difference
between the dynamics in patch 1 and patch 2.

Looking at r=0.01

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.01]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.01]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.01]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.03,0.06}],PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.01]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.01]},{gg,1,1000,10}],Joined->True,PlotStyle->Dashing[{0.005,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.01]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.05,0.05}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.01]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->{Dashing[{0.06,0.06}],Thickness[0.01]},PlotRange->All],
		ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.01]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->{Dashing[{0.01,0.01}],Thickness[0.01]},PlotRange->All],
		ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.5]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.25,0.1,0.1,0.1,0.0]},{gg,1,1000,10}],Joined->True,PlotStyle->{Thickness[0.01]}],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0]/.p0->0.01*x1/.allele[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0],{t,0,1000},PlotStyle->Dashing[{0.06,0.06}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.01]/.p0->0.01,{t,0,1000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.01]/.p0->0.01*x1/.allele[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.0],{t,0,1000},PlotStyle->Dashing[{0.01,0.01}]],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.5]/.p0->0.01,{t,0,1000}]
]

eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.01]

eig3[0.5,0.4,0.5,0.1,0.1,0.25,0.5,0.5,0.5]

Two-locus start - Weak selection in patch 2 (high mig)

x1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
x3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
x4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-x1/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y1[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y2[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.99*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];
y3[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*(1-y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r]);
y4[0,m1_,m2_,s_,h11_,h12_,h22_,t_,k11_,k12_,k22_,r_]:=0.01*y2/.allele[h11,h12,s,k11,k12,t,m1,m2,r];

We have to use different starting frequencies, because the polymorphism is lost.

allele[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]

Looking at r=0

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]]

eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]

allele[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]

x1[0,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,100}],Joined->True,PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]/.p0->0.01,{t,0,100},PlotStyle->Dashing[{0.01,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]/.p0->0.01,{t,0,100},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All]
]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.0]},{gg,1,70}],Joined->True,PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.5,0.5}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solutr0[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.0]/.p0->0.01,{t,0,72},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All]
]

Looking at r=0.5

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,100}],Joined->True,PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.5]},{gg,1,100}],Joined->True,PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.5]/.p0->0.01,{t,0,100},PlotStyle->Dashing[{0.01,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.5]/.p0->0.01,{t,0,100},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All],
	PlotRange->{Automatic,{0,0.02}}
]

Looking at r=0.01

For[gg=1,gg<1001,gg++,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]]

Show[
	ListPlot[Table[{gg,x3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+x4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,20}],Joined->True,PlotRange->All],
	ListPlot[Table[{gg,y3[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]+y4[gg,0.5,0.5,0.5,0.5,0.4,0.3,0.1,0.1,0.1,0.1,0.01]},{gg,1,20}],Joined->True,PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->eig3[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.01]/.p0->0.01,{t,0,100},PlotStyle->Dashing[{0.01,0.01}],PlotRange->All],
	Plot[p0*Exp[s*t]/(1-p0+p0*Exp[s*t])/.s->solut[0.5,0.4,0.5,0.1,0.1,0.1,0.5,0.5,0.01]/.p0->0.01,{t,0,100},PlotStyle->Dashing[{0.04,0.01}],PlotRange->All],
	PlotRange->{Automatic,{0,0.02}}
]