Mean fitness in the two-locus haploid selection model

In this Mathematica file, we determine the relationship between mean fitness and epistasis in a two-locus model with mutation. The method is based on a Taylor Series approximation for the mean fitness, keeping terms to order epistasis and epistasis^2.

With several pairs of interacting loci, the epistasis^2 term can be used to predict how variability in the genetic interactions between different pairs of loci will affect mean fitness.

The recursions

Off[General::"spell1"]

(ORDER OF EVENTS: Selection, Recombination, Mutation)

For two selected loci:

Frequencey Genotype Fitness
x1 AB w11 = 1
x2 Ab w12 = 1+s2
x3 aB w21 = 1+s1
x4 ab w22 = (1+s1)(1+s2)+eps

Selection:

jx1 := w11*x1
jx2 := w12*x2
jx3 := w21*x3
jx4 := w22*x4

Fitness:=(jx1+jx2+jx3+jx4);

sx1 := 1/Fitness*jx1
sx2 := 1/Fitness*jx2
sx3 := 1/Fitness*jx3
sx4 := 1/Fitness*jx4

Recombination

nx1 := (sx1-r*(sx1*sx4-sx2*sx3));
nx2 := (sx2+r*(sx1*sx4-sx2*sx3));
nx3 := (sx3+r*(sx1*sx4-sx2*sx3));
nx4 := (sx4-r*(sx1*sx4-sx2*sx3));

Mutation:

newx1 := (1-mu1)*(1-mu2)*nx1
newx2 := (mu2*(1-mu1)*nx1+(1-mu1)*nx2)
newx3 := (mu1*(1-mu2)*nx1+(1-mu2)*nx3)
newx4 := (mu1*mu2*nx1+mu1*nx2+mu2*nx3+nx4)

Factor[newx1+newx2+newx3+newx4]

1

Transformation to allele frequencies and disequilibria from chromosome frequencies:

subtwoloci={x1->(1-p1)*(1-p2)+dis,x2->(1-p1)*p2-dis,x3->p1*(1-p2)-dis,x4->p1*p2+dis};

Fitness definitions for haploid selection model (s will be assumed to be negative):

subfit={w11->1,w12->1+s2,w21->1+s1,w22->(1+s1)(1+s2)+eps};

Two locus analysis

In the following, we will get approximations for the disequilibrium, allele frequency,
and mean fitness to order (eps)^2.

From the following equation (which equals zero at equilibrium), we get an
equation for the disequilibrium at equilibrium:

Factor[(newx1+newx2)-(x1+x2)/.subtwoloci/.subfit]

(-(dis*eps) - mu1 + dis*eps*p1 + mu1*p1 - eps*p1*p2 + eps*p1^2*p2 - p1*s1 + p1^2*s1 - dis*s2 + dis*mu1*s2 -
mu1*p2*s2 + mu1*p1*p2*s2 - dis*s1*s2 + dis*p1*s1*s2 - p1*p2*s1*s2 + p1^2*p2*s1*s2)/
(1 + dis*eps + eps*p1*p2 + p1*s1 + p2*s2 + dis*s1*s2 + p1*p2*s1*s2)

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

{-(((-1 + p1)*(mu1 + eps*p1*p2 + p1*s1 + mu1*p2*s2 + p1*p2*s1*s2))/(-eps + eps*p1 - s2 + mu1*s2 - s1*s2 + p1*s1*s2))}

DsSolve=((1 - p1)*(mu1 + eps*p1*p2 + p1*s1 + mu1*p2*s2 + p1*p2*s1*s2))/(-eps + eps*p1 - s2 + mu1*s2 - s1*s2 + p1*s1*s2);

From the following equations (which equal zero at equilibrium), we get
two equations for the two unknown allele frequencies at equilibrium:

Factor[(newx1+newx3)-(x1+x3)/.subtwoloci/.subfit/.dis->DsSolve]

Factor[newx1-x1/.subtwoloci/.subfit/.dis->DsSolve]

These two equations equal zero. We can drop the denominators to make the calculations simpler:

