Deleterious Mutations, Variable Epistatic Interactions, and the Evolution of Recombination

Sarah P. Otto & Marcus W. Feldman

Local stability analysis

HAPLOID MODEL

The recursions

sx1 := (1-mu)^2*x1
sx2 := (mu*(1-mu)*x1+(1-mu)*x2)
sx3 := (mu*(1-mu)*x1+(1-mu)*x3)
sx4 := (mu^2*x1+mu*x2+mu*x3+x4)
sx5 := (1-mu)^2*x5
sx6 := (mu*(1-mu)*x5+(1-mu)*x6)
sx7 := (mu*(1-mu)*x5+(1-mu)*x7)
sx8 := (mu^2*x5+mu*x6+mu*x7+x8)

newx1 := (sx1-r1*(sx1*sx4-sx2*sx3)-(r2+R-2*r2*R)*(sx1*sx6-sx2*sx5)-
	R*(sx1*sx7-sx3*sx5)-(R-r2*R)*(sx1*sx8-sx4*sx5)+
	r2*R*(sx3*sx6-sx2*sx7)+r2*(sx2*sx7-sx1*sx8))

newx2 := W1*(sx2+r1*(sx1*sx4-sx2*sx3)+(r2+R-2*r2*R)*(sx1*sx6-sx2*sx5)-
	R*(sx2*sx8-sx4*sx6)+(R-r2*R)*(sx3*sx6-sx2*sx7)-
	r2*R*(sx1*sx8-sx4*sx5)+r2*(sx1*sx8-sx2*sx7))

newx3 := W1*(sx3+r1*(sx1*sx4-sx2*sx3)-(r2+R-2*r2*R)*(sx3*sx8-sx4*sx7)+
	R*(sx1*sx7-sx3*sx5)-(R-r2*R)*(sx3*sx6-sx2*sx7)+
	r2*R*(sx1*sx8-sx4*sx5)-r2*(sx3*sx6-sx4*sx5))

newx4 := W2*(sx4-r1*(sx1*sx4-sx2*sx3)+(r2+R-2*r2*R)*(sx3*sx8-sx4*sx7)+
	R*(sx2*sx8-sx4*sx6)+(R-r2*R)*(sx1*sx8-sx4*sx5)-
	r2*R*(sx3*sx6-sx2*sx7)-r2*(sx4*sx5-sx3*sx6))

newx5 := (sx5-r3*(sx5*sx8-sx6*sx7)+(r2+R-2*r2*R)*(sx1*sx6-sx2*sx5)+
	R*(sx1*sx7-sx3*sx5)+(R-r2*R)*(sx1*sx8-sx4*sx5)-
	r2*R*(sx3*sx6-sx2*sx7)-r2*(sx4*sx5-sx3*sx6))

newx6 := W1*(sx6+r3*(sx5*sx8-sx6*sx7)-(r2+R-2*r2*R)*(sx1*sx6-sx2*sx5)+
	R*(sx2*sx8-sx4*sx6)-(R-r2*R)*(sx3*sx6-sx2*sx7)+r2*R*(sx1*sx8-sx4*sx5)-
	r2*(sx3*sx6-sx4*sx5))

newx7 := W1*(sx7+r3*(sx5*sx8-sx6*sx7)+(r2+R-2*r2*R)*(sx3*sx8-sx4*sx7)-
	R*(sx1*sx7-sx3*sx5)+(R-r2*R)*(sx3*sx6-sx2*sx7)-r2*R*(sx1*sx8-sx4*sx5)+
	r2*(sx1*sx8-sx2*sx7))

newx8 := W2*(sx8-r3*(sx5*sx8-sx6*sx7)-(r2+R-2*r2*R)*(sx3*sx8-sx4*sx7)-
	R*(sx2*sx8-sx4*sx6)-(R-r2*R)*(sx1*sx8-sx4*sx5)+r2*R*(sx3*sx6-sx2*sx7)+
	r2*(sx2*sx7-sx1*sx8))

Fitness:=(newx1+newx2+newx3+newx4+newx5+newx6+newx7+newx8);

Showing Remark 2

A1=Simplify[Wb^2*(x1*x4-x2*x3)-(newx4*newx1-newx2*newx3)/.{x5->0,x6->0,x7->0,x8->0}]

A2=Simplify[A1/.x1*x4->D+x2*x3/.x4->1-x2-x3-x1/.x3->x2
/.{x1->D+(1-pm)^2,x2->-D+pm*(1-pm)}]

A3=Factor[A2-(Wb^2-W2*(1-mu)^2*(1-r1))*D]

A3/Factor[((1-mu)*(mu*x1+x2)+r1*(1-mu)^2*D)^2/.x4->1-x2-x3-x1/.x3->x2
/.{x1->D+(1-pm)^2,x2->-D+pm*(1-pm)}]

Getting the  stability matrix

F1=Factor[Factor[Series[newx5/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row1={Factor[Coefficient[Expand[F1],x5]],Factor[Coefficient[Expand[F1],x6]],
Factor[Coefficient[Expand[F1],x7]],Factor[Coefficient[Expand[F1],x8]]}

F2=Factor[Factor[Series[newx6/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row2={Factor[Coefficient[Expand[F2],x5]],Factor[Coefficient[Expand[F2],x6]],
Factor[Coefficient[Expand[F2],x7]],Factor[Coefficient[Expand[F2],x8]]}

F3=Factor[Factor[Series[newx7/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row3={Factor[Coefficient[Expand[F3],x5]],Factor[Coefficient[Expand[F3],x6]],
Factor[Coefficient[Expand[F3],x7]],Factor[Coefficient[Expand[F3],x8]]}

F4=Factor[Factor[Series[newx8/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row4={Factor[Coefficient[Expand[F4],x5]],Factor[Coefficient[Expand[F4],x6]],
Factor[Coefficient[Expand[F4],x7]],Factor[Coefficient[Expand[F4],x8]]}

The  matrix

Row1={
(-1 + mu)^2*(1 - R*x2 - r2*x2 + 2*R*r2*x2 - R*x3 - R*x4 - r2*x4 + R*r2*x4), 
(-1 + mu)^2*(R*x1 + r2*x1 - 2*R*r2*x1 + r2*x3 - R*r2*x3),
(-1 + mu)^2*R*(x1 + r2*x2), 
(-1 + mu)^2*R*(1 - r2)*x1};

Row2={
(-1 + mu)*W1*(-mu - R*x2 + mu*R*x2 - r2*x2 + mu*r2*x2 + 2*R*r2*x2 -
	2*mu*R*r2*x2 + mu*R*x3 + mu*R*x4 - r2*x4 + mu*r2*x4 + R*r2*x4 - mu*R*r2*x4), 
(-1 + mu)*W1*(-1 + R*x1 - mu*R*x1 + r2*x1 - mu*r2*x1 - 2*R*r2*x1 +
	2*mu*R*r2*x1 + R*x3 + r2*x3 - mu*r2*x3 - R*r2*x3 + mu*R*r2*x3 + R*x4),
(-1 + mu)*R*W1*(-(mu*x1) - x2 + r2*x2 - mu*r2*x2), 
(-1 + mu)*R*W1*(-(mu*x1) - r2*x1 + mu*r2*x1 - x2)};

Row3={
(-1 + mu)*W1*(-mu + mu*R*x2 + mu*r2*x2 - 2*mu*R*r2*x2 - R*x3 + mu*R*x3 +
	mu*R*x4 + mu*r2*x4 - R*r2*x4 - mu*R*r2*x4),
(-1 + mu)*W1*(-(mu*R*x1) - mu*r2*x1 + 2*mu*R*r2*x1 - R*x3 - mu*r2*x3 +
	R*r2*x3 + mu*R*r2*x3),
(-1 + mu)*W1*(-1 + R*x1 - mu*R*x1 + R*x2 + r2*x2 - R*r2*x2 - mu*R*r2*x2 +
	R*x4 + r2*x4 - 2*R*r2*x4),
(-1 + mu)*W1*(-(mu*R*x1) - r2*x1 + R*r2*x1 + mu*R*r2*x1 - R*x3 - r2*x3 + 
	2*R*r2*x3)};

Row4={
W2*(mu^2 + mu*R*x2 - mu^2*R*x2 + mu*r2*x2 - mu^2*r2*x2 - 2*mu*R*r2*x2 +
	2*mu^2*R*r2*x2 + mu*R*x3 - mu^2*R*x3 + R*x4 - mu^2*R*x4 + mu*r2*x4 -
	mu^2*r2*x4 - R*r2*x4 + mu^2*R*r2*x4), 
W2*(mu - mu*R*x1 + mu^2*R*x1 - mu*r2*x1 + mu^2*r2*x1 + 2*mu*R*r2*x1 -
	2*mu^2*R*r2*x1 - mu*r2*x3 + mu^2*r2*x3 + R*r2*x3 - mu^2*R*r2*x3 +
	R*x4 - mu*R*x4), 
W2*(mu - mu*R*x1 + mu^2*R*x1 + r2*x2 - mu*r2*x2 - R*r2*x2 + mu^2*R*r2*x2 +
	R*x4 - mu*R*x4 + r2*x4 - mu*r2*x4 - 2*R*r2*x4 + 2*mu*R*r2*x4), 
W2*(1 - R*x1 + mu^2*R*x1 - r2*x1 + mu*r2*x1 + R*r2*x1 - mu^2*R*r2*x1 -
	R*x2 + mu*R*x2 - R*x3 + mu*R*x3 - r2*x3 + mu*r2*x3 + 2*R*r2*x3 -
	2*mu*R*r2*x3)};

StabMatrix=Factor[{Row1,Row2,Row3,Row4}/.{x4->1-x1-x2-x3}/.x3->x2];

SMatrix[mu_,W1_,W2_,r1_,r2_,R_,pm_,Ds_]=Factor[{Row1,Row2,Row3,Row4}/.{x4->1-x1-x2-x3}/.x3->x2]/.
	{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)};

Series of matrix manipulations

Here we will manipulate the above stability matrix to get it into a simpler format.

The value of the characteristic polynomial at lambda=1 is equal to:

m = Wb*IdentityMatrix[4] - StabMatrix;

where we have multiplied the characteristic polynomial by the mean fitness^4 (Wb^4).

tm = Transpose[m];

Multiply columns by x_i (need to divide final result by x1*x2*x3*x4)

m2 = Transpose[{tm[[1]]*x1,tm[[2]]*x2,tm[[3]]*x3,tm[[4]]*x4}];

tm2=Transpose[m2];

Add all columns to the first column

m3 = Transpose[{tm2[[1]]+tm2[[2]]+tm2[[3]]+tm2[[4]],
			tm2[[2]],tm2[[3]],tm2[[4]]}]/.x4->1-x2-x3-x1/.x3->x2/.
		{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)};

Make subsutitutions for Wb=newxi/xi in the first column

m4 = Factor[{{m3[[1]][[1]]/.Wb->Factor[(newx1/x1)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[1]][[2]],m3[[1]][[3]],m3[[1]][[4]]},
	{m3[[2]][[1]]/.Wb->Factor[(newx2/x2)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[2]][[2]],m3[[2]][[3]],m3[[2]][[4]]},
	{m3[[3]][[1]]/.Wb->Factor[(newx3/x3)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[3]][[2]],m3[[3]][[3]],m3[[3]][[4]]},
	{m3[[4]][[1]]/.Wb->Factor[(newx4/x4)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[4]][[2]],m3[[4]][[3]],m3[[4]][[4]]}}]

Transpose[m4];

Add row 1 to rows 2 and 3 and subtract it from row 4.

m5 = Factor[{m4[[1]],W1*m4[[1]]+m4[[2]],W1*m4[[1]]+m4[[3]],-W2*m4[[1]]+m4[[4]]}]

The stability matrix now equals a constant (C1) times a 3x3 matrix

C1=m5[[1,1]]

m6 = m5[[{2,3,4},{2,3,4}]]

tm6=Transpose[m6];

Divide the first column by x2 ,  the second column by x3, and the last column by x4.
(Need to multiply the final result by x2*x3*x4)

m7=Factor[Transpose[{tm6[[1]]/(Ds-pm+pm^2),tm6[[2]]/(Ds-pm+pm^2),tm6[[3]]/(Ds+pm^2)}]]

tm7=Transpose[m7];

Subtract column 1 from column 2.  Add column 1 to column 3.

m8=Factor[Transpose[{tm7[[1]],tm7[[2]]-tm7[[1]],tm7[[3]]+tm7[[1]]}]]

Add Row1*W2 +Row2*W2 +Row3*W1.
(Will need to divide the final result by W2)

m9={Factor[m8[[1]]*W2+m8[[2]]*W2+m8[[3]]*W1],m8[[2]],m8[[3]]}

Add Row2/W1 to Row3/W2.
(Will need to multiply the final result by W2)

m10={m9[[1]],m9[[2]],Factor[m9[[3]]/W2+m9[[2]]/W1]}

tm10=Transpose[m10];

Add the first colum to the second column.  Substract the third column from the first.

mFinal=Transpose[{Factor[tm10[[1]]-tm10[[3]]],Factor[tm10[[1]]+tm10[[2]]],tm10[[3]]}]

The determinant of this equals

Det1=Factor[Det[mFinal]];

which can be written as

Det2=Factor[((
		-R^2*(1-mu)^2*(1-2*r2)*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2))
		-R*(2*(W-W1*(1-mu))*(1-r2)*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))
			+r2*(1-mu)^2*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2)))
		+(W-W1*(1-mu))*K)/.
{K->(W-W1*(1-mu))*(W-W2)-r2*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))}/.
{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.W->Wb)];

Factor[Det1/Det2]

Bringing all the terms together, we need to multiply this by 
	Ds*(1 - mu)^2*(r2-r1)/(x1*Wb^4)
to get equation 9.

The characteristic polynomial at lambda=1

Lambda1=Ds*(1 - mu)^2*(r2-r1)/(x1*Wb^4)*(
		-R^2*(1-mu)^2*(1-2*r2)*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2))
		-R*(2*(W-W1*(1-mu))*(1-r2)*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))
			+r2*(1-mu)^2*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2)))
		+(W-W1*(1-mu))*K)/.
{K->(W-W1*(1-mu))*(W-W2)-r2*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))}/.
{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.W->Wb;

Value at R=1/2

Factor[Lambda1/.R->1/2]

To get this equation into the form of equation 11, we have to prove that F1 equals F2 where

F1= Factor[(2*Wb^2 - 3*Wb*W1 + 3*mu*Wb*W1 + 2*pm*Wb*W1 - 2*mu*pm*Wb*W1 - 2*mu*Wb*W2 
 - 2*pm*Wb*W2 + 2*mu*pm*Wb*W2 + W1*W2 - mu*W1*W2)/x1]

F2 = Factor[ Wb* ((1-2*W1+W2)*Wb*(1 - pm)+(1-W1)*(Wb-W2))/(1 - pm)]

We can use this equation to solve for the equilibrium level of disequilibrium...

DsSolve = Ds /. Solve[
Factor[(newx1+newx2)-(newx1+newx2+newx3+newx4)*(x1+x2)
	/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
	/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}]==0,	Ds]

Using the above value for D, we note that the following equation will equal zero and we
will solve it for mu...

Factor[(newx1+newx2+newx3+newx4)-Wb
	/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
	/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.Ds->DsSolve]

muSolve = mu/.Solve[%==0,mu][[1]]

Substitute mu in F1 for the above

F1new = Factor[F1/.mu->muSolve]

Now, x1 equals

Factor[newx1/Wb/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.Ds->DsSolve]

and replacing mu in the above equation we get that x1 equals

Factor[%/.mu->muSolve]

Replacing x1 in F1new with the above we get

Factor[F1new/.x1->%]

Factor[%-F2]

which equals F2, concluding the proof that the two ways of writing the characteristic polynomial at R=1/2 are equivalent.

Solution with weak epistasis

When epistasis is  is zero

We develop three equations that equal zero and solve them for the allele frequency (pm=x3+x3), the disequilibrium (Ds) and the mean fitness (Wb).

A1=Factor[(newx1)-Wb*x1/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}
/.{W2->W1^2}]

A2=Factor[(newx1+newx2)-Wb*(x1+x2)/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}
/.{W2->W1^2}]

A3=Factor[(newx1*newx4-newx2*newx3)-Wb^2*(x1*x4-x2*x3)/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}
/.{W2->W1^2}]

Simplify[Solve[{A1==0,A2==0,Ds==0},{pm,Wb,Ds}]]

It is this last factor that is biologically feasible.

When epistasis is  is small

We develop three equations that equal zero.  We write the allele frequency (pm=x3+x3), the disequilibrium (Ds) and the mean fitness (Wb) as
	Wb = (1 - mu)^2+eps*Wbf + O(eps^2)
	Ds = Dsf*eps + O(eps^2)
	pm = (mu*W1)/(1 - mu - W1 + mu*W1)+pmf*eps + O(eps^2)
We then solve the three equations that equal zero for Wb,Dsf, and pmf .

A1=(Factor[(newx1)-Wb*x1/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}
/.{Wb -> (1 - mu)^2+eps*Wbf, Ds -> Dsf*eps, pm -> (mu*W1)/(1 - mu - W1 + mu*W1)+pmf*eps}
/.{W2->W1^2+eps}]/eps)/.eps->0

A2=(Factor[(newx1+newx2)-Wb*(x1+x2)/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}
/.{Wb -> (1 - mu)^2+eps*Wbf, Ds -> Dsf*eps, pm -> (mu*W1)/(1 - mu - W1 + mu*W1)+pmf*eps}
/.{W2->W1^2+eps}]/eps)/.eps->0

A3=(Factor[(newx1*newx4-newx2*newx3)-Wb^2*(x1*x4-x2*x3)/.
{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}
/.{Wb -> (1 - mu)^2+eps*Wbf, Ds -> Dsf*eps, pm -> (mu*W1)/(1 - mu - W1 + mu*W1)+pmf*eps}
/.{W2->W1^2+eps}]/eps)/.eps->0

Sub1=Simplify[Solve[{A1[[3]]==0,A2[[3]]==0,A3==0},{pmf,Wbf,Dsf}]]

OLD VERSION WAS:

Simplify[Solve[{A1[[3]]==0,A2[[3]]==0,A3==0},{pmf,Wbf,Dsf}]]

WbSmall = (1 - mu)^2+Factor[eps*Wbf/.Sub1]
DsSmall = Factor[Dsf*eps/.Sub1]
pmSmall = (mu*W1)/(1 - mu - W1 + mu*W1)+Factor[pmf*eps/.Sub1]

The characteristic polynomial at lambda=1 with weak epistasis

Lambda1

WbSmall=(1 - mu)^2 - (eps*(1 - mu)^2*mu^2*r1)/((1 - W1)^2*((1-mu)^2 - W1^2 + r1*W1^2));
DsSmall=(eps*mu^2*(1 - mu - W1)^2)/((1 - mu)^2*(1 - W1)^4*((1-mu)^2 - W1^2 + r1*W1^2));
pmSmall=(mu*W1)/(1 - mu - W1 + mu*W1) + (eps*mu^2*(1 - mu - W1 + r1*W1))/((1 - mu)*(1 - W1)^3*((1-mu)^2 - W1^2 + r1*W1^2));

Which can be rewritten as

Sm2=(eps*mu^2*(r2 - r1)*(1 - mu - W1)^2*
( (1-mu)^2 - W1^2 (1-r2) + R*W1^2*(2-3*r2) - R^2*W1^2*(1-2*r2) ))/
((1 - mu)^4*(1 - W1)^2*((1-mu)^2 - W1^2 + r1*W1^2))

which equals equation (10)

Factor[Sm1-Sm2]

NUMERICAL ANALYSES OF HAPLOID MODEL

Various factors we will need

The leading eigenvalue

LAM1[mu_,W1_,W2_,r1_,r2_,R_,pm_,Ds_,Wb_]=Ds*(1 - mu)^2*(r2-r1)/(x1*Wb^4)*(
		-R^2*(1-mu)^2*(1-2*r2)*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2))
		-R*(2*(W-W1*(1-mu))*(1-r2)*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))
			+r2*(1-mu)^2*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2)))
		+(W-W1*(1-mu))*K)/.
{K->(W-W1*(1-mu))*(W-W2)-r2*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))}/.
{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.W->Wb;

The disequilibrium at the initial equilibrium

Factor[(newx1+newx2)-(newx1+newx2+newx3+newx4)*(x1+x2)/.{x5->0,x6->0,x7->0,x8->0}
/.x4->1-x1-x2-x3/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}]

Factor[Ds/.Solve[%==0,Ds]]

DsSol=((-1 + pm)*(-pm + 2*mu*pm - mu^2*pm + pm^2 - 2*mu*pm^2 + mu^2*pm^2 + mu*W1 - mu^2*W1 + pm*W1 - 
4*mu*pm*W1 + 3*mu^2*pm*W1 - 2*pm^2*W1 + 4*mu*pm^2*W1 - 2*mu^2*pm^2*W1 + mu^2*W2 + 2*mu*pm*W2 - 
2*mu^2*pm*W2 + pm^2*W2 - 2*mu*pm^2*W2 + mu^2*pm^2*W2))/((-1 + mu)^2*(-1 + r1)*(pm + W1 - 2*pm*W1 - W2 + pm*W2));

DsSolve[mu_,W1_,W2_,r1_,pm_]=((-1 + pm)*(-pm + 2*mu*pm - mu^2*pm + pm^2 - 2*mu*pm^2 + mu^2*pm^2 + mu*W1 - mu^2*W1 + pm*W1 - 
4*mu*pm*W1 + 3*mu^2*pm*W1 - 2*pm^2*W1 + 4*mu*pm^2*W1 - 2*mu^2*pm^2*W1 + mu^2*W2 + 2*mu*pm*W2 - 
2*mu^2*pm*W2 + pm^2*W2 - 2*mu*pm^2*W2 + mu^2*pm^2*W2))/((-1 + mu)^2*(-1 + r1)*(pm + W1 - 2*pm*W1 - W2 + pm*W2));

The allele frequency at the initial equilibrium

Factor[(newx1)-(newx1+newx2+newx3+newx4)*(x1)/.{x5->0,x6->0,x7->0,x8->0}
/.x4->1-x1-x2-x3/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.Ds->DsSol]

The allele frequency at the initial equilibrium

PmSolve[mu_,W1_,W2_,r1_,pm_]=(pm*W1 - 3*mu*pm*W1 + 3*mu^2*pm*W1 - mu^3*pm*W1 - 2*pm^2*W1 + 6*mu*pm^2*W1 - 6*mu^2*pm^2*W1 + 
2*mu^3*pm^2*W1 + pm^3*r1*W1 - 3*mu*pm^3*r1*W1 + 3*mu^2*pm^3*r1*W1 - mu^3*pm^3*r1*W1 - mu*W1^2 + 
2*mu^2*W1^2 - mu^3*W1^2 - pm*W1^2 + 6*mu*pm*W1^2 - 9*mu^2*pm*W1^2 + 4*mu^3*pm*W1^2 + 2*pm^2*W1^2 - 
8*mu*pm^2*W1^2 + 10*mu^2*pm^2*W1^2 - 4*mu^3*pm^2*W1^2 + pm^2*r1*W1^2 - 3*mu*pm^2*r1*W1^2 + 
3*mu^2*pm^2*r1*W1^2 - mu^3*pm^2*r1*W1^2 - 2*pm^3*r1*W1^2 + 6*mu*pm^3*r1*W1^2 - 6*mu^2*pm^3*r1*W1^2 + 
2*mu^3*pm^3*r1*W1^2 - pm*W2 + 4*mu*pm*W2 - 5*mu^2*pm*W2 + 2*mu^3*pm*W2 + 2*pm^2*W2 - 6*mu*pm^2*W2 + 
6*mu^2*pm^2*W2 - 2*mu^3*pm^2*W2 + pm^2*r1*W2 - 4*mu*pm^2*r1*W2 + 5*mu^2*pm^2*r1*W2 - 
2*mu^3*pm^2*r1*W2 - 2*pm^3*r1*W2 + 6*mu*pm^3*r1*W2 - 6*mu^2*pm^3*r1*W2 + 2*mu^3*pm^3*r1*W2 + 
mu*W1*W2 - 4*mu^2*W1*W2 + 3*mu^3*W1*W2 - 6*mu*pm*W1*W2 + 15*mu^2*pm*W1*W2 - 9*mu^3*pm*W1*W2 + 
6*mu*pm^2*W1*W2 - 12*mu^2*pm^2*W1*W2 + 6*mu^3*pm^2*W1*W2 + pm*r1*W1*W2 - 5*mu*pm*r1*W1*W2 + 
7*mu^2*pm*r1*W1*W2 - 3*mu^3*pm*r1*W1*W2 - 6*pm^2*r1*W1*W2 + 21*mu*pm^2*r1*W1*W2 - 
24*mu^2*pm^2*r1*W1*W2 + 9*mu^3*pm^2*r1*W1*W2 + 6*pm^3*r1*W1*W2 - 18*mu*pm^3*r1*W1*W2 + 
18*mu^2*pm^3*r1*W1*W2 - 6*mu^3*pm^3*r1*W1*W2 + mu*W1^2*W2 - mu^2*W1^2*W2 + pm*W1^2*W2 - 
4*mu*pm*W1^2*W2 + 3*mu^2*pm*W1^2*W2 - 2*pm^2*W1^2*W2 + 4*mu*pm^2*W1^2*W2 - 2*mu^2*pm^2*W1^2*W2 - 
mu*r1*W1^2*W2 + 2*mu^2*r1*W1^2*W2 - mu^3*r1*W1^2*W2 - pm*r1*W1^2*W2 + 6*mu*pm*r1*W1^2*W2 - 
9*mu^2*pm*r1*W1^2*W2 + 4*mu^3*pm*r1*W1^2*W2 + 3*pm^2*r1*W1^2*W2 - 11*mu*pm^2*r1*W1^2*W2 + 
13*mu^2*pm^2*r1*W1^2*W2 - 5*mu^3*pm^2*r1*W1^2*W2 - 2*pm^3*r1*W1^2*W2 + 6*mu*pm^3*r1*W1^2*W2 - 
6*mu^2*pm^3*r1*W1^2*W2 + 2*mu^3*pm^3*r1*W1^2*W2 + mu^2*W2^2 - 2*mu^3*W2^2 + pm*W2^2 - 2*mu*pm*W2^2 - 
3*mu^2*pm*W2^2 + 4*mu^3*pm*W2^2 - 2*pm^2*W2^2 + 2*mu*pm^2*W2^2 + 2*mu^2*pm^2*W2^2 - 
2*mu^3*pm^2*W2^2 - pm*r1*W2^2 + 4*mu*pm*r1*W2^2 - 5*mu^2*pm*r1*W2^2 + 2*mu^3*pm*r1*W2^2 + 
3*pm^2*r1*W2^2 - 10*mu*pm^2*r1*W2^2 + 11*mu^2*pm^2*r1*W2^2 - 4*mu^3*pm^2*r1*W2^2 - 2*pm^3*r1*W2^2 + 
6*mu*pm^3*r1*W2^2 - 6*mu^2*pm^3*r1*W2^2 + 2*mu^3*pm^3*r1*W2^2 - mu*W1*W2^2 + 2*mu^2*W1*W2^2 - 
pm*W1*W2^2 + 5*mu*pm*W1*W2^2 - 4*mu^2*pm*W1*W2^2 + 2*pm^2*W1*W2^2 - 4*mu*pm^2*W1*W2^2 + 
2*mu^2*pm^2*W1*W2^2 + mu*r1*W1*W2^2 - 2*mu^2*r1*W1*W2^2 + mu^3*r1*W1*W2^2 + pm*r1*W1*W2^2 - 
5*mu*pm*r1*W1*W2^2 + 7*mu^2*pm*r1*W1*W2^2 - 3*mu^3*pm*r1*W1*W2^2 - 2*pm^2*r1*W1*W2^2 + 
7*mu*pm^2*r1*W1*W2^2 - 8*mu^2*pm^2*r1*W1*W2^2 + 3*mu^3*pm^2*r1*W1*W2^2 + pm^3*r1*W1*W2^2 - 
3*mu*pm^3*r1*W1*W2^2 + 3*mu^2*pm^3*r1*W1*W2^2 - mu^3*pm^3*r1*W1*W2^2);

Factor[(newx1+newx2+newx3+newx4)/.{x5->0,x6->0,x7->0,x8->0}
/.x4->1-x1-x2-x3/.x3->x2
/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.Ds->DsSol]

WbSolve[mu_,W1_,W2_,r1_,pm_]=(W1 - mu*W1 - pm*W1 + mu*pm*W1 - W2 + 2*mu*W2 + 2*pm*W2 - 2*mu*pm*W2 - mu*W1*W2 - pm*W1*W2 +
mu*pm*W1*W2)/(pm + W1 - 2*pm*W1 - W2 + pm*W2);

Factor[(newx1+newx2+newx3+newx4)/.{x5->0,x6->0,x7->0,x8->0}/.W1->1-s/.W2->(1-s)^2+eps
/.x4->1-x1-x2-x3/.x3->x2/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}]

Mean fitness increases if Ds  is negative and eps+s^2 is still positive.
(this assumes that recombination is increased AT the old equilibrium so that recombination has an immediate effect on the

Factor[D[%/.Ds->Dis[r1],r1]]

Mean fitness in the completely epistatic model is 1-mu:

Factor[WbSolve[mu,1,W2,r1,pm]]

Factor[WbSolve[mu,W1,W1^2,r1,pm]]

Factor[WbSolve[mu,1-mu,(1-mu)^2,r1,pm]]

Factor[WbSolve[mu,1,W2,r1,pm]]

A numerical solution for the eigenvalue

Eigensoln is the approximate solution for the leading eigenvalue

Eigen[mu_,W1_,W2_,r1_,r2_,R_,pm_,Ds_,Wb_,lam_]:=
	Det[SMatrix[mu,W1,W2,r1,r2,R,pm,Ds]-Wb*lam*IdentityMatrix[4]]

EigenSoln[mu_,W1_,W2_,r1_,r2_,R_]:=
	Module[{ZZ1,ZZ2,ZZ3,ZZ4},ZZ1=
	If[AA1=Re[pm/.N[Solve[PmSolve[mu,W1,W2,r1,pm]==0,pm],30]];
		AA1[[1]]>0&&AA1[[1]]<0.2,AA1[[1]],
			If[AA1[[2]]>0&&AA1[[2]]<0.2,AA1[[2]],
				If[AA1[[3]]>0&&AA1[[3]]<0.2,AA1[[3]],0]]];
	ZZ2=N[DsSolve[mu,W1,W2,r1,ZZ1],30];
	ZZ3=N[WbSolve[mu,W1,W2,r1,ZZ1],30];
	ZZ4=Max[lam/.Solve[Eigen[mu,W1,W2,r1,r2,R,ZZ1,ZZ2,ZZ3,lam]==0,lam]]]

Solution is the approximate value of the characteristic polynomial at lambda=1

Solution[mu_,W1_,W2_,r1_,r2_,R_]:=
	Module[{ZZ1,ZZ2,ZZ3,ZZ4},ZZ1=
	If[AA1=Re[pm/.N[Solve[PmSolve[mu,W1,W2,r1,pm]==0,pm],30]];
		AA1[[1]]>0&&AA1[[1]]<0.2,AA1[[1]],
			If[AA1[[2]]>0&&AA1[[2]]<0.2,AA1[[2]],
				If[AA1[[3]]>0&&AA1[[3]]<0.2,AA1[[3]],0]]];
	ZZ2=N[DsSolve[mu,W1,W2,r1,ZZ1],30];
	ZZ3=N[WbSolve[mu,W1,W2,r1,ZZ1],30];
	ZZ4=LAM1[mu,W1,W2,r1,r2,R,ZZ1,ZZ2,ZZ3]]

r1=0.1		r2=1/10+/-1/100  mu=0.001

Figure 1.   W1=0.975 W2=0.928

Q1=Plot[EigenSoln[1/1000,975/1000,928/1000,1/10,2/10,R]-1,
{R,0.0,0.5},PlotRange->{Automatic,Automatic}]

Q2=Plot[EigenSoln[1/1000,975/1000,928/1000,1/10,0/10,R]-1,
{R,0.0,0.5},PlotRange->{Automatic,Automatic}]

Q3=Plot[EigenSoln[1/1000,975/1000,928/1000,1/10,11/100,R]-1,
{R,0.0,0.5},PlotRange->{Automatic,Automatic}]

Q4=Plot[EigenSoln[1/1000,975/1000,928/1000,1/10,9/100,R]-1,
{R,0.0,0.5},PlotRange->{Automatic,Automatic}]

Show[Q3,Q4,
Graphics[Text[FontForm["Modifier that decreases recombination",{"Times",12}],{0.28,0.6*10^(-7)}]],
Graphics[Text[FontForm["Modifier that increases recombination",{"Times",12}],{0.28,-0.6*10^(-7)}]],
PlotRange->{{0,0.5},{-1.2*10^(-7),1.2*10^(-7)}},
Frame->True,
FrameLabel->{FontForm["R",{"Times",12}],
	FontForm["l-1",{"Symbol",12}]}]

Figure 2.   W1=0.794 W2=0554

E3=Plot[EigenSoln[1/1000,794/1000,554/1000,1/10,11/100,R]-1,
{R,0.0,0.5},PlotRange->{Automatic,Automatic}]

E4=Plot[EigenSoln[1/1000,794/1000,554/1000,1/10,9/100,R]-1,
{R,0.0,0.5},PlotRange->{Automatic,Automatic}]

Show[E3,E4,
Graphics[Text[FontForm["Modifier that increases recombination",{"Times",12}],{0.28,0.2*10^(-7)}]],
Graphics[Text[FontForm["Modifier that decreases recombination",{"Times",12}],{0.28,-0.2*10^(-7)}]],
PlotRange->{{0,0.5},{-4*10^(-8),4*10^(-8)}},
Frame->True,
FrameLabel->{FontForm["R",{"Times",12}],
	FontForm["l-1",{"Symbol",12}]}]

r1=0.2		r2=0.3 		mu=0.001

Figure 3.  W1=9/10

P1=Plot[EigenSoln[1/1000,9/10,W2,2/10,3/10,0],
{W2,0.01,0.9},PlotRange->{Automatic,Automatic}]

P2=Plot[EigenSoln[1/1000,9/10,W2,2/10,3/10,0.01]-1,
{W2,0.01,0.9},PlotRange->{Automatic,Automatic}]

P3=Plot[EigenSoln[1/1000,9/10,W2,2/10,3/10,0.1]-1,
{W2,0.01,0.9},PlotRange->{Automatic,Automatic}]

P4=Plot[EigenSoln[1/1000,9/10,W2,2/10,3/10,0.5]-1,
{W2,0.01,0.9},PlotRange->{Automatic,Automatic}]

Plot[EigenSoln[1/1000,9/10,W2,2/10,3/10,0.1]-1,
{W2,0.65,0.652},PlotRange->{Automatic,Automatic}]

Plot[EigenSoln[1/1000,9/10,W2,2/10,3/10,0.5]-1,
{W2,0.77,0.81},PlotRange->{Automatic,Automatic}]

Show[P1,P2,P3,P4,
Graphics[Text[FontForm["R=0",{"Times",12}],{0.5,1.7*10^(-6)}]],
Graphics[Text[FontForm["R=0.01",{"Times",12}],{0.2,1*10^(-6)}]],
Graphics[Text[FontForm["R=0.1",{"Times",12}],{0.2,-1.6*10^(-6)}]],
Graphics[Text[FontForm["R=0.5",{"Times",12}],{0.5,-2*10^(-6)}]],
PlotRange->{{0,1.01},{-4*10^(-6),2.2*10^(-6)}},
Frame->True,
FrameLabel->{FontForm["W2",{"Times",12}],
	FontForm["l-1",{"Symbol",12}]}]

Figure 4.  W1=99/100

P1=Plot[EigenSoln[1/1000,99/100,W2,2/10,3/10,0]-1,
{W2,0.9,0.9899},PlotRange->{Automatic,Automatic}]

P4=Plot[EigenSoln[1/1000,99/100,W2,2/10,3/10,0.01]-1,
{W2,0.9,0.9899},PlotRange->{{0.9,1},Automatic}]

P2=Plot[EigenSoln[1/1000,99/100,W2,2/10,3/10,0.1]-1,
{W2,0.9,0.9899},PlotRange->{Automatic,Automatic}]

P3=Plot[EigenSoln[1/1000,99/100,W2,2/10,3/10,0.5]-1,
{W2,0.9,0.9899},PlotRange->{Automatic,Automatic}]

Plot[EigenSoln[1/1000,99/100,W2,2/10,3/10,0.1]-1,
{W2,0.9788,0.979},PlotRange->{Automatic,Automatic}]

Plot[EigenSoln[1/1000,99/100,W2,2/10,3/10,0.5]-1,
{W2,0.9797,0.9801},PlotRange->{Automatic,Automatic}]

Show[P1,P2,P3,P4,
Graphics[Text[FontForm["R=0",{"Times",12}],{0.91,3.7*10^(-6)}]],
Graphics[Text[FontForm["R=0.01",{"Times",12}],{0.91,-6.2*10^(-6)}]],
Graphics[Text[FontForm["R=0.1",{"Times",12}],{0.92,-8*10^(-6)}]],
Graphics[Text[FontForm["R=0.5",{"Times",12}],{0.945,-8*10^(-6)}]],
PlotRange->{{0.9,0.99},{-1*10^(-5),5*10^(-6)}},
Frame->True,
FrameLabel->{FontForm["W2",{"Times",12}],
	FontForm["l-1",{"Symbol",12}]}]

Show[P1,P2,P3,P4,
Graphics[Text[FontForm["R=0",{"Times",12}],{0.955,2*10^(-6)}]],
Graphics[Text[FontForm["R=0.01",{"Times",12}],{0.955,-0.7*10^(-6)}]],
Graphics[Text[FontForm["R=0.1",{"Times",12}],{0.955,-2.7*10^(-6)}]],
Graphics[Text[FontForm["R=0.5",{"Times",12}],{0.955,-4.0*10^(-6)}]],
PlotRange->{{0.95,0.981},{-4.5*10^(-6),2.5*10^(-6)}},
Frame->True,
FrameLabel->{FontForm["W2",{"Times",12}],
	FontForm["l-1",{"Symbol",12}]}]

HAPLOID MODEL - 
Recursions as in text (same answers since the order of rec and mut doesn't matter)

The recursions

sx1 := (x1-r1*(x1*x4-x2*x3)-(r2+R-2*r2*R)*(x1*x6-x2*x5)-
	R*(x1*x7-x3*x5)-(R-r2*R)*(x1*x8-x4*x5)+
	r2*R*(x3*x6-x2*x7)+r2*(x2*x7-x1*x8))
sx2 := (x2+r1*(x1*x4-x2*x3)+(r2+R-2*r2*R)*(x1*x6-x2*x5)-
	R*(x2*x8-x4*x6)+(R-r2*R)*(x3*x6-x2*x7)-
	r2*R*(x1*x8-x4*x5)+r2*(x1*x8-x2*x7))
sx3 := (x3+r1*(x1*x4-x2*x3)-(r2+R-2*r2*R)*(x3*x8-x4*x7)+
	R*(x1*x7-x3*x5)-(R-r2*R)*(x3*x6-x2*x7)+
	r2*R*(x1*x8-x4*x5)-r2*(x3*x6-x4*x5))
sx4 := (x4-r1*(x1*x4-x2*x3)+(r2+R-2*r2*R)*(x3*x8-x4*x7)+
	R*(x2*x8-x4*x6)+(R-r2*R)*(x1*x8-x4*x5)-
	r2*R*(x3*x6-x2*x7)-r2*(x4*x5-x3*x6))
sx5 := (x5-r3*(x5*x8-x6*x7)+(r2+R-2*r2*R)*(x1*x6-x2*x5)+
	R*(x1*x7-x3*x5)+(R-r2*R)*(x1*x8-x4*x5)-
	r2*R*(x3*x6-x2*x7)-r2*(x4*x5-x3*x6))
sx6 := (x6+r3*(x5*x8-x6*x7)-(r2+R-2*r2*R)*(x1*x6-x2*x5)+
	R*(x2*x8-x4*x6)-(R-r2*R)*(x3*x6-x2*x7)+r2*R*(x1*x8-x4*x5)-
	r2*(x3*x6-x4*x5))
sx7 := (x7+r3*(x5*x8-x6*x7)+(r2+R-2*r2*R)*(x3*x8-x4*x7)-
	R*(x1*x7-x3*x5)+(R-r2*R)*(x3*x6-x2*x7)-r2*R*(x1*x8-x4*x5)+
	r2*(x1*x8-x2*x7))
sx8 := (x8-r3*(x5*x8-x6*x7)-(r2+R-2*r2*R)*(x3*x8-x4*x7)-
	R*(x2*x8-x4*x6)-(R-r2*R)*(x1*x8-x4*x5)+r2*R*(x3*x6-x2*x7)+
	r2*(x2*x7-x1*x8))

newx1 := (1-mu)^2*sx1
newx2 := W1*(mu*(1-mu)*sx1+(1-mu)*sx2)
newx3 := W1*(mu*(1-mu)*sx1+(1-mu)*sx3)
newx4 := W2*(mu^2*sx1+mu*sx2+mu*sx3+sx4)
newx5 := (1-mu)^2*sx5
newx6 := W1*(mu*(1-mu)*sx5+(1-mu)*sx6)
newx7 := W1*(mu*(1-mu)*sx5+(1-mu)*sx7)
newx8 := W2*(mu^2*sx5+mu*sx6+mu*sx7+sx8)

Fitness:=(newx1+newx2+newx3+newx4+newx5+newx6+newx7+newx8);

Getting the  stability matrix

F1=Factor[Factor[Series[newx5/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row1={Factor[Coefficient[Expand[F1],x5]],Factor[Coefficient[Expand[F1],x6]],
Factor[Coefficient[Expand[F1],x7]],Factor[Coefficient[Expand[F1],x8]]}

F2=Factor[Factor[Series[newx6/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row2={Factor[Coefficient[Expand[F2],x5]],Factor[Coefficient[Expand[F2],x6]],
Factor[Coefficient[Expand[F2],x7]],Factor[Coefficient[Expand[F2],x8]]}

F3=Factor[Factor[Series[newx7/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row3={Factor[Coefficient[Expand[F3],x5]],Factor[Coefficient[Expand[F3],x6]],
Factor[Coefficient[Expand[F3],x7]],Factor[Coefficient[Expand[F3],x8]]}

F4=Factor[Factor[Series[newx8/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}]

Row4={Factor[Coefficient[Expand[F4],x5]],Factor[Coefficient[Expand[F4],x6]],
Factor[Coefficient[Expand[F4],x7]],Factor[Coefficient[Expand[F4],x8]]}

The  matrix

StabMatrix=Factor[{Row1,Row2,Row3,Row4}/.{x4->1-x1-x2-x3}/.x3->x2];

SMatrix[mu_,W1_,W2_,r1_,r2_,R_,pm_,Ds_]=Factor[{Row1,Row2,Row3,Row4}/.{x4->1-x1-x2-x3}/.x3->x2]/.
	{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)};

Series of matrix manipulations

Here we will manipulate the above stability matrix to get it into a simpler format.

The value of the characteristic polynomial at lambda=1 is equal to:

m = Wb*IdentityMatrix[4] - StabMatrix;

where we have multiplied the characteristic polynomial by the mean fitness^4 (Wb^4).

tm = Transpose[m];

Multiply columns by x_i (need to divide final result by x1*x2*x3*x4)

m2 = Transpose[{tm[[1]]*x1,tm[[2]]*x2,tm[[3]]*x3,tm[[4]]*x4}];

tm2=Transpose[m2];

Add all columns to the first column

m3 = Transpose[{tm2[[1]]+tm2[[2]]+tm2[[3]]+tm2[[4]],
			tm2[[2]],tm2[[3]],tm2[[4]]}]/.x4->1-x2-x3-x1/.x3->x2/.
		{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)};

Make subsutitutions for Wb=newxi/xi in the first column

m4 = Factor[{{m3[[1]][[1]]/.Wb->Factor[(newx1/x1)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[1]][[2]],m3[[1]][[3]],m3[[1]][[4]]},
	{m3[[2]][[1]]/.Wb->Factor[(newx2/x2)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[2]][[2]],m3[[2]][[3]],m3[[2]][[4]]},
	{m3[[3]][[1]]/.Wb->Factor[(newx3/x3)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[3]][[2]],m3[[3]][[3]],m3[[3]][[4]]},
	{m3[[4]][[1]]/.Wb->Factor[(newx4/x4)
			/.{x5->0,x6->0,x7->0,x8->0}/.x4->1-x2-x3-x1/.x3->x2
			/.{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}],m3[[4]][[2]],m3[[4]][[3]],m3[[4]][[4]]}}]

Transpose[m4];

Add row 1 to rows 2 and 3 and subtract it from row 4.

m5 = Factor[{m4[[1]],W1*m4[[1]]+m4[[2]],W1*m4[[1]]+m4[[3]],-W2*m4[[1]]+m4[[4]]}]

The stability matrix now equals a constant (C1) times a 3x3 matrix

C1=m5[[1,1]]

m6 = m5[[{2,3,4},{2,3,4}]]

tm6=Transpose[m6];

Divide the first column by x2 ,  the second column by x3, and the last column by x4.
(Need to multiply the final result by x2*x3*x4)

m7=Factor[Transpose[{tm6[[1]]/(Ds-pm+pm^2),tm6[[2]]/(Ds-pm+pm^2),tm6[[3]]/(Ds+pm^2)}]]

tm7=Transpose[m7];

Subtract column 1 from column 2.  Add column 1 to column 3.

m8=Factor[Transpose[{tm7[[1]],tm7[[2]]-tm7[[1]],tm7[[3]]+tm7[[1]]}]]

Add Row1*W2 +Row2*W2 +Row3*W1.
(Will need to divide the final result by W2)

m9={Factor[m8[[1]]*W2+m8[[2]]*W2+m8[[3]]*W1],m8[[2]],m8[[3]]}

Add Row2/W1 to Row3/W2.
(Will need to multiply the final result by W2)

m10={m9[[1]],m9[[2]],Factor[m9[[3]]/W2+m9[[2]]/W1]}

tm10=Transpose[m10];

Add the first colum to the second column.  Substract the third column from the first.

mFinal=Transpose[{Factor[tm10[[1]]-tm10[[3]]],Factor[tm10[[1]]+tm10[[2]]],tm10[[3]]}]

The determinant of this equals

Det1=Factor[Det[mFinal]];

which can be written as

Det2=Factor[((
		-R^2*(1-mu)^2*(1-2*r2)*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2))
		-R*(2*(W-W1*(1-mu))*(1-r2)*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))
			+r2*(1-mu)^2*(W1^2*(W2-W)+2*W*W1*(W1-W2)*(x1+x2)))
		+(W-W1*(1-mu))*K)/.
{K->(W-W1*(1-mu))*(W-W2)-r2*(1-mu)*(W1*(W2-W)+W*(W1-W2)*(x1+x2))}/.
{x1->Ds+(1-pm)^2,x2->-Ds+pm*(1-pm)}/.W->Wb)];

Factor[Det1/Det2]

Bringing all the terms together, we need to multiply this by 
	Ds*(1 - mu)^2*(r2-r1)/(x1*Wb^4)
to get equation 9.

DIPLOID MODEL

The recursions

xstar1 :=(x1*(x1*w11+x2*w12+x3*w13+x4*w14+x5*w11+x6*w12+x7*w13+x8*w14)-
	r1*(x1*x4*w14-x2*x3*w23)-(r2+R-2*r2*R)*(x1*x6*w12-x2*x5*w12)-R*(
	x1*x7*w13-x3*x5*w13)-(R-r2*R)*(x1*x8*w14-x4*x5*w14)+
	r2*R*(x3*x6*w23-x2*x7*w23)+r2*(x2*x7*w23-x1*x8*w14));

xstar2 := (x2*(x1*w12+x2*w22+x3*w23+x4*w24+x5*w12+x6*w22+x7*w23+x8*w24)+
	r1*(x1*x4*w14-x2*x3*w23)+(r2+R-2*r2*R)*(x1*x6*w12-x2*x5*w12)-R*(
	x2*x8*w24-x4*x6*w24)+(R-r2*R)*(x3*x6*w23-x2*x7*w23)-
	r2*R*(x1*x8*w14-x4*x5*w14)+r2*(x1*x8*w14-x2*x7*w23));

xstar3 := (x3*(x1*w13+x2*w23+x3*w33+x4*w34+x5*w13+x6*w23+x7*w33+x8*w34)+
	r1*(x1*x4*w14-x2*x3*w23)-(r2+R-2*r2*R)*(x3*x8*w34-x4*x7*w34)+R*(
	x1*x7*w13-x3*x5*w13)-(R-r2*R)*(x3*x6*w23-x2*x7*w23)+
	r2*R*(x1*x8*w14-x4*x5*w14)-r2*(x3*x6*w23-x4*x5*w14));

xstar4 := (x4*(x1*w14+x2*w24+x3*w34+x4*w44+x5*w14+x6*w24+x7*w34+x8*w44)-
	r1*(x1*x4*w14-x2*x3*w23)+(r2+R-2*r2*R)*(x3*x8*w34-x4*x7*w34)+R*(
	x2*x8*w24-x4*x6*w24)+(R-r2*R)*(x1*x8*w14-x4*x5*w14)-
	r2*R*(x3*x6*w23-x2*x7*w23)-r2*(x4*x5*w14-x3*x6*w23));

xstar5 := (x5*(x1*w11+x2*w12+x3*w13+x4*w14+x5*w11+x6*w12+x7*w13+x8*w14)-
	r3*(x5*x8*w14-x6*x7*w23)+(r2+R-2*r2*R)*(x1*x6*w12-x2*x5*w12)+R*(
	x1*x7*w13-x3*x5*w13)+(R-r2*R)*(x1*x8*w14-x4*x5*w14)-
	r2*R*(x3*x6*w23-x2*x7*w23)-r2*(x4*x5*w14-x3*x6*w23));

xstar6 := (x6*(x1*w12+x2*w22+x3*w23+x4*w24+x5*w12+x6*w22+x7*w23+x8*w24)+
	r3*(x5*x8*w14-x6*x7*w23)-(r2+R-2*r2*R)*(x1*x6*w12-x2*x5*w12)+R*(
	x2*x8*w24-x4*x6*w24)-(R-r2*R)*(x3*x6*w23-x2*x7*w23)+
	r2*R*(x1*x8*w14-x4*x5*w14)-r2*(x3*x6*w23-x4*x5*w14));

xstar7 := (x7*(x1*w13+x2*w23+x3*w33+x4*w34+x5*w13+x6*w23+x7*w33+x8*w34)+
	r3*(x5*x8*w14-x6*x7*w23)+(r2+R-2*r2*R)*(x3*x8*w34-x4*x7*w34)-R*(
	x1*x7*w13-x3*x5*w13)+(R-r2*R)*(x3*x6*w23-x2*x7*w23)-
	r2*R*(x1*x8*w14-x4*x5*w14)+r2*(x1*x8*w14-x2*x7*w23));

xstar8 := (x8*(x1*w14+x2*w24+x3*w34+x4*w44+x5*w14+x6*w24+x7*w34+x8*w44)-
	r3*(x5*x8*w14-x6*x7*w23)-(r2+R-2*r2*R)*(x3*x8*w34-x4*x7*w34)-R*(
	x2*x8*w24-x4*x6*w24)-(R-r2*R)*(x1*x8*w14-x4*x5*w14)+
	r2*R*(x3*x6*w23-x2*x7*w23)+r2*(x2*x7*w23-x1*x8*w14));

newx1 := (1-mu)^2*xstar1
newx2 := (mu*(1-mu)*xstar1+(1-mu)*xstar2)
newx3 := (mu*(1-mu)*xstar1+(1-mu)*xstar3)
newx4 := (mu^2*xstar1+mu*xstar2+mu*xstar3+xstar4)
newx5 := (1-mu)^2*xstar5
newx6 := (mu*(1-mu)*xstar5+(1-mu)*xstar6)
newx7 := (mu*(1-mu)*xstar5+(1-mu)*xstar7)
newx8 := (mu^2*xstar5+mu*xstar6+mu*xstar7+xstar8)

Fitness=Factor[(newx1+newx2+newx3+newx4+newx5+newx6+newx7+newx8)]
startFit=Factor[newx1+newx2+newx3+newx4/.{x5->0,x6->0,x7->0,x8->0}]

Getting the  stability matrix

F1=Factor[Series[newx5/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}];

Row1={Factor[Coefficient[Expand[F1],x5]],Factor[Coefficient[Expand[F1],x6]],
Factor[Coefficient[Expand[F1],x7]],Factor[Coefficient[Expand[F1],x8]]}

F2=Factor[Series[newx6/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}];

Row2={Factor[Coefficient[Expand[F2],x5]],Factor[Coefficient[Expand[F2],x6]],
Factor[Coefficient[Expand[F2],x7]],Factor[Coefficient[Expand[F2],x8]]}

F3=Factor[Series[newx7/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}];

Row3={Factor[Coefficient[Expand[F3],x5]],Factor[Coefficient[Expand[F3],x6]],
Factor[Coefficient[Expand[F3],x7]],Factor[Coefficient[Expand[F3],x8]]}

F4=Factor[Series[newx8/.
{x5->ee*ex5,x6->ee*ex6,x7->ee*ex7,x8->ee*ex8},{ee,0,1}]/.
{ex5->x5/ee,ex6->x6/ee,ex7->x7/ee,ex8->x8/ee}];

Row4={Factor[Coefficient[Expand[F4],x5]],Factor[Coefficient[Expand[F4],x6]],
Factor[Coefficient[Expand[F4],x7]],Factor[Coefficient[Expand[F4],x8]]}

The  matrix

Row1={(-1 + mu)^2*(w11*x1 + w12*x2 - R*w12*x2 - r2*w12*x2 + 2*R*r2*w12*x2 + w13*x3 - R*w13*x3 + w14*x4 - R*w14*x4 - 
	r2*w14*x4 + R*r2*w14*x4),
	(-1 + mu)^2*(R*w12*x1 + r2*w12*x1 - 2*R*r2*w12*x1 + r2*w23*x3 - R*r2*w23*x3), 
	(-1 + mu)^2*R*(w13*x1 + r2*w23*x2),
	(-1 + mu)^2*R*(1 - r2)*w14*x1};

Row2={(-1 + mu)*(-(mu*w11*x1) - mu*w12*x2 - R*w12*x2 + mu*R*w12*x2 - r2*w12*x2 + mu*r2*w12*x2 + 2*R*r2*w12*x2 - 
	2*mu*R*r2*w12*x2 - mu*w13*x3 + mu*R*w13*x3 - mu*w14*x4 + mu*R*w14*x4 - r2*w14*x4 + mu*r2*w14*x4 + 
	R*r2*w14*x4 - mu*R*r2*w14*x4),
	(-1 + mu)*(-(w12*x1) + R*w12*x1 - mu*R*w12*x1 + r2*w12*x1 - mu*r2*w12*x1 - 
	2*R*r2*w12*x1 + 2*mu*R*r2*w12*x1 - w22*x2 - w23*x3 + R*w23*x3 + r2*w23*x3 - mu*r2*w23*x3 - R*r2*w23*x3 + 
	mu*R*r2*w23*x3 - w24*x4 + R*w24*x4),
	(-1 + mu)*R*(-(mu*w13*x1) - w23*x2 + r2*w23*x2 - mu*r2*w23*x2), 
	(-1 + mu)*R*(-(mu*w14*x1) - r2*w14*x1 + mu*r2*w14*x1 - w24*x2)};

Row3={(-1 + mu)*(-(mu*w11*x1) - mu*w12*x2 + mu*R*w12*x2 + mu*r2*w12*x2 - 2*mu*R*r2*w12*x2 - mu*w13*x3 - R*w13*x3 + 
	mu*R*w13*x3 - mu*w14*x4 + mu*R*w14*x4 + mu*r2*w14*x4 - R*r2*w14*x4 - mu*R*r2*w14*x4), 
	(-1 + mu)*(-(mu*R*w12*x1) - mu*r2*w12*x1 + 2*mu*R*r2*w12*x1 - R*w23*x3 - mu*r2*w23*x3 + R*r2*w23*x3 + 
	mu*R*r2*w23*x3),
	(-1 + mu)*(-(w13*x1) + R*w13*x1 - mu*R*w13*x1 - w23*x2 + R*w23*x2 + r2*w23*x2 - R*r2*w23*x2 - 
	mu*R*r2*w23*x2 - w33*x3 - w34*x4 + R*w34*x4 + r2*w34*x4 - 2*R*r2*w34*x4), 
	(-1 + mu)*(-(mu*R*w14*x1) - r2*w14*x1 + R*r2*w14*x1 + mu*R*r2*w14*x1 - R*w34*x3 - r2*w34*x3 + 2*R*r2*w34*x3)};

Row4={mu^2*w11*x1 + mu^2*w12*x2 + mu*R*w12*x2 - mu^2*R*w12*x2 + mu*r2*w12*x2 - mu^2*r2*w12*x2 - 2*mu*R*r2*w12*x2 + 
	2*mu^2*R*r2*w12*x2 + mu^2*w13*x3 + mu*R*w13*x3 - mu^2*R*w13*x3 + mu^2*w14*x4 + R*w14*x4 - mu^2*R*w14*x4 + 
	mu*r2*w14*x4 - mu^2*r2*w14*x4 - R*r2*w14*x4 + mu^2*R*r2*w14*x4, 
	mu*w12*x1 - mu*R*w12*x1 + mu^2*R*w12*x1 - mu*r2*w12*x1 + mu^2*r2*w12*x1 + 2*mu*R*r2*w12*x1 - 2*mu^2*R*r2*w12*x1 + 
	mu*w22*x2 + mu*w23*x3 - mu*r2*w23*x3 + mu^2*r2*w23*x3 + R*r2*w23*x3 - mu^2*R*r2*w23*x3 + mu*w24*x4 + R*w24*x4 - 
	mu*R*w24*x4, mu*w13*x1 - mu*R*w13*x1 + mu^2*R*w13*x1 + mu*w23*x2 + r2*w23*x2 - mu*r2*w23*x2 - R*r2*w23*x2 + 
	mu^2*R*r2*w23*x2 + mu*w33*x3 + mu*w34*x4 + R*w34*x4 - mu*R*w34*x4 + r2*w34*x4 - mu*r2*w34*x4 - 2*R*r2*w34*x4 + 
	2*mu*R*r2*w34*x4, w14*x1 - R*w14*x1 + mu^2*R*w14*x1 - r2*w14*x1 + mu*r2*w14*x1 + R*r2*w14*x1 - 
	mu^2*R*r2*w14*x1 + w24*x2 - R*w24*x2 + mu*R*w24*x2 + w34*x3 - R*w34*x3 + mu*R*w34*x3 - r2*w34*x3 + 
	mu*r2*w34*x3 + 2*R*r2*w34*x3 - 2*mu*R*r2*w34*x3 + w44*x4};

StabMatrix=Factor[{Row1,Row2,Row3,Row4}];

Series of matrix manipulations

Here we will manipulate the above stability matrix to get it into a simpler format.

The value of the characteristic polynomial at lambda=1 is equal to:

m = Wb*IdentityMatrix[4] - StabMatrix;

where we have multiplied the characteristic polynomial by the mean fitness^4 (Wb^4).

Need to divide by Wb^4

tm = Transpose[m];

Multiply columns by x_i (need to divide final result by x1*x2*x3*x4)

m2 = Transpose[{tm[[1]]*x1,tm[[2]]*x2,tm[[3]]*x3,tm[[4]]*x4}];

tm2=Transpose[m2];

Add all columns to the first column

m3 = Transpose[{tm2[[1]]+tm2[[2]]+tm2[[3]]+tm2[[4]],
			tm2[[2]],tm2[[3]],tm2[[4]]}];

Make subsutitutions for Wb=newxi/xi in the first column

m4 = Factor[{{m3[[1]][[1]]/.Wb->Factor[(newx1/x1)
			/.{x5->0,x6->0,x7->0,x8->0}],m3[[1]][[2]],m3[[1]][[3]],m3[[1]][[4]]},
	{m3[[2]][[1]]/.Wb->Factor[(newx2/x2)
			/.{x5->0,x6->0,x7->0,x8->0}],m3[[2]][[2]],m3[[2]][[3]],m3[[2]][[4]]},
	{m3[[3]][[1]]/.Wb->Factor[(newx3/x3)
			/.{x5->0,x6->0,x7->0,x8->0}],m3[[3]][[2]],m3[[3]][[3]],m3[[3]][[4]]},
	{m3[[4]][[1]]/.Wb->Factor[(newx4/x4)
			/.{x5->0,x6->0,x7->0,x8->0}],m3[[4]][[2]],m3[[4]][[3]],m3[[4]][[4]]}}]

Transpose[m4];

Add row 1 to rows 2 and 3 and subtract it from row 4.

m5 = Factor[{m4[[1]],m4[[1]]+m4[[2]],m4[[1]]+m4[[3]],-m4[[1]]+m4[[4]]}]

The stability matrix now equals a constant (F1) times a 3x3 matrix

F1=m5[[1,1]]

Need to multiply by (-1 + mu)^2*(-r1 + r2)*(-(w23*x2*x3) + w14*x1*x4)

m6 = m5[[{2,3,4},{2,3,4}]];

tm6=Transpose[m6];

Divide the first column by x2 ,  the second column by x3, and the last column by x4.
(Need to multiply the final result by x2*x3*x4)

m7=Factor[Transpose[{tm6[[1]]/x2,tm6[[2]]/x3,tm6[[3]]/x4}]]

Add all rows to row 3.

m8=Factor[{m7[[1]],m7[[2]],m7[[1]]+m7[[2]]+m7[[3]]}]

Rewrite m8 in terms of the marginal allele fitnesses

m9=Factor[m8/.w22*x2->W2dot-w12*x1-w23*x3-w24*x4/.w33*x3->W3dot-w13*x1-w23*x2-w34*x4
/.w11*x1->W1dot-w12*x2-w13*x3-w14*x4/.w44*x4->W4dot-w14*x1-w24*x2-w34*x3]

Make symmetry assumptions

m10=Factor[m9/.{x3->x2,W3dot->W2dot,w13->w12,w23->w22,w33->w22,w34->w24}]

tm10=Transpose[m10];

Subtract column 1 from 2 and column 1 *(Wb-W4dot) to column 3*(Wb-W2dot).
Need to divide final solution by (Wb-W2dot)

m11=Factor[Transpose[{tm10[[1]],tm10[[2]]-tm10[[1]],
	tm10[[3]]*(Wb-W2dot)-tm10[[1]]*(Wb-W4dot)}]]

This is now a 2x2 matrix times (Wb-W2dot).
Need to multiply final solution by (Wb-W2dot)

m12=m11[[{1,2},{2,3}]];

Add first column to second column

m13=Factor[{m12[[1]],m12[[2]]+m12[[1]]}]

The stability matrix at lambda=1 is m13 times
(1-mu)^2*(r2-r1)*(w14*x1*x4-w23*x2*x3)/(Wb^4*x1)

The simplified 2x2 stability matrix

SimpMatrix={{-Wb + W2dot - mu*W2dot - R*w12*x1 + mu*R*w12*x1 - 2*R*w22*x2 + 2*mu*R*w22*x2 - R*w24*x4 + mu*R*w24*x4, 
-Wb^2 + Wb*W2dot - mu*Wb*W2dot + Wb*W4dot - W2dot*W4dot + mu*W2dot*W4dot - R*Wb*w14*x1 + mu*R*Wb*w14*x1 + 
R*w14*W2dot*x1 - mu*R*w14*W2dot*x1 - R*Wb*w22*x2 + mu*R*Wb*w22*x2 - R*Wb*w24*x2 + mu*R*Wb*w24*x2 + 
R*W2dot*w24*x2 - mu*R*W2dot*w24*x2 + R*w22*W4dot*x2 - mu*R*w22*W4dot*x2 - R*Wb*w24*x4 + mu*R*Wb*w24*x4 + 
R*w24*W4dot*x4 - mu*R*w24*W4dot*x4}, {(-1 + mu)*(-1 + 2*R)*r2*(w12*x1 + 2*w22*x2 + w24*x4), 
-Wb^2 + Wb*W2dot - mu*Wb*W2dot + Wb*W4dot - W2dot*W4dot + mu*W2dot*W4dot + R*Wb*w12*x1 - mu*R*Wb*w12*x1 + 
r2*Wb*w12*x1 - mu*r2*Wb*w12*x1 - 2*R*r2*Wb*w12*x1 + 2*mu*R*r2*Wb*w12*x1 - 2*R*Wb*w14*x1 + 2*mu*R*Wb*w14*x1 - 
r2*Wb*w14*x1 + mu*r2*Wb*w14*x1 + 2*R*r2*Wb*w14*x1 - 2*mu*R*r2*Wb*w14*x1 + 2*R*w14*W2dot*x1 - 
2*mu*R*w14*W2dot*x1 + r2*w14*W2dot*x1 - mu*r2*w14*W2dot*x1 - 2*R*r2*w14*W2dot*x1 + 2*mu*R*r2*w14*W2dot*x1 - 
R*w12*W4dot*x1 + mu*R*w12*W4dot*x1 - r2*w12*W4dot*x1 + mu*r2*w12*W4dot*x1 + 2*R*r2*w12*W4dot*x1 - 
2*mu*R*r2*w12*W4dot*x1 + r2*Wb*w22*x2 - mu*r2*Wb*w22*x2 - 2*R*r2*Wb*w22*x2 + 2*mu*R*r2*Wb*w22*x2 - 
2*R*Wb*w24*x2 + 2*mu*R*Wb*w24*x2 - r2*Wb*w24*x2 + mu*r2*Wb*w24*x2 + 2*R*r2*Wb*w24*x2 - 2*mu*R*r2*Wb*w24*x2 + 
2*R*W2dot*w24*x2 - 2*mu*R*W2dot*w24*x2 + r2*W2dot*w24*x2 - mu*r2*W2dot*w24*x2 - 2*R*r2*W2dot*w24*x2 + 
2*mu*R*r2*W2dot*w24*x2 - r2*w22*W4dot*x2 + mu*r2*w22*W4dot*x2 + 2*R*r2*w22*W4dot*x2 - 2*mu*R*r2*w22*W4dot*x2 - 
R*Wb*w24*x4 + mu*R*Wb*w24*x4 + R*w24*W4dot*x4 - mu*R*w24*W4dot*x4}};

The stability matrix at lambda=1 is SimpMatrix times
(1-mu)^2*(r2-r1)*(w14*x1*x4-w23*x2*x3)/(Wb^4*x1)

Working with the determinant of SimpMatrix

Define

L1=(Wb-(1-mu)*W2dot)*(Wb-W4dot)-r2*(1-mu)*((W2dot-Wb)*(x2*w24+x1*w14)+(W4dot-Wb)*(x2*w22+x4*w24))

L2=Factor[((W2dot-Wb)*(x2*w24+x1*w14)+(W4dot-Wb)*(x2*w22+x4*w24))]

L3=Factor[(-2*(Wb-W2dot)*(x1*w14+x2*w24)+(Wb-W4dot)*(w12*x1-x4*w24))
(w12*x1 + 2*w22*x2 + w24*x4)]

R=0 coefficient of SimpMatrix

Part0=Factor[SeriesCoefficient[Series[Det[SimpMatrix],{R,0,0}],0]]

Factor[Part0/L1]

R coefficient of SimpMatrix

Part1=Factor[SeriesCoefficient[Series[Det[SimpMatrix],{R,0,1}],1]]

Factor[Part1+2*(1-mu)*(Wb-W2dot*(1-mu))*(1-r2)*((W2dot-Wb)*(x2*w24+x1*w14)+(W4dot-Wb)*(x2*w22+x4*w24))]

Factor[Part1-(-2*(1-mu)*(Wb-W2dot*(1-mu))*(1-r2)*L2-(1-mu)^2*r2*L3)]

R^2 coefficient of SimpMatrix

Part2=Factor[SeriesCoefficient[Series[Det[SimpMatrix],{R,0,2}],2]]

Factor[Part2/L3]

Overall we have

Lambda1=(1-mu)^2*(r2-r1)*(w14*x1*x4-w23*x2*x3)/(Wb^4*x1)*(
	(Wb - W2dot + mu*W2dot)*L1
	-R*(2*(1-mu)*(Wb-W2dot*(1-mu))*(1-r2)*L2+(1-mu)^2*r2*L3)
	-R^2*(1-mu)^2*(1-2*r2)*L3);

Which gives equation (14)

Analysis at R=1/2

Rhalf=Factor[Lambda1/.R->1/2]

A1=Rhalf[[7]]

A2=Factor[Rhalf/A1/.w22*x2->(W2dot-w12*x1-w24*x4)/2]

A3=Factor[-2*(Wb-(1-mu)*W2dot)*(Wb-W4dot)]

Factor[A1-A3]

Factor[A1-(A3+L3*(1-mu)/W2dot)/.w22*x2->(W2dot-w12*x1-w24*x4)/2]

Overall we have that at R=1/2, the characteristic polynomial at lambda=1 is

(2*(Wb-(1-mu)*W2dot)*(Wb-W4dot)-L3*(1-mu)/W2dot)*
	(1-mu)^2*(r2 - r1)*(2*Wb - W2dot + mu*W2dot)*(w14*x1*x4-w23*x2*x3)/(4*Wb^4*x1);

which is given in the appendix.

If mutation is much weaker than selection ...

Factor[(2*(Wb-(1-mu)*W2dot)*(Wb-W4dot)-L3*(1-mu)/W2dot)/.
	{Wb->1,x2->0,x4->0,x1->1,W2dot->w12,W4dot->w14,mu->0}]

which is the same condition as found in the haploid case.  At R=0 D will determine the outcome.  Same at R=1/2, unless 2-3*W12+W12*W14<0.