Chapter 12: The Evolution of Dominance

Inputs (Enter first)

SP[x_] := Simplify[PowerExpand[x]]

Off[General::"spell"]

Off[General::"spell1"]

<< Graphics`Graphics`

Model set up (Enter second)

Example: Evolution of Dominance

	In 1930, Fisher noted that the “pronounced tendency of the mutant gene to be recessive, to the gene of wild type from which it arises, calls for explanation.” Fisher proposed an explanation in which the degree of dominance is an evolved characteristic. He supposed that there is a second locus that modifies the degree of dominance of alleles at the locus of interest and argued that modifier alleles increasing the dominance of the wildtype allele at the first locus would be favoured. Let’s build a mathematical model to explore Fisher’s hypothesis in more detail. In modeling the evolution of dominance, we assume that there are two loci, one of which is at mutation-selection balance, the other of which modifies the degree of dominance of the wild-type allele at the first locus (e.g., a regulatory gene controlling expression levels). A mutant allele is then introduced at the modifier locus, and the conditions for its invasion are determined. 
	Let’s use A to denote the locus under selection, and B to denote the modifier locus. We use A1 to denote the selectively favoured allele at the A locus, and A2 the deleterious mutant. Specifically, we suppose that the fitness of an A1A1 individual is standardized to equal one, while the fitness of an A2A2 individual is 1 – s. The fitness of an A1A2 heterozygote is assumed to depend on the alleles carried at the modifier locus, which determines the level of dominance of the wildtype allele:
	
	Fitness of A1A2 in B1B1 = 1- h11
	Fitness of A1A2 in B1B2 = 1- h12
	Fitness of A1A2 in B2B2 = 1- h22
	
We assume a life cycle that entails random union of gametes, diploid selection, meiosis with recombination at rate r, and finally mutation at rate μ to produce the next generation of gametes.  Starting with  gamete frequencies:

	x1=A1 B1
	x2 = A1 B2
	x3 = A2 B1
	x4 = A2 B2

The frequencies of the gametes in the next generation are:

After random mating, selection, and meiosis:

w[1, 
          1] = w[1, 
            2] = w[2, 1] = w[2, 2] = 
              1;
w[3, 
          3] = w[3, 
            4] = w[4, 3] = w[4, 4] = 
              1 - s;
w[3, 1] = w[1, 3] = 1 -  h11 s;
w[3, 
          2] = w[2, 
            3] = w[1, 4] = w[4, 1] = 
              1  -  h12 s;
w[2, 4] = w[4, 2] = 1 - h22 s;

wbar = 
        Sum[x[i] x[j] w[i, j], {i, 1, 4}, {j, 1, 
            4}] ;
dis = 
        x[1] x[4] (w[1, 4] + w[4, 1])/(2 wbar) - 
          x[2] x[3] (w[2, 3] + w[3, 2])/(2 wbar);

x1p = (1 - μ) (x[1] Sum[
                  x[j] (w[1, j] + w[j, 1])/(2 wbar), {j, 1, 
                    4}]  -  r dis);
x2p = (1 - μ) (x[2] Sum[
                  x[j] (w[2, j] + w[j, 2])/(2 wbar), {j, 1, 
                    4}]  +  r dis);
x3p = μ (x[1] Sum[
                    x[j] (w[1, j] + w[j, 1])/(2 wbar), {j, 1, 
                      4}]  -  
                r dis) + 
(x[3] Sum[
                  x[j] (w[3, j] + w[j, 3])/(2 wbar), {j, 1, 
                    4}]  +  r dis);
x4p = μ (x[2] Sum[
                    x[j] (w[2, j] + w[j, 2])/(2 wbar), {j, 1, 
                      4}]  +  
                r dis) + 
(x[4] Sum[
                  x[j] (w[4, j] + w[j, 4])/(2 wbar), {j, 1, 
                    4}]  -  r dis);

The frequencies correctly sum to one:

Factor[x1p + x2p + x3p + x4p]

The equilibrium conditions results in a quadratic equation:

eqcond = 
      Factor[x3p - x[3] /. x[2] -> 0 /. x[4] -> 0 /. 
          x[1] -> 1 - x[3]]

The stability matrix:

mat = {{D[x1p, x[1]], D[x1p, x[2]], D[x1p, x[3]], 
                D[x1p, x[4]]}, 
{D[x2p, x[1]], 
                D[x2p, x[2]], D[x2p, x[3]], 
                D[x2p, x[4]]}, 
{D[x3p, x[1]], 
                D[x3p, x[2]], D[x3p, x[3]], 
                D[x3p, x[4]]}, 
{D[x4p, x[1]], 
                D[x4p, x[2]], D[x4p, x[3]], D[x4p, x[4]]}} /. {x[2] -> 
                0, x[4] -> 0} // Simplify;
MatrixForm[mat]

The above does not have a block triangular form because the mutant recursions x2p and x4p are separated by x3p.  Reordering these recursions as x1p ,x3p, x2p and x4p gives a block upper triangular matrix:

mat = {{D[x1p, x[1]], D[x1p, x[3]], D[x1p, x[2]], 
                D[x1p, x[4]]}, 
{D[x3p, x[1]], 
                D[x3p, x[3]], D[x3p, x[2]], 
                D[x3p, x[4]]}, 
{D[x2p, x[1]], 
                D[x2p, x[3]], D[x2p, x[2]], 
                D[x2p, x[4]]}, 
{D[x4p, x[1]], 
                D[x4p, x[3]], D[x4p, x[2]], D[x4p, x[4]]}} /. {x[2] -> 
                0, x[4] -> 0} // Simplify;
MatrixForm[mat]

with the mutant submatrix:

submat = {{D[x2p, x[2]], 
                D[x2p, x[4]]}, 
{D[x4p, x[2]], 
                D[x4p, x[4]]}} /. {x[2] -> 0, x[4] -> 0} // 
          Simplify;
MatrixForm[submat]

If the leading eigenvalue of this matrix, λ > 1, then the new modifier allele B2 spreads through the population.  We use three different approaches to show that modifier alleles can only spread that decrease h, thereby increasing the level of dominance of the wildtype allele.

In the following, let's figure out how sperm and accessory protein allocation should evolve through the successive appearance and fixation of mutations.

Approach 1:  Directly calculating the eigenvalue

The equilibrium condition is a quadratic, which  can be solved (although the equilibrium is cumbersome to work with).  Make sure to have entered the inputs in the first section above, because we will use SP[x_] = Simplify[PowerExpand[x]] repeatedly to expand and evaluate the terms.

SP[Solve[eqcond == 0, x[3]]]

Let's set mu = 0 to determine which eigenvalue is relevant:

SP[% /. μ -> 0]

Because the equilibrium frequency of the deleterious allele A2 should be near zero, we choose the solution near zero:

subEQ = {x[
              3] -> ((h11 s (1 + μ) - (Sqrt[s (4 μ - 8 h11 μ + (h11^2) s (((1 + μ))^2))]))/(2 (-1 + 2 h11) s))};

We can then plug this into the mutant submatrix and solve for its eigenvalues:

SP[submat /. x[1] -> 1 - x[3] /. subEQ]

Calculating the eigenvalues takes a while...and gives a long and complicated answer:

SP[Eigenvalues[%]]

The above is uninterpretable.  Let's take a Taylor Series of it with respect to μ, assuming that the mutation rate is small.

SP[Normal[Series[%, {μ, 0, 1}]]]

Mathematica often can simplify things further using Simplify[PowerExpand[x]] (sometimes it can help to enter the following command multiple times in a row):

eigensoln = SP[%]

Setting h12 = h11 indicates that the second solution is the relevant one (at least with our version of Mathematica), because it gives an eigenvalue of one when the modifier allele B2 has no effect on dominance:

SP[eigensoln /. h12 -> h11]

Collect[] allows us to gather together those terms involving μ:

Collect[eigensoln[[2]], {μ}, Factor]

which can be written in a way that it is easier to interpret the sign of the eigenvalue:

1 + (((h11 - h12) (2 r (1 - h12 s) + h12 s) μ)/(h11 (r (1 - h12 s) + h12 s)))

Thus, any modifier allele B2 that causes h12 < h11 will have an eigenvalue greater than one and so spread within the population.  Also observe that the leading eigenvalue is very near one when  is small.

This first approach is direct and fairly intuitive.  The disadvantages are:  (1) it takes Mathematica longer to calculate, (2) it can be harder to get the equations to simplify, and (3) it won't work in cases where the equilibrium condition cannot be solved for the equilibrium.

Approach 2:  Peforming a perturbation analysis

Alternatively a perturbation analysis can be used (Box 5.1). We already found in Chapter 8 that x[3] is approximately (μ/(h11 s)), but let's verify this fact.  First we solve eqcond = 0 for the constant order term by writing x[3] as a power series expandion in the mutation rate:  x[3] = x30+μ x31+(μ^2)x32:

SP[Normal[
        Series[eqcond /. 
            x[3] -> x30 + μ x31 + (μ^2) x32, {μ, 0, 
            0}]]]

Of the three potential solutions, the one at x30 -> 0 is the relevant equilibrium when we assume that the deleterious allele A2 is kept rare by selection:

Solve[% == 0, x30]

The next order term must then solve eqcond = 0 to linear order in the mutation rate:

SP[Normal[
        Series[eqcond /. 
              x[3] -> x30 + μ x31 + (μ^2) x32 /. 
            x30 -> 0, {μ, 0, 1}]]]

Solve[% == 0, x31]

which demonstrates that  x[3] = x30+μ x31+(μ^2)x32+... = (μ/(h11 s))to linear order in μ.

Next, we perform a Taylor series of the Det(M - I) = 0, which is the characteristic polynomial that the eigenvalue must solve.  Now, we write the eigenvalue λ as a power series expandion in the mutation rate:  λ = λ0+μ λ1+(μ^2)λ2 .  To be safe, we also write x[3] = (μ/(h11 s)) + (μ^2) x32, where we have not solved for the next term (x32) in the power series for the equilibrium, because it should be irrelevant, but we include it to make sure that it drops out.

To constant order, we have two potential solutions:

SP[Normal[
        Series[Det[submat - λ*IdentityMatrix[2] /. 
                  x[1] -> 1 - x[3] /. 
                x[3] -> (μ/(h11 s)) + (μ^2) 
                      x32 /. λ -> λ0 + μ λ1 + (μ^2) λ2], {μ, 0, 0}]]]

Simplify[Solve[% == 0, λ0]]

The relevant solution is λ0->1.  The next order term must then solve Det(M - I) = 0 to linear order in the mutation rate:

The next order term must then solve eqcond = 0 to linear order in the mutation rate:

SP[Normal[
        Series[Det[submat - λ*IdentityMatrix[2] /. 
                    x[1] -> 1 - x[3] /. 
                  x[3] -> (μ/(h11 s)) + (μ^2) 
                        x32 /. λ -> λ0 + μ λ1 + (μ^2) λ2 /. λ0 -> 1], {μ, 0, 
            1}]]]

Simplify[Solve[% == 0, λ1]]

Thus, to linear order in the mutation rate, the leading eigenvalue is λ=1+μ*\(\(\((h11 - h12)\)\ \((h12\ s + 2\ r\ \((1 - h12\ s)\))\)\)\/\(h11\ \((h12\ s + r\ \((1 - h12\ s)\))\)\)\), exactly as found in the previous section.  Again, we can conclude that mutant alleles B2 can spread if they reduce h12 relative to h11, causing the wildtype allele to become more dominant.

A perturbation approach is powerful and useful in many situations.  The main disadvantage is that it requires you to identify some parameter (such as the mutation rate) that is small.

Approach 3:  Using Recipe 12.5 to obtain the eigenvalue

Alternatively Box 13.1 can be used to write the eigenvalue as λ=1+(h12-h11)*(Subscript[((d(λ/dh12))), h12 = h11]).  First we must calculate the left and right eigenvectors of submat when h12 = h11:

submat /. x[1] -> 1 - x[3] /. h12 -> h11 // 
      Factor

These eigenvectors aren't going to be very useful, however, unless we evaluate them at the equilibrium.  Rather than solving the equilibrium condition for x[3], however, let's solve it for some term that appears linearly in eqcond:

Factor[eqcond]

μ is a good choice because we might think that x[3] and μ  provide similar information about the frequency of mutations.  h11 would be a bad choice, because we want to know whether a modifier that alters h12 away from h11 spreads, so we want to keep h11 in our calculations.

eqcond2 gives a relationship that must be satisfied by μ at the equilibrium:

eqcond2 = SP[Solve[% == 0, μ]] // Flatten

The mutant submatrix then becomes:

submatE = submat /. x[1] -> 1 - x[3] /. eqcond2 // 
        Factor

An eigenvalue correctly reduces to one when h12=h11:

Simplify[
      PowerExpand[Eigenvalues[submatE /. h12 -> h11]]]

At this point, we can find the left and right eigenvectors associated with the eigenvalue of one of the mutant submatrix when h12 = h11, as required by equation 13.1.1.  We do this by arbitrarily choosing the first element of the eigenvectors as one and finding the second element (recall that an eigenvector can be any length):

SP[Solve[(submatE /. h12 -> h11) . {1, u2} == {1, u2}, u2]] // 
      Flatten

Giving a right eigenvector of:

uvec = {1, x[3]/(1 - x[3])};

SP[Solve[{1, v2} . (submatE /. h12 -> h11) == {1, v2}, v2]] // 
      Flatten

Giving a left eigenvector of:

vvec = {1, (r (-1 + h11 s) (1 + 
                    2 h11 s (-1 + x[3]) x[3] - 
                    s x[3]^2))/(-s (-1 + 
                      h11 s x[3]) (-x[3] + 
                      h11 (-1 + 2 x[3])) + 
                r (-1 + 2 (h11^2) (s^2) (-1 + x[3]) x[3] + 
                      s x[3]^2 + 
                      h11 s (1 + 
                            2 x[3] - (2 + 
                                  s) x[3]^2)))};

Next, we have to calculate (Subscript[(Det((dM/dh12))), h12 = h11])and multiply on the left by vvec and on the right by uvec:

solution = 
      Factor[(vvec . (D[submatE, h12] /. h12 -> h11) . 
              uvec)/(vvec . uvec)]

The sign of this is not obvious yet, so let's rewrite it.  This can probably be done in many ways.  We just collected the parts involving r and not involving r in the numerator and denominator:

Collect[-2 r - h11 s + 2 h11 r s - s x[3] + 2 h11 s x[3], 
      r, Factor]

Collect[-r - h11 s + h11 r s - s x[3] + 3 h11 s x[3] + 
        2 h11 r s x[3] + (h11^2) (s^2) x[3] - 2 (h11^2) r (s^2) x[3] + 
        s x[3]^2 - 2 h11 s x[3]^2 + r s x[3]^2 - 
        2 h11 r s x[3]^2 + h11 (s^2) x[3]^2 - 
        3 (h11^2) (s^2) x[3]^2 - h11 r (s^2) x[3]^2 + 
        2 (h11^2) r (s^2) x[3]^2 - h11 (s^2) x[3]^3 + 
        2 (h11^2) (s^2) x[3]^3, r, Factor]

Thus, we can rewrite "solution" as:

(s (-1 + x[3]) x[
              3] (2 r (-1 + h11 s) + 
                s (-h11 - x[3] + 
                      2 h11 x[3])))/(s (-1 + 
                  x[3]) (-h11 - x[3] + 2 h11 x[3]) (-1 + 
                  h11 s x[3]) + 
            r (-1 + h11 s) (1 - 2 h11 s x[3] - s x[3]^2 + 
                  2 h11 s x[3]^2));

Next, we carefully change the sign of things like -1 + x[3], which is negative, to 1 - x[3], which is positive.  This helps to evaluate the sign of this equation.

-(s (1 - x[3]) x[
                3] (2 r (1 - h11 s) + 
                  s (h11 (1 - x[3]) + 
                        x[3] (1 - h11))))/(s (1 - 
                  x[3]) (h11 (1 - x[3]) + 
                  x[3] (1 - h11)) (1 - h11 s x[3]) + 
            r (1 - h11 s) (1 - 2 h11 s x[3] - s x[3]^2 + 
                  2 h11 s x[3]^2));

Just to make sure, this does still equal "solution":

Factor[%/solution]

The only term whose sign is not obvious is (1 - 2 h11 s x[3] - s x[3]^2 + 2 h11 s x[3]^2).  It can help to rewrite this as (1 - 2 h11 s x[3] (1 - x[3]) - 
              s x[3]^2) .   Because x[3] equals q, the frequency of the deleterious allele, we can write this as (1 - 2 p q  h11 s - s (q^2)).  At this point, you might recognize this as the mean fitness, ((p^2) + 
          2 p q (1 -  h11 s) +  (q^2) (1 - s)) .  Admittedly, this takes insight into the model (and some lateral thinking).

Using the mean fitness:

wbar = ((1 - x[3])^2 + (1 - h11 s) 
              2 x[3] (1 - x[3]) + (1 - s) x[3]^2);

we can then rewrite solution as an equation whose sign is negative (because every term is positive except for the "-" in front):

-(s (1 - x[3]) x[
                3] (2 r (1 - h11 s) + 
                  s (h11 (1 - x[3]) + 
                        x[3] (1 - h11))))/(s (1 - 
                  x[3]) (h11 (1 - x[3]) + 
                  x[3] (1 - h11)) (1 - h11 s x[3]) + 
            r (1 - h11 s) wbar);

Factor[%/solution]

The above gives us (Subscript[((d(λ/dh12))), h12 = h11]), which we can plug into λ=1+(h12-h11)*(Subscript[((d(λ/dh12))), h12 = h11]).  As (Subscript[((d(λ/dh12))), h12 = h11])is negative, the leading eigenvalue will be greater than one and the B2 allele will spread when h12 < h11.  Yet again, we find that dominance evolves such that the wildtype allele A1 becomes more dominant than the mutant A2 allele.

Using Recipe 12.5 seems rather complicated in this model.  It is important to point out, however, that we have not used the fact that the mutation rate is small in the above calculations.  Rather, we have only assumed that the modifier allele B2 has only a small effect so that the eigenvalue can be approximated by 1+(h12-h11)*(Subscript[((d(λ/dh12))), h12 = h11])(i.e., terms such as (((h12 - h11))^2) have been ignored).

To compare the result of Approach 3 to that obtained from Approaches 1 and 2, we take the Taylor series of (Subscript[((d(λ/dh12))), h12 = h11]) with respect to x[3], which is small if the mutation rate,  μ, is small:

Normal[Series[solution, {x[3], 0, 1}]]

which indicates that the eigenvalue is λ=1-(h12-h11)*\(TraditionalForm\`\(s\ \((2\ \(r(1 - h11\ s)\) + h11\ s)\)\ x[3]\)\/\(r(1 - h11\ s) + h11\ s\)\), exactly as found in the previous sections once we substitute the approximate equilibrium,  x[3] = (μ/(h11 s)).

Other approaches

The above methods are not the only ones that can be used to analyze this and similar models.  For example, the graphical method outlined in Sup. Mat. 8.1 can be used to determine when the leading eigenvalue is greater than one.  Alternatively, you could take the direct Approach 1, but solve the equilibrium condition for μ rather than x[3] (as we did in Approach 3 above).  Or you could perform a perturbation analysis assuming that (h12 - h11) is small rather than μ is small.  If you are curious, try one (or more) of these approaches to see for yourself which one(s) work best and what assumptions are required for each.

As demonstrated by this lab, there are often alternative routes to the same main result.  Not all of them will work in any particular model, so it is worth knowing that there are other methods if the first method you try leads to a dead end (or a hopelessly complicated end).