eqn1[p1_,p2_,eps_]=(eps*mu1 - eps*mu2 - eps*mu1*p1 + eps*mu2*p1 - eps*mu1*p2 + eps*mu2*p2 + eps*mu1*p1*p2 - eps*mu2*p1*p2 + mu1*s1 -
mu1*mu2*s1 + eps*p1*s1 - mu1*p1*s1 - eps*mu2*p1*s1 + mu1*mu2*p1*s1 - eps*p1^2*s1 + eps*mu2*p1^2*s1 + p1*s1^2 -
mu2*p1*s1^2 - p1^2*s1^2 + mu2*p1^2*s1^2 - mu2*s2 + mu1*mu2*s2 - eps*p2*s2 + eps*mu1*p2*s2 + mu2*p2*s2 -
mu1*mu2*p2*s2 + eps*p2^2*s2 - eps*mu1*p2^2*s2 + mu1*s1*s2 - mu2*s1*s2 - mu1*p1*s1*s2 + mu1*mu2*p1*s1*s2 +
mu2*p2*s1*s2 - mu1*mu2*p2*s1*s2 + p1*s1^2*s2 - mu2*p1*s1^2*s2 - p1^2*s1^2*s2 + mu2*p1^2*s1^2*s2 - p2*s2^2 +
mu1*p2*s2^2 + p2^2*s2^2 - mu1*p2^2*s2^2 - p2*s1*s2^2 + mu1*p2*s1*s2^2 + p2^2*s1*s2^2 - mu1*p2^2*s1*s2^2);

eqn2[p1_,p2_,eps_]= (eps^2*mu2 - eps*mu1*mu2 - eps^2*mu1*mu2 + eps*mu1^2*mu2 - 2*eps^2*mu2*p1 + eps*mu1*mu2*p1 + 2*eps^2*mu1*mu2*p1 -
eps*mu1^2*mu2*p1 + eps^2*mu2*p1^2 - eps^2*mu1*mu2*p1^2 - eps^2*mu2*p2 + eps^2*mu1*mu2*p2 + eps^2*mu2*p1*p2 -
eps^2*mu1*mu2*p1*p2 - eps*mu1*r - eps^2*mu1*r + eps*mu1^2*r + eps*mu1*mu2*r + eps^2*mu1*mu2*r - eps*mu1^2*mu2*r +
eps*mu1*p1*r + 2*eps^2*mu1*p1*r - eps*mu1^2*p1*r - eps*mu1*mu2*p1*r - 2*eps^2*mu1*mu2*p1*r + eps*mu1^2*mu2*p1*r -
eps^2*mu1*p1^2*r + eps^2*mu1*mu2*p1^2*r + eps^2*mu1*p2*r - eps^2*mu1*mu2*p2*r - eps^2*p1*p2*r - eps^2*mu1*p1*p2*r +
eps^2*mu2*p1*p2*r + eps^2*mu1*mu2*p1*p2*r + eps^2*p1^2*p2*r - eps^2*mu2*p1^2*p2*r - eps*mu2*p1*s1 +
eps*mu1*mu2*p1*s1 + eps*mu2*p1^2*s1 - eps*mu1*mu2*p1^2*s1 - eps*mu1*r*s1 + eps*mu1*mu2*r*s1 - eps*p1*r*s1 -
eps^2*p1*r*s1 + 3*eps*mu1*p1*r*s1 + eps*mu2*p1*r*s1 + eps^2*mu2*p1*r*s1 - 3*eps*mu1*mu2*p1*r*s1 + eps*p1^2*r*s1 +
2*eps^2*p1^2*r*s1 - 2*eps*mu1*p1^2*r*s1 - eps*mu2*p1^2*r*s1 - 2*eps^2*mu2*p1^2*r*s1 + 2*eps*mu1*mu2*p1^2*r*s1 -
eps^2*p1^3*r*s1 + eps^2*mu2*p1^3*r*s1 - eps*p1*r*s1^2 + eps*mu2*p1*r*s1^2 + 2*eps*p1^2*r*s1^2 -
2*eps*mu2*p1^2*r*s1^2 - eps*p1^3*r*s1^2 + eps*mu2*p1^3*r*s1^2 + eps*mu1*s2 - mu1^2*s2 - eps*mu1^2*s2 + mu1^3*s2 +
2*eps*mu2*s2 - mu1*mu2*s2 - 4*eps*mu1*mu2*s2 + 2*mu1^2*mu2*s2 + 2*eps*mu1^2*mu2*s2 - mu1^3*mu2*s2 - eps*mu1*p1*s2 +
eps*mu1^2*p1*s2 - 2*eps*mu2*p1*s2 + 4*eps*mu1*mu2*p1*s2 - 2*eps*mu1^2*mu2*p1*s2 + eps^2*p2*s2 - 2*eps*mu1*p2*s2 -
eps^2*mu1*p2*s2 + 2*eps*mu1^2*p2*s2 - 2*eps*mu2*p2*s2 + 3*eps*mu1*mu2*p2*s2 - eps*mu1^2*mu2*p2*s2 - eps^2*p1*p2*s2 +
eps^2*mu1*p1*p2*s2 + eps*mu2*p1*p2*s2 - eps*mu1*mu2*p1*p2*s2 - eps^2*p2^2*s2 + eps^2*mu1*p2^2*s2 - mu1*r*s2 -
2*eps*mu1*r*s2 + 2*mu1^2*r*s2 + 2*eps*mu1^2*r*s2 - mu1^3*r*s2 + mu1*mu2*r*s2 + 2*eps*mu1*mu2*r*s2 -
2*mu1^2*mu2*r*s2 - 2*eps*mu1^2*mu2*r*s2 + mu1^3*mu2*r*s2 + 2*eps*mu1*p1*r*s2 - 2*eps*mu1^2*p1*r*s2 -
2*eps*mu1*mu2*p1*r*s2 + 2*eps*mu1^2*mu2*p1*r*s2 + eps*mu1*p2*r*s2 - eps*mu1^2*p2*r*s2 - eps*mu1*mu2*p2*r*s2 +
eps*mu1^2*mu2*p2*r*s2 - eps*p1*p2*r*s2 + eps*mu1*p1*p2*r*s2 + eps*mu2*p1*p2*r*s2 - eps*mu1*mu2*p1*p2*r*s2 +
2*eps*mu2*s1*s2 - mu1*mu2*s1*s2 - 2*eps*mu1*mu2*s1*s2 + mu1^2*mu2*s1*s2 + eps*p1*s1*s2 - 2*mu1*p1*s1*s2 -
eps*mu1*p1*s1*s2 + 2*mu1^2*p1*s1*s2 - mu2*p1*s1*s2 - 4*eps*mu2*p1*s1*s2 + 3*mu1*mu2*p1*s1*s2 +
4*eps*mu1*mu2*p1*s1*s2 - 2*mu1^2*mu2*p1*s1*s2 - eps*p1^2*s1*s2 + eps*mu1*p1^2*s1*s2 + 2*eps*mu2*p1^2*s1*s2 -
2*eps*mu1*mu2*p1^2*s1*s2 - 2*eps*mu2*p2*s1*s2 + 2*eps*mu1*mu2*p2*s1*s2 - 2*eps*p1*p2*s1*s2 + 2*eps*mu1*p1*p2*s1*s2 +
2*eps*mu2*p1*p2*s1*s2 - 2*eps*mu1*mu2*p1*p2*s1*s2 - 2*mu1*r*s1*s2 - 2*eps*mu1*r*s1*s2 + 2*mu1^2*r*s1*s2 +
2*mu1*mu2*r*s1*s2 + 2*eps*mu1*mu2*r*s1*s2 - 2*mu1^2*mu2*r*s1*s2 - p1*r*s1*s2 - 2*eps*p1*r*s1*s2 + 4*mu1*p1*r*s1*s2 +
6*eps*mu1*p1*r*s1*s2 - 3*mu1^2*p1*r*s1*s2 + mu2*p1*r*s1*s2 + 2*eps*mu2*p1*r*s1*s2 - 4*mu1*mu2*p1*r*s1*s2 -
6*eps*mu1*mu2*p1*r*s1*s2 + 3*mu1^2*mu2*p1*r*s1*s2 + 2*eps*p1^2*r*s1*s2 - 4*eps*mu1*p1^2*r*s1*s2 -
2*eps*mu2*p1^2*r*s1*s2 + 4*eps*mu1*mu2*p1^2*r*s1*s2 + eps*mu1*p2*r*s1*s2 - eps*mu1*mu2*p2*r*s1*s2 -
eps*p1*p2*r*s1*s2 - eps*mu1*p1*p2*r*s1*s2 + eps*mu2*p1*p2*r*s1*s2 + eps*mu1*mu2*p1*p2*r*s1*s2 +
eps*p1^2*p2*r*s1*s2 - eps*mu2*p1^2*p2*r*s1*s2 - mu2*p1*s1^2*s2 + mu1*mu2*p1*s1^2*s2 - p1^2*s1^2*s2 +
mu1*p1^2*s1^2*s2 + mu2*p1^2*s1^2*s2 - mu1*mu2*p1^2*s1^2*s2 - mu1*r*s1^2*s2 + mu1*mu2*r*s1^2*s2 - 2*p1*r*s1^2*s2 -
2*eps*p1*r*s1^2*s2 + 4*mu1*p1*r*s1^2*s2 + 2*mu2*p1*r*s1^2*s2 + 2*eps*mu2*p1*r*s1^2*s2 - 4*mu1*mu2*p1*r*s1^2*s2 +
2*p1^2*r*s1^2*s2 + 4*eps*p1^2*r*s1^2*s2 - 3*mu1*p1^2*r*s1^2*s2 - 2*mu2*p1^2*r*s1^2*s2 - 4*eps*mu2*p1^2*r*s1^2*s2 +
3*mu1*mu2*p1^2*r*s1^2*s2 - 2*eps*p1^3*r*s1^2*s2 + 2*eps*mu2*p1^3*r*s1^2*s2 - p1*r*s1^3*s2 + mu2*p1*r*s1^3*s2 +
2*p1^2*r*s1^3*s2 - 2*mu2*p1^2*r*s1^3*s2 - p1^3*r*s1^3*s2 + mu2*p1^3*r*s1^3*s2 + mu1*s2^2 - 2*mu1^2*s2^2 +
mu1^3*s2^2 + mu2*s2^2 - 3*mu1*mu2*s2^2 + 3*mu1^2*mu2*s2^2 - mu1^3*mu2*s2^2 + 2*eps*p2*s2^2 - 2*mu1*p2*s2^2 -
3*eps*mu1*p2*s2^2 + 2*mu1^2*p2*s2^2 + eps*mu1^2*p2*s2^2 - mu2*p2*s2^2 + 2*mu1*mu2*p2*s2^2 - mu1^2*mu2*p2*s2^2 -
eps*p1*p2*s2^2 + eps*mu1*p1*p2*s2^2 - 2*eps*p2^2*s2^2 + 2*eps*mu1*p2^2*s2^2 - mu1*r*s2^2 + 2*mu1^2*r*s2^2 -
mu1^3*r*s2^2 + mu1*mu2*r*s2^2 - 2*mu1^2*mu2*r*s2^2 + mu1^3*mu2*r*s2^2 + mu1*s1*s2^2 - mu1^2*s1*s2^2 +
2*mu2*s1*s2^2 - 4*mu1*mu2*s1*s2^2 + 2*mu1^2*mu2*s1*s2^2 + p1*s1*s2^2 - 3*mu1*p1*s1*s2^2 + 2*mu1^2*p1*s1*s2^2 -
2*mu2*p1*s1*s2^2 + 4*mu1*mu2*p1*s1*s2^2 - 2*mu1^2*mu2*p1*s1*s2^2 + 2*eps*p2*s1*s2^2 - 2*mu1*p2*s1*s2^2 -
2*eps*mu1*p2*s1*s2^2 + 2*mu1^2*p2*s1*s2^2 - 2*mu2*p2*s1*s2^2 + 3*mu1*mu2*p2*s1*s2^2 - mu1^2*mu2*p2*s1*s2^2 -
2*p1*p2*s1*s2^2 - 2*eps*p1*p2*s1*s2^2 + 2*mu1*p1*p2*s1*s2^2 + 2*eps*mu1*p1*p2*s1*s2^2 + mu2*p1*p2*s1*s2^2 -
mu1*mu2*p1*p2*s1*s2^2 - 2*eps*p2^2*s1*s2^2 + 2*eps*mu1*p2^2*s1*s2^2 - 2*mu1*r*s1*s2^2 + 2*mu1^2*r*s1*s2^2 +
2*mu1*mu2*r*s1*s2^2 - 2*mu1^2*mu2*r*s1*s2^2 - p1*r*s1*s2^2 + 4*mu1*p1*r*s1*s2^2 - 3*mu1^2*p1*r*s1*s2^2 +
mu2*p1*r*s1*s2^2 - 4*mu1*mu2*p1*r*s1*s2^2 + 3*mu1^2*mu2*p1*r*s1*s2^2 + mu2*s1^2*s2^2 - mu1*mu2*s1^2*s2^2 +
p1*s1^2*s2^2 - mu1*p1*s1^2*s2^2 - 2*mu2*p1*s1^2*s2^2 + 2*mu1*mu2*p1*s1^2*s2^2 - p1^2*s1^2*s2^2 +
mu1*p1^2*s1^2*s2^2 + mu2*p1^2*s1^2*s2^2 - mu1*mu2*p1^2*s1^2*s2^2 - mu2*p2*s1^2*s2^2 + mu1*mu2*p2*s1^2*s2^2 -
2*p1*p2*s1^2*s2^2 + 2*mu1*p1*p2*s1^2*s2^2 + mu2*p1*p2*s1^2*s2^2 - mu1*mu2*p1*p2*s1^2*s2^2 - mu1*r*s1^2*s2^2 +
mu1*mu2*r*s1^2*s2^2 - 2*p1*r*s1^2*s2^2 + 4*mu1*p1*r*s1^2*s2^2 + 2*mu2*p1*r*s1^2*s2^2 - 4*mu1*mu2*p1*r*s1^2*s2^2 +
2*p1^2*r*s1^2*s2^2 - 3*mu1*p1^2*r*s1^2*s2^2 - 2*mu2*p1^2*r*s1^2*s2^2 + 3*mu1*mu2*p1^2*r*s1^2*s2^2 - p1*r*s1^3*s2^2 +
mu2*p1*r*s1^3*s2^2 + 2*p1^2*r*s1^3*s2^2 - 2*mu2*p1^2*r*s1^3*s2^2 - p1^3*r*s1^3*s2^2 + mu2*p1^3*r*s1^3*s2^2 +
p2*s2^3 - 2*mu1*p2*s2^3 + mu1^2*p2*s2^3 - p2^2*s2^3 + mu1*p2^2*s2^3 + 2*p2*s1*s2^3 - 3*mu1*p2*s1*s2^3 +
mu1^2*p2*s1*s2^3 - p1*p2*s1*s2^3 + mu1*p1*p2*s1*s2^3 - 2*p2^2*s1*s2^3 + 2*mu1*p2^2*s1*s2^3 + p2*s1^2*s2^3 -
mu1*p2*s1^2*s2^3 - p1*p2*s1^2*s2^3 + mu1*p1*p2*s1^2*s2^3 - p2^2*s1^2*s2^3 + mu1*p2^2*s1^2*s2^3);

These are both solved in the absence of epistasis by setting p1=-mu1/s1 and p2=-mu2/s2
[RECALL that selection is assumed to be negative.]:

Factor[DsSolve/.p1->-mu1/s1/.p2->-mu2/s2/.eps->0]

0

Factor[eqn1[-mu1/s1,-mu2/s2,0]]

0

Factor[eqn2[-mu1/s1,-mu2/s2,0]]

0

The mean fitness of the population is:

Factor[Fitness]

w11 x1+w12 x2+w21 x3+w22 x4

Simplify[Fitness/.subfit/.subtwoloci]

eps p1 p2+(1+p1 s1) (1+p2 s2)+dis (eps+s1 s2)

eqn3[p1_,p2_,eps_]=Factor[Fitness/.subtwoloci/.subfit/.dis->DsSolve]

-(((-1 + mu1)*(-eps + eps*p1 - s2 - eps*p2*s2 - s1*s2 - p2*s2^2 - p2*s1*s2^2))/
(-eps + eps*p1 - s2 + mu1*s2 - s1*s2 + p1*s1*s2))

We will do a Taylor Series of this equation for the mean fitness around eps=0,
keeping linear and quadratic terms in eps:

Simplify[Normal[Series[eqn3[p1[eps],p2[eps],eps],{eps,0,2}]]]

((-1 + mu1)*((1 + s1)*(-1 + mu1 + s1*(-1 + p1[0]))*(1 + s2*p2[0]) +
(eps*(-((1 + s1)*s2*(1 + s2*p2[0])*(-1 + p1[0] + s1*s2*p1'[0])) +
s2*(-1 + mu1 + s1*(-1 + p1[0]))*(1 - p1[0] + s2*p2[0] + s2^2*p2'[0] +
s1*s2^2*p2'[0])))/s2^2 +
(eps^2*((-1 + p1[0] + s1*s2*p1'[0])*
(-1 + p1[0] - s2*p2[0] - s2^2*p2'[0] - s1*s2^2*p2'[0]) -
((1 + s1)*(1 + s2*p2[0])*(-2 - 2*p1[0]^2 + 2*(-1 + mu1 + s1)*s2*p1'[0] -
s1*s2^2*(2*s1*p1'[0]^2 + p1''[0] - mu1*p1''[0] +
s1*p1''[0]) + p1[0]*
(4 - 2*s1*s2*p1'[0] + s1^2*s2^2*p1''[0])))/
(2*(-1 + mu1 + s1*(-1 + p1[0]))) +
(s2*(-1 + mu1 + s1*(-1 + p1[0]))*(-2*p1'[0] +
s2*(2*p2'[0] + (1 + s1)*s2*p2''[0])))/2))/s2^2))/
(-1 + mu1 + s1*(-1 + p1[0]))^2

This Taylor Series depends on terms like p1'[0] and p1''[0], the first and second
derivative with respect to eps of the allele frequency functions. These can be
determined by differentiation of the two equations (eqn1 and eqn2) giving p1 and p2 implicitly.

The first derivatives with respect to eps:

a1=Factor[D[eqn1[p1[eps],p2[eps],eps],eps]/.eps->0/.{p1[0]->-mu1/s1,p2[0]->-mu2/s2}]

a2=Factor[D[eqn2[p1[eps],p2[eps],eps],eps]/.eps->0/.{p1[0]->-mu1/s1,p2[0]->-mu2/s2}]

Since eqn1 and eqn2 equal zero, their derivatives give us functions thac can be solved
for p1'[0] and p2'[0]:

Simplify[Solve[{a1==0,a2==0},{p1'[0],p2'[0]}]]

{{p1'[0] ->
(mu1*mu2*(mu1 + s1 - r*(1 + s1)))/(s1^2*s2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)),
p2'[0] ->
(mu1*mu2*(mu2 + s2 - r*(1 + s2)))/(s1*s2^2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2))}}

subepszero={p1[0]->-mu1/s1,p2[0]->-mu2/s2,
p1'[0] ->
(mu1*mu2*(mu1 + s1 - r*(1 + s1)))/(s1^2*s2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)),
p2'[0] ->
(mu1*mu2*(mu2 + s2 - r*(1 + s2)))/(s1*s2^2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2))};

Next we find the double derivatives with respect to eps:

a3=Factor[D[eqn1[p1[eps],p2[eps],eps],eps,eps]/.eps->0/.subepszero]

a4=Factor[D[eqn2[p1[eps],p2[eps],eps],eps,eps]/.eps->0/.subepszero]

These give two equations in two unknowns (p1''[0] and p2''[0]), which can be solved:

a5=Solve[a3==0,p1''[0]];

Factor[a4/.%]

a6=Simplify[Solve[%==0,p2''[0] ]]

{{p2''[0] ->
(-2*mu1*mu2*(mu1^2*(-1 + mu2)*(mu2^2 + r*(r - s2 + r*s2) + mu2*(-3*r + r^2 + s2 - 2*r*s2 + r^2*s2)) +
(mu2 - r + s2 - r*s2)*(mu2^2*r + (-1 + r)*s1*s2*(r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2) -
mu2*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2)) -
mu1*(mu2^3*(1 + r) + mu2^2*(-3*r - 2*r^2 + s1 + r*s1 - 2*r^2*s1 + 2*s2 - 2*r^2*s2 + 2*s1*s2 - 2*r^2*s1*s2) -
(r - s2 + r*s2)*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2) +
mu2*(3*r^2 + r^3 - 3*r*s1 + 2*r^2*s1 + r^3*s1 - 4*r*s2 + 2*r^2*s2 + 2*r^3*s2 - 4*r*s1*s2 + 2*r^2*s1*s2 +
2*r^3*s1*s2 + s2^2 - r*s2^2 - r^2*s2^2 + r^3*s2^2 + 2*s1*s2^2 - 3*r*s1*s2^2 + r^3*s1*s2^2))))/
(s1^2*s2^3*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^3)}}

a7=Solve[a4==0,p2''[0]];

Factor[a3/.%]

a8=Simplify[Solve[%==0,p1''[0] ]]

{{p1''[0] ->
(-2*mu1*mu2*(mu1^3*(-1 + mu2)*(mu2 - r) - (r - s1 + r*s1)*
(mu2^2*r + (-1 + r)*s1*s2*(r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2) -
mu2*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2)) +
mu1^2*(-2*r^2 + 2*r*s1 - 2*r^2*s1 + mu2^2*(-1 - 3*r + r^2 + s1 - 2*r*s1 + r^2*s1) + r*s2 - r^2*s2 + s1*s2 -
r^2*s1*s2 + mu2*(3*r + 2*r^2 - 2*s1 + 2*r^2*s1 - s2 - r*s2 + 2*r^2*s2 - 2*s1*s2 + 2*r^2*s1*s2)) +
mu1*(r^3 - 2*r^2*s1 + 2*r^3*s1 + r*s1^2 - 2*r^2*s1^2 + r^3*s1^2 + mu2^2*(3*r - s1 + r*s1) - r^2*s2 +
r^3*s2 - r*s1*s2 - r^2*s1*s2 + 2*r^3*s1*s2 + 2*s1^2*s2 - 3*r*s1^2*s2 + r^3*s1^2*s2 + s1*s2^2 -
2*r*s1*s2^2 + r^2*s1*s2^2 + s1^2*s2^2 - 2*r*s1^2*s2^2 + r^2*s1^2*s2^2 -
mu2*(3*r^2 + r^3 - 4*r*s1 + 2*r^2*s1 + 2*r^3*s1 + s1^2 - r*s1^2 - r^2*s1^2 + r^3*s1^2 - 3*r*s2 +
2*r^2*s2 + r^3*s2 - 4*r*s1*s2 + 2*r^2*s1*s2 + 2*r^3*s1*s2 + 2*s1^2*s2 - 3*r*s1^2*s2 + r^3*s1^2*s2))))/
(s1^3*s2^2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^3)}}

Clear[a1,a2,a3,a4,a5,a6,a7,a8]

subepszero={p1[0]->-mu1/s1,p2[0]->-mu2/s2,
p1'[0] ->
(mu1*mu2*(mu1 + s1 - r*(1 + s1)))/(s1^2*s2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)),
p2'[0] ->
(mu1*mu2*(mu2 + s2 - r*(1 + s2)))/(s1*s2^2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)),

p1''[0] ->
(-2*mu1*mu2*(mu1^3*(-1 + mu2)*(mu2 - r) - (r - s1 + r*s1)*
(mu2^2*r + (-1 + r)*s1*s2*(r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2) -
mu2*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2)) +
mu1^2*(-2*r^2 + 2*r*s1 - 2*r^2*s1 + mu2^2*(-1 - 3*r + r^2 + s1 - 2*r*s1 + r^2*s1) + r*s2 - r^2*s2 + s1*s2 -
r^2*s1*s2 + mu2*(3*r + 2*r^2 - 2*s1 + 2*r^2*s1 - s2 - r*s2 + 2*r^2*s2 - 2*s1*s2 + 2*r^2*s1*s2)) +
mu1*(r^3 - 2*r^2*s1 + 2*r^3*s1 + r*s1^2 - 2*r^2*s1^2 + r^3*s1^2 + mu2^2*(3*r - s1 + r*s1) - r^2*s2 +
r^3*s2 - r*s1*s2 - r^2*s1*s2 + 2*r^3*s1*s2 + 2*s1^2*s2 - 3*r*s1^2*s2 + r^3*s1^2*s2 + s1*s2^2 -
2*r*s1*s2^2 + r^2*s1*s2^2 + s1^2*s2^2 - 2*r*s1^2*s2^2 + r^2*s1^2*s2^2 -
mu2*(3*r^2 + r^3 - 4*r*s1 + 2*r^2*s1 + 2*r^3*s1 + s1^2 - r*s1^2 - r^2*s1^2 + r^3*s1^2 - 3*r*s2 +
2*r^2*s2 + r^3*s2 - 4*r*s1*s2 + 2*r^2*s1*s2 + 2*r^3*s1*s2 + 2*s1^2*s2 - 3*r*s1^2*s2 + r^3*s1^2*s2))))/
(s1^3*s2^2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^3),

p2''[0] ->
(-2*mu1*mu2*(mu1^2*(-1 + mu2)*(mu2^2 + r*(r - s2 + r*s2) + mu2*(-3*r + r^2 + s2 - 2*r*s2 + r^2*s2)) +
(mu2 - r + s2 - r*s2)*(mu2^2*r + (-1 + r)*s1*s2*(r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2) -
mu2*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2)) -
mu1*(mu2^3*(1 + r) + mu2^2*(-3*r - 2*r^2 + s1 + r*s1 - 2*r^2*s1 + 2*s2 - 2*r^2*s2 + 2*s1*s2 - 2*r^2*s1*s2) -
(r - s2 + r*s2)*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2) +
mu2*(3*r^2 + r^3 - 3*r*s1 + 2*r^2*s1 + r^3*s1 - 4*r*s2 + 2*r^2*s2 + 2*r^3*s2 - 4*r*s1*s2 + 2*r^2*s1*s2 +
2*r^3*s1*s2 + s2^2 - r*s2^2 - r^2*s2^2 + r^3*s2^2 + 2*s1*s2^2 - 3*r*s1*s2^2 + r^3*s1*s2^2))))/
(s1^2*s2^3*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^3)};

To zero order in eps, fitness equals:

Factor[eqn3[p1[eps],p2[eps],eps]/.eps->0/.subepszero]

(-1+mu1) (-1+mu2)

The first order term in eps of the mean fitness equals:

Factor[D[eqn3[p1[eps],p2[eps],eps],eps]/.eps->0/.subepszero]

-(((-1 + mu1)*mu1*(-1 + mu2)*mu2*r)/(s1*s2*(-mu1 - mu2 + mu1*mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)))

The second order term in eps of the mean fitness equals:

Factor[D[eqn3[p1[eps],p2[eps],eps],eps,eps]/.eps->0/.subepszero]

(2*(-1 + mu1)*mu1*(-1 + mu2)*mu2*r*(-(mu1^2*mu2) - mu1*mu2^2 + mu1^2*mu2^2 + mu1^2*r + mu1*mu2*r - mu1^2*mu2*r +
mu2^2*r - mu1*mu2^2*r - mu1*r^2 - mu2*r^2 + 2*mu1*mu2*r^2 - mu1*mu2*s1 + mu1*r*s1 + mu2*r*s1 - mu1*mu2*r*s1 -
mu1*r^2*s1 - mu2*r^2*s1 + 2*mu1*mu2*r^2*s1 - mu1*mu2*s2 + mu1*r*s2 + mu2*r*s2 - mu1*mu2*r*s2 - mu1*r^2*s2 -
mu2*r^2*s2 + 2*mu1*mu2*r^2*s2 + mu1*s1*s2 + mu2*s1*s2 - 2*mu1*mu2*s1*s2 - r*s1*s2 + r^2*s1*s2 - mu1*r^2*s1*s2 -
mu2*r^2*s1*s2 + 2*mu1*mu2*r^2*s1*s2 + s1^2*s2 - 2*r*s1^2*s2 + r^2*s1^2*s2 + s1*s2^2 - 2*r*s1*s2^2 +
r^2*s1*s2^2 + s1^2*s2^2 - 2*r*s1^2*s2^2 + r^2*s1^2*s2^2))/
(s1^2*s2^2*(-mu1 - mu2 + mu1*mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^3)

Factor[(%/(mu1 mu2))/.mu1->0/.mu2->0]

(2*(-1 + r)*r)/(s1*s2*(r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^2)

This is always negative and becomes more negative as the variability in epistasis
increases (contributing to E[eps^2]).
Thus mean fitness is expected to go down as the variability in epistasis goes up.

(Note that you need to keep up to O(eps^2) throughout!
The sign changes if p1 and p2 are solved only to order eps.)

Finally, let's look at the dependence of the average growth rate, or log(mean fitness), on eps:

Simplify[Log[Fitness]/.subfit/.subtwoloci]

Log[eps p1 p2+(1+p1 s1) (1+p2 s2)+dis (eps+s1 s2)]

eqn4[p1_,p2_,eps_]=Simplify[Log[Fitness]/.subtwoloci/.subfit/.dis->DsSolve]

Log[-(((-1 + mu1)*(eps*(-1 + p1 - p2*s2) - (1 + s1)*s2*(1 + p2*s2)))/(eps*(-1 + p1) + (-1 + mu1 - s1 + p1*s1)*s2))]

To zero order in eps, log fitness equals:

Simplify[eqn4[p1[eps],p2[eps],eps]/.eps->0/.subepszero]

Log[(-1+mu1) (-1+mu2)]

The first order term in eps of the log mean fitness equals:

Simplify[D[eqn4[p1[eps],p2[eps],eps],eps]/.eps->0/.subepszero]

-((mu1*mu2*r)/(s1*s2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)))

To leading order in the mutation rate:

(mu1 mu2)*Simplify[(%/(mu1 mu2))/.mu1->0/.mu2->0]

(mu1*mu2*r)/(s1*s2*(s1 + s2 + s1*s2 - r*(1 + s1)*(1 + s2)))

The second order term in eps of the log mean fitness equals:

Simplify[D[eqn4[p1[eps],p2[eps],eps],eps,eps]/.eps->0/.subepszero]

-((mu1*mu2*r*(mu1^2*(-1 + mu2)*(mu2*(-2 + r) + 2*r) +
mu1*(mu2^2*(2 + r) + mu2*(-2*r - 3*r^2 + 2*s1 + r*s1 - 3*r^2*s1 + 2*s2 + r*s2 - 3*r^2*s2 + 4*s1*s2 -
r*s1*s2 - 3*r^2*s1*s2) + 2*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2)) -
2*(mu2^2*r + (-1 + r)*s1*s2*(r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2) -
mu2*(r^2 - r*s1 + r^2*s1 - r*s2 + r^2*s2 - s1*s2 + r^2*s1*s2))))/
(s1^2*s2^2*(mu1*(-1 + mu2) - mu2 + r - s1 + r*s1 - s2 + r*s2 - s1*s2 + r*s1*s2)^3))

To leading order in the mutation rate:

(mu1 mu2)*Simplify[(%/(mu1 mu2))/.mu1->0/.mu2->0]

(2*mu1*mu2*(-1 + r)*r)/(s1*s2*(s1 + s2 + s1*s2 - r*(1 + s1)*(1 + s2))^2)

This is always negative.
Thus log(mean fitness) is expected to go down as variability in epistasis goes up.

The denominator (phi) can also be written as

Simplify[(1-(1-r)*(1+s1)*(1+s2))/
(-s1-s2-s1 s2+r (1+s1) (1+s2))]

1