Generating the Figures in Chapter 9 of Otto and Day

Inputs (Enter first)

<< Graphics`Graphics`

<< Graphics`MultipleListPlot`

<< Graphics`Arrow`

Off[General::"spell"]

Off[General::"spell1"]

Needs["Graphics`PlotField`"]

Figure 9.1 & 9.2 - Change of bases to eigenvectors

matA = {{ 4, 4}, {1, 12}};

dmat = {{2/3, 0}, {0, 4/3}}

mat = matA . dmat . Inverse[matA]

Eigensystem[mat]

axisplot is just used to generate the proper axes and plot range:

PlotRangeVectors = {{0, 5}, {0, 12.5}};

vec1 = {4, 1};
vec2 = {4, 12};
vec3 = {2, 1};

DisplayTogether[
      ListPlot[Table[MatrixPower[mat, t] . vec3, {t, 0, 10}], 
        PlotStyle -> PointSize[0.018]], 
      Graphics[{Arrow[{0, 0}, vec1], 
          Arrow[{0, 0}, vec2]}], 
PlotRange -> 
        PlotRangeVectors]

Inverse[matA] . vec3

DisplayTogether[
      ListPlot[Table[
          MatrixPower[dmat, t] . Inverse[matA] . vec3, {t, 0, 10}], 
        PlotStyle -> 
          PointSize[0.018]], 
PlotRange -> {{0, 
            1}, {0, 1}}]

Figure 9.4 - Kimura's model

mat = {{1 - a - 2  b, a, b, b}, {a, 1 - a - 2  b, b, b}, {b, b, 
            1 - a - 2  b, a}, {b, b, a, 1 - a - 2  b}};

Eigensystem[mat]

amat = Transpose[Eigenvectors[mat]]
dmat = Eigenvalues[mat]^t*IdentityMatrix[4]

gensol = amat . dmat . Inverse[amat] . {x10, x20, x30, x40}

Starting from an adenine, the number of A, G, C and T vary over time according to:

numsoln = (gensol /. a ->  10^(-7) /. 
            b ->  10^(-8) /. {x10 -> 1, x20 -> 0, 
            x30 -> 0, x40 -> 0})

The probability that two sequences that were both A at time 0 differ at time t is:

probdiffer = 
        1 - (numsoln[[1]]^2 + numsoln[[2]]^2 + numsoln[[3]]^2 + 
              numsoln[[4]]^2);

Given that they differ, the fraction that differ by a transition is:

probtransition = (2 numsoln[[1]]*numsoln[[2]] + 
              2 numsoln[[3]]*numsoln[[4]])/
          probdiffer;

Given that they differ, the fraction that differ by a transversion is:

probtransversion = (2 numsoln[[1]]*numsoln[[3]] + 
              2 numsoln[[1]]*numsoln[[4]] + 
              2 numsoln[[2]]*numsoln[[3]] + 
              2 numsoln[[2]]*numsoln[[4]])/
          probdiffer;

The sum of these two probabilities is one, because we conditioned on the fact that they differed:

Simplify[probtransition + probtransversion]

endtime = 50000000. ;

DisplayTogether[
Plot[
        probdiffer, {t, 0, endtime}], 
Plot[
        probtransition, {t, 0, endtime}], 
Plot[
        probtransversion, {t, 0, 
          endtime}], 
PlotRange -> {{0, 
            1.05 endtime}, {0, 
            1}}
]

Figure 9.5 - Density-dependent natural selection

w[n_, r_, α_] := 1 + r (1 - α n)

eqn1 = 
      n[t] (p[t]^2 w[n[t], r, αAA] + 
            2 p[t] (1 - p[t]) w[n[t], 
                r, αAa] + (1 - p[t])^2 w[n[t], r, αaa]); 
    eqn2 = p[t + 1] = 
        p[t]*(p[t] w[n[t], r, αAA] + (1 - p[t]) w[n[t], 
                    r, αAa])/(p[t]^2 w[n[t], r, αAA] + 
                2 p[t] (1 - p[t]) w[n[t], 
                    r, αAa] + (((1 - p[t]))^2) w[n[t], 
                    r, αaa]);

Factor[eqn1 - n[t]]

Solve[% == 0, n[t]]

Factor[αaa - 2 αaa p[t] + 
        2 αAa p[t] + αaa p[t]^2 - 
        2 αAa p[t]^2 + αAA p[t]^2 - (αaa (1 - 
                    p[t])^2 + 
            2 αAa p[
                t] (1 - p[t]) + αAA p[t]^2)]

subn = {n[t] -> 
            (1/((αaa (1 - p[t])^2 + 2 αAa p[t] (1 - p[t]) + αAA p[t]^2)))};

Factor[eqn2 - p[t] /. subn]

This is similar to the non-frequency dependent model of selection multiplied by -r:

subW = {Waa ->  αaa, WAa -> αAa, 
          WAA -> αAA};

(WAA p[t]^2 + WAa p[t] (1 - p[t]))/(WAA p[t]^2 + 
                WAa 2 p[t] (1 - p[t]) + Waa (1 - p[t])^2) - 
          p[t] /. subW // Factor

%%%/% // Factor

subsmall = {αaa ->  αbar (1 + 
                  daa eps), αAa ->  αbar (1 + 
                  dAa eps), αAA ->  αbar (1 + 
                  dAA eps)};

Factor[Normal[Series[eqn1 /. subsmall, {eps, 0, 0}]]]

Factor[Normal[Series[eqn2 /. subsmall, {eps, 0, 0}]]]

(There's no choice of constant fitnesses that will make them equal, however.)

Figure (a):

n[1] = 1800; p[1] :=  .05; 
    papprox[1] = 0.05;
αAA =  .0004; αAa =  .00045; αaa =  .0005;
r = 0.5;
{1 + r, Exp[2]} // N
{αAA, αAa} // N
tmax = 200;
For[t = 1, 
      t ≤ tmax, t++, 
n[t + 1] = 
        n[t] (p[t]^2 w[n[t], r, αAA] + 
              2 p[t] (1 - p[t]) w[n[t], 
                  r, αAa] + (1 - p[t])^2 w[n[t], 
                  r, αaa]); 
      p[t + 1] = 
        p[t]*(p[t] w[n[t], r, αAA] + (1 - p[t]) w[n[t], 
                    r, αAa])/(p[t]^2 w[n[t], r, αAA] + 
                2 p[t] (1 - p[t]) w[n[t], 
                    r, αAa] + (((1 - p[t]))^2) w[n[t], 
                    r, αaa]); 
papprox[t + 1] = 
        papprox[t] + (-r) ((1 - papprox[t]) papprox[
                    t] ((1 - 
                            papprox[t]) (αAa  - αaa) + 
                      papprox[
                          t] (αAA  - αAa)))/(αaa (((1 - papprox[t]))^2) + 
                  2 αAa (1 - papprox[t]) papprox[
                      t] + αAA papprox[t]^2)]

plot1 = 
      DisplayTogether[
ListPlot[Table[n[i], {i, 1, tmax}], 
          PlotJoined -> True, 
          PlotStyle -> Thickness[0.01]], 
ListPlot[
          Table[3600 p[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[1/(αaa (1 - papprox[i])^2 + 
                  2 αAa papprox[
                      i] (1 - 
                        papprox[i]) + αAA papprox[i]^2), {i, 1, 
              tmax}], PlotJoined -> True, 
          PlotStyle -> {Thickness[0.01], 
              Dashing[{ .02}]}], 
ListPlot[
          Table[3600 papprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{ .02}]], 
PlotRange -> {{0, 
              160}, {0, 3600}}
]

Figure (b):

n[1] = 1800; p[1] :=  .05; 
    papprox[1] = 0.05;
αAA =  .0004; αAa =  .00045; αaa =  .0005;
r = 0.05;
{1 + r, Exp[2]} // N
{αAA, αAa} // N
tmax = 2000;
For[t = 1, 
      t ≤ tmax, t++, 
n[t + 1] = 
        n[t] (p[t]^2 w[n[t], r, αAA] + 
              2 p[t] (1 - p[t]) w[n[t], 
                  r, αAa] + (1 - p[t])^2 w[n[t], 
                  r, αaa]); 
      p[t + 1] = 
        p[t]*(p[t] w[n[t], r, αAA] + (1 - p[t]) w[n[t], 
                    r, αAa])/(p[t]^2 w[n[t], r, αAA] + 
                2 p[t] (1 - p[t]) w[n[t], 
                    r, αAa] + (((1 - p[t]))^2) w[n[t], 
                    r, αaa]); 
papprox[t + 1] = 
        papprox[t] + (-r) ((1 - papprox[t]) papprox[
                    t] ((1 - 
                            papprox[t]) (αAa  - αaa) + 
                      papprox[
                          t] (αAA  - αAa)))/(αaa (((1 - papprox[t]))^2) + 
                  2 αAa (1 - papprox[t]) papprox[
                      t] + αAA papprox[t]^2)]

plot2 = 
      DisplayTogether[
ListPlot[Table[n[i], {i, 1, tmax}], 
          PlotJoined -> True, 
          PlotStyle -> Thickness[0.01]], 
ListPlot[
          Table[3600 p[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[1/(αaa (1 - papprox[i])^2 + 
                  2 αAa papprox[
                      i] (1 - 
                        papprox[i]) + αAA papprox[i]^2), {i, 1, 
              tmax}], PlotJoined -> True, 
          PlotStyle -> {Thickness[0.01], 
              Dashing[{ .02}]}], 
ListPlot[
          Table[3600 papprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{ .015}]], 
PlotRange -> {{0, 
              1600}, {0, 3600}}
]

Figure (c):

n[1] = 300; p[1] :=  .05; papprox[1] = 0.05;
αAA =  .0004; αAa =  .00045; αaa =  .0005;
r = 0.05;
{1 + r, Exp[2]} // N
{αAA, αAa} // N
tmax = 2000;
For[t = 1, 
      t ≤ tmax, t++, 
n[t + 1] = 
        n[t] (p[t]^2 w[n[t], r, αAA] + 
              2 p[t] (1 - p[t]) w[n[t], 
                  r, αAa] + (1 - p[t])^2 w[n[t], 
                  r, αaa]); 
      p[t + 1] = 
        p[t]*(p[t] w[n[t], r, αAA] + (1 - p[t]) w[n[t], 
                    r, αAa])/(p[t]^2 w[n[t], r, αAA] + 
                2 p[t] (1 - p[t]) w[n[t], 
                    r, αAa] + (((1 - p[t]))^2) w[n[t], 
                    r, αaa]); 
papprox[t + 1] = 
        papprox[t] + (-r) ((1 - papprox[t]) papprox[
                    t] ((1 - 
                            papprox[t]) (αAa  - αaa) + 
                      papprox[
                          t] (αAA  - αAa)))/(αaa (((1 - papprox[t]))^2) + 
                  2 αAa (1 - papprox[t]) papprox[
                      t] + αAA papprox[t]^2)]

plot3 = 
      DisplayTogether[
ListPlot[Table[n[i], {i, 1, tmax}], 
          PlotJoined -> True, 
          PlotStyle -> Thickness[0.01]], 
ListPlot[
          Table[3600 p[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[1/ (αaa (1 - papprox[i])^2 + 
                  2 αAa papprox[
                      i] (1 - 
                        papprox[i]) + αAA papprox[i]^2), {i, 1, 
              tmax}], PlotJoined -> True, 
          PlotStyle -> {Thickness[0.01], 
              Dashing[{ .02}]}], 
ListPlot[
          Table[3600 papprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{ .02}]], 
PlotRange -> {{0, 
              1600}, {0, 3600}}
]

Figure (d):

n[1] = 1400; p[1] :=  .05; 
    papprox[1] = 0.05;
αAA =  .00015; αAa =  .00045; αaa =  .00075;
r = 0.05;
{1 + r, Exp[2]} // N
{αAA, αAa} // N
tmax = 500;
For[t = 1, 
      t ≤ tmax, t++, 
n[t + 1] = 
        n[t] (p[t]^2 w[n[t], r, αAA] + 
              2 p[t] (1 - p[t]) w[n[t], 
                  r, αAa] + (1 - p[t])^2 w[n[t], 
                  r, αaa]); 
      p[t + 1] = 
        p[t]*(p[t] w[n[t], r, αAA] + (1 - p[t]) w[n[t], 
                    r, αAa])/(p[t]^2 w[n[t], r, αAA] + 
                2 p[t] (1 - p[t]) w[n[t], 
                    r, αAa] + (((1 - p[t]))^2) w[n[t], 
                    r, αaa]); 
papprox[t + 1] = 
        papprox[t] + (-r) ((1 - papprox[t]) papprox[
                    t] ((1 - 
                            papprox[t]) (αAa  - αaa) + 
                      papprox[
                          t] (αAA  - αAa)))/(αaa (((1 - papprox[t]))^2) + 
                  2 αAa (1 - papprox[t]) papprox[
                      t] + αAA papprox[t]^2)]

plot4 = 
      DisplayTogether[
ListPlot[Table[n[i], {i, 1, tmax}], 
          PlotJoined -> True, 
          PlotStyle -> Thickness[0.01]], 
ListPlot[
          Table[10000 p[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[1/(αaa (1 - papprox[i])^2 + 
                  2 αAa papprox[
                      i] (1 - 
                        papprox[i]) + αAA papprox[i]^2), {i, 1, 
              tmax}], PlotJoined -> True, 
          PlotStyle -> {Thickness[0.01], 
              Dashing[{ .02}]}], 
ListPlot[
          Table[10000 papprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{ .02}]], 
PlotRange -> {{0, 
              300}, {0, 10000}}
]

Display["graph1.eps", plot1, "EPS"]

Display["graph2.eps", plot2, "EPS"]

Display["graph3.eps", plot3, "EPS"]

Display["graph4.eps", plot4, "EPS"]

Figure 9.6 & 9.7 -  Two-locus model

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);
eqn1 = 
        x[1] Sum[
              x[j] (w[1, j] + w[j, 1])/(2 wbar), {j, 1, 4}]  -  
          r dis;
eqn2 = 
        x[2] Sum[
              x[j] (w[2, j] + w[j, 2])/(2 wbar), {j, 1, 4}]  +  
          r dis;
eqn3 = 
        x[3] Sum[
              x[j] (w[3, j] + w[j, 3])/(2 wbar), {j, 1, 4}]  +  
          r dis;
eqn4 = 
        x[4] Sum[
              x[j] (w[4, j] + w[j, 4])/(2 wbar), {j, 1, 4}]  -  
          r dis;

w[i_, j_] = w[i] w[j];

backtransform = {pA -> x[1] + x[2], pB -> x[1] + x[3], 
          d -> x[1]*x[4] - x[2]*x[3]};

Simplify[
      Solve[{pA == x[1] + x[2], pB == x[1] + x[3], 
          d == x[1]*x[4] - x[2]*x[3], 
          x[1] + x[2] + x[3] + x[4] == 1}, {x[1], x[2], x[3], 
          x[4]}]]

transform = {x[4] -> d + (-1 + pA) (-1 + pB), 
          x[2] -> -d + pA - pA pB, 
          x[3] -> -d + pB - pA pB, 
          x[1] -> d + pA pB};

wbarhaploid = 
      Simplify[Sum[x[i]  w[i], {i, 1, 4}] /. transform]

subfit = {w[1] -> 1 + s1*eps, w[2] -> 1 + s2*eps, 
          w[3] -> 1 + s3*eps, w[4] -> 1 + s4*eps};

eqnpA = Simplify[eqn1 + eqn2 /. transform]

Simplify[%*wbarhaploid]

eqnpB = Simplify[eqn1 + eqn3 /. transform]

Simplify[%*wbarhaploid]

eqnD = Simplify[eqn1 eqn4 - eqn2 eqn3 /. transform]

Simplify[%*wbarhaploid^2]

Factor[% -  (1 - 
              r) (w[1] w[
                  4] (pA pB + d) ( (1 - pA) (1 - pB) + 
                    d) - w[2] w[3] ((1 - pA) pB  - d) ( 
                  pA (1 - pB) - d))]

At the quasi-equilibrium, the disequilibrium is zero to O(1), in which case the allele frequency dynamics are:

Factor[eqnpA /. d -> 0]

Factor[% - (( pA WA)/(WA pA + Wa (1 - pA))) /. 
          WA -> pB w[1] + (1 - pB) w[2] /. 
        Wa -> pB w[3] + (1 - pB) w[4]]

The same picture emerges if we also substitute in for the siftnesses:

Factor[
      Normal[Series[eqnpA - pA /. d -> d*eps /. subfit, {eps, 0, 
            1}]]]

Factor[
      Normal[Series[% - pA (1 - pA) (WA -  Wa) /. 
                WA -> pB w[1] + (1 - pB) w[2] /. 
              Wa -> pB w[3] + (1 - pB) w[4] /. subfit, {eps, 0, 
            2}]]]

Here we do a perturbation analysis, setting d=d0+d1*eps+d2*eps^2 and obtaining successively more accurate approximations:

Factor[
        Normal[Series[eqnD - d /. 
                  w[4] -> (epi*eps + w[2] w[3])/w[1] /. subfit /. 
              d -> d0 + d1*eps, {eps, 0, 0}]]] /. 
      eps -> 1

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

Factor[
        Normal[Series[eqnD - d /. 
                    w[4] -> (epi*eps + w[2] w[3])/w[1] /. 
                  subfit /. d -> d0 + d1*eps /. d0 -> 0, {eps, 
              0, 1}]]] /. eps -> 1

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

Assuming weaker epistasis:

Factor[
        Normal[Series[eqnD - d /. 
                    w[4] -> (epi*eps^2 + w[2] w[3])/w[1] /. 
                  subfit /. d -> d0 + d1*eps + d2*eps^2 /. 
              d0 -> 0, {eps, 0, 1}]]] /. eps -> 1

Factor[
        Normal[Series[eqnD - d /. 
                      w[4] -> (epi*eps^2 + w[2] w[3])/w[1] /. 
                    subfit /. d -> d0 + d1*eps + d2*eps^2 /. 
                d0 -> 0 /. d1 -> 0, {eps, 0, 2}]]] /. 
      eps -> 1

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

dQLE = ((epi (1 - pA) pA (1 - pB) pB (1 - r))/(r ));

Factor[
      Normal[Series[eqnD - dQLE /. d -> delta + dQLE /. 
                epi -> epi*eps /. 
              w[4] -> (epi*eps + w[2] w[3])/w[1] /. subfit, {eps, 
            0, 0}]]]

If we wanted a more accurate prediction for the allele frequency dynamics, we could use

Factor[
      Normal[Series[eqnpA - pA /. d -> dQLE /. 
              epi -> w[1] w[4] - w[2] w[3] /. subfit, {eps, 0, 
            1}]]]

Simplify[
      SeriesCoefficient[
        Series[eqnpA - pA /. d -> d1*eps /. 
              epi -> w[1] w[4] - w[2] w[3] /. subfit, {eps, 0, 2}], 
        2]]

Simplify[% /. d1 -> 0]

Simplify[%% - %]

pA[1] = 0.05; pB[1] = 0.05; pAapprox[1] = 0.05; pBapprox[1] = 0.05; 
    d[1] = 0;
w[1] = 1.05; w[2] = 1.01; w[3] = 1.01; 
    w[4] = 1;
r = 0.1;
tmax = 400;
For[t = 1, 
      t ≤ 
        tmax, t++, 
whap[
          t] = (pA[t] pB[t] + d[t]) w[
              1] + (pA[t] (1 - pB[t]) - d[t]) w[
              2] + ((1 - pA[t]) pB[t] - d[t]) w[
              3] + ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
              4]; 
pA[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + (pA[t] (1 - pB[t]) - d[t]) w[2])/
          whap[t]; 
pB[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + ((1 - pA[t]) pB[t] - d[t]) w[3])/
          whap[t]; 
d[
          t + 1] = (1 - 
              r) ((pA[t] pB[t] + d[t]) w[
                    1] ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
                    4] - (pA[t] (1 - pB[t]) - d[t]) w[
                    2] ((1 - pA[t]) pB[t] - d[t]) w[3])/
            whap[t]^2; 
pAapprox[t + 1] = 
        pAapprox[t] + 
          pAapprox[
              t] (1 - 
                pAapprox[
                  t]) (pBapprox[
                    t] (w[1] - w[3]) + (1 - pBapprox[t]) (w[2] - 
                      w[4])); 
pBapprox[t + 1] = 
        pBapprox[t] + 
          pBapprox[t] (1 - pBapprox[t]) ( 
              pAapprox[
                    t] (w[1] - w[2]) + (1 - pAapprox[t]) (w[3] - 
                      w[4]))
]

plot1a = 
      DisplayTogether[
ListPlot[Table[pA[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[pAapprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{0.02}]], 
PlotRange -> {{0, 
              tmax}, {0, 1}}
]

plot1b = 
      DisplayTogether[
ListPlot[Table[d[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[(((w[1] w[4] - w[2] w[3]) (1 - pAapprox[i]) pAapprox[i] (1 - pBapprox[i]) pBapprox[i] (1 - r))/(r )), {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> {Dashing[{0.02}]}], 
PlotRange -> {{0, tmax}, {0, 0.021}}
]

pA[1] = 0.05; pB[1] = 0.05; pAapprox[1] = 0.05; pBapprox[1] = 0.05; 
    d[1] = 0;
w[1] = 1.05; w[2] = 1.01; w[3] = 1.01; 
    w[4] = 1;
r = 0.5;
tmax = 400;
For[t = 1, 
      t ≤ 
        tmax, t++, 
whap[
          t] = (pA[t] pB[t] + d[t]) w[
              1] + (pA[t] (1 - pB[t]) - d[t]) w[
              2] + ((1 - pA[t]) pB[t] - d[t]) w[
              3] + ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
              4]; 
pA[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + (pA[t] (1 - pB[t]) - d[t]) w[2])/
          whap[t]; 
pB[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + ((1 - pA[t]) pB[t] - d[t]) w[3])/
          whap[t]; 
d[
          t + 1] = (1 - 
              r) ((pA[t] pB[t] + d[t]) w[
                    1] ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
                    4] - (pA[t] (1 - pB[t]) - d[t]) w[
                    2] ((1 - pA[t]) pB[t] - d[t]) w[3])/
            whap[t]^2; 
pAapprox[t + 1] = 
        pAapprox[t] + 
          pAapprox[
              t] (1 - 
                pAapprox[
                  t]) (pBapprox[
                    t] (w[1] - w[3]) + (1 - pBapprox[t]) (w[2] - 
                      w[4])); 
pBapprox[t + 1] = 
        pBapprox[t] + 
          pBapprox[t] (1 - pBapprox[t]) ( 
              pAapprox[
                    t] (w[1] - w[2]) + (1 - pAapprox[t]) (w[3] - 
                      w[4]))
]

plot2a = 
      DisplayTogether[
ListPlot[Table[pA[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[pAapprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{0.02}]], 
PlotRange -> {{0, 
              tmax}, {0, 1}}
]

plot2b = 
      DisplayTogether[
ListPlot[Table[d[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[(((w[1] w[4] - w[2] w[3]) (1 - pAapprox[i]) pAapprox[i] (1 - pBapprox[i]) pBapprox[i] (1 - r))/(r )), {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> {Dashing[{0.02}]}], 
PlotRange -> {{0, tmax}, {0, 0.0021}}
]

pA[1] = 0.05; pB[1] = 0.05; pAapprox[1] = 0.05; pBapprox[1] = 0.05; 
    d[1] = 0;
w[1] = 1.005; w[2] = 1.001; w[3] = 1.001; 
    w[4] = 1;
r = 0.1;
tmax = 5000;
For[t = 1, 
      t ≤ 
        tmax, t++, 
whap[
          t] = (pA[t] pB[t] + d[t]) w[
              1] + (pA[t] (1 - pB[t]) - d[t]) w[
              2] + ((1 - pA[t]) pB[t] - d[t]) w[
              3] + ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
              4]; 
pA[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + (pA[t] (1 - pB[t]) - d[t]) w[2])/
          whap[t]; 
pB[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + ((1 - pA[t]) pB[t] - d[t]) w[3])/
          whap[t]; 
d[
          t + 1] = (1 - 
              r) ((pA[t] pB[t] + d[t]) w[
                    1] ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
                    4] - (pA[t] (1 - pB[t]) - d[t]) w[
                    2] ((1 - pA[t]) pB[t] - d[t]) w[3])/
            whap[t]^2; 
pAapprox[t + 1] = 
        pAapprox[t] + 
          pAapprox[t] (1 - pAapprox[t]) ( 
              pB[t] (w[1] - w[3]) + (1 - pB[t]) (w[2] - 
                      w[4])); 
pBapprox[t + 1] = 
        pBapprox[t] + 
          pBapprox[t] (1 - pBapprox[t]) ( 
              pA[t] (w[1] - w[2]) + (1 - pA[t]) (w[3] - 
                      w[4]))
]

plot3a = 
      DisplayTogether[
ListPlot[Table[pA[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[pAapprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{0.02}]], 
PlotRange -> {{0, 
              tmax}, {0, 1}}
]

plot3b = 
      DisplayTogether[
ListPlot[Table[d[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[(((w[1] w[4] - w[2] w[3]) (1 - pAapprox[i]) pAapprox[i] (1 - pBapprox[i]) pBapprox[i] (1 - r))/(r )), {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> {Dashing[{0.02}]}], 
PlotRange -> {{0, tmax}, {0, 0.0021}}
]

pA[1] = 0.05; pB[1] = 0.05; pAapprox[1] = 0.05; pBapprox[1] = 0.05; 
    d[1] = 0;
w[1] = 1.0015; w[2] = 1.001; w[3] = 1.001; 
    w[4] = 1;
r = 0.1;
tmax = 10000;
For[t = 1, 
      t ≤ 
        tmax, t++, 
whap[
          t] = (pA[t] pB[t] + d[t]) w[
              1] + (pA[t] (1 - pB[t]) - d[t]) w[
              2] + ((1 - pA[t]) pB[t] - d[t]) w[
              3] + ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
              4]; 
pA[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + (pA[t] (1 - pB[t]) - d[t]) w[2])/
          whap[t]; 
pB[
          t + 1] = ((pA[t] pB[t] + d[t]) w[
                  1] + ((1 - pA[t]) pB[t] - d[t]) w[3])/
          whap[t]; 
d[
          t + 1] = (1 - 
              r) ((pA[t] pB[t] + d[t]) w[
                    1] ((1 - pA[t]) (1 - pB[t]) + d[t]) w[
                    4] - (pA[t] (1 - pB[t]) - d[t]) w[
                    2] ((1 - pA[t]) pB[t] - d[t]) w[3])/
            whap[t]^2; 
pAapprox[t + 1] = 
        pAapprox[t] + 
          pAapprox[t] (1 - pAapprox[t]) ( 
              pB[t] (w[1] - w[3]) + (1 - pB[t]) (w[2] - 
                      w[4])); 
pBapprox[t + 1] = 
        pBapprox[t] + 
          pBapprox[t] (1 - pBapprox[t]) ( 
              pA[t] (w[1] - w[2]) + (1 - pA[t]) (w[3] - 
                      w[4]))
]

plot4a = 
      DisplayTogether[
ListPlot[Table[pA[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[pAapprox[i], {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> 
            Dashing[{0.02}]], 
PlotRange -> {{0, 
              tmax}, {0, 1}}
]

plot4b = 
      DisplayTogether[
ListPlot[Table[d[i], {i, 1, tmax}], 
          PlotJoined -> True], 
ListPlot[
          Table[(((w[1] w[4] - w[2] w[3]) (1 - pAapprox[i]) pAapprox[i] (1 - pBapprox[i]) pBapprox[i] (1 - r))/(r )), {i, 1, tmax}], PlotJoined -> True, 
          PlotStyle -> {Dashing[{0.02}]}], 
PlotRange -> {{0, tmax}, {0, -0.00041}}
]

Display["graph.eps", plot1b, "EPS"]

Display["graph2.eps", plot2b, "EPS"]

Display["graph4.eps", plot3b, "EPS"]

Display["graph3.eps", plot4b, "EPS"]

Figure 9.2.1 - Complex plane

PlotRangeVectors = {{-4, 4}, {-4, 4}};

vec1 = {3, 2};

360/(2 Pi) ArcTan[2/3] // N

DisplayTogether[ListPlot[{vec1}, PlotRange -> PlotRangeVectors], 
      Graphics[{Arrow[{0, 0}, vec1]}], 
PlotRange -> 
        PlotRangeVectors]

Figure 9.3.1 - Cyclic dynamics

solncomplex =  {n1[
              t] -> (E^((1/2) (a + d) t)) (n20 (( b Sin[beta t])/beta) + 
                  n10 ( 
                      Cos[beta t] + (((a - d) Sin[beta t])/(2 beta)))), 
          n2[t] -> (E^((1/2) (a + d) t)) (n10 ((c Sin[beta t])/beta) + 
                  n20 ( 
                      Cos[beta t] + (((-a + d) Sin[beta t])/(2 beta))))};

sub = {Gp -> Random[], Gz -> Random[], c -> Random[], 
          b -> Random[], B -> Random[], 
          a -> Random[]};
4 Pi/Sqrt[-(a - d)^2 - 4 b c] /. {a -> - Gz c, 
          b -> Gz a  -  b B, c -> -B c, 
          d -> -Gp b  +  a B} /. 
      sub
(a + d)/2 /. {a -> - Gz c, b -> Gz a  -  b B, 
          c -> -B c, d -> -Gp b  +  a B} /. 
      sub
sub

sub = {Gp -> 0.8, Gz -> 0.15, c -> 0.45, 
          b -> 0.3, B -> 0.32, 
          a -> 0.95};
4 Pi/Sqrt[-(a - d)^2 - 4 b c] /. {a -> - Gz c, 
          b -> Gz a  -  b B, c -> -B c, 
          d -> -Gp b  +  a B} /. 
      sub
(a + d)/2 /. {a -> - Gz c, b -> Gz a  -  b B, 
          c -> -B c, d -> -Gp b  +  a B} /. 
      sub
sub

{- Gz c, Gz a  -  b B, -B c, -Gp b  +  a B} /. 
      sub

range = 3; 
    DisplayTogether[
      ParametricPlot[{n1[t], n2[t]} /. solncomplex /. 
                  n10 -> 1 /. n20 -> 0 /. 
              beta -> 
                Sqrt[-(a - d)^2 - 4 b c]/
                  2 /. {a -> - Gz c, b -> Gz a  -  b B, 
              c -> -B c, d -> -Gp b  +  a B} /. 
          sub, {t, 0, 500}, 
        PlotRange -> {{-range, range}, {-range, 
              range}}], 
ListPlot[
        Table[{n1[t], n2[t]} /. solncomplex /. n10 -> 1 /. 
                  n20 -> 0 /. 
                beta -> 
                  Sqrt[-(a - d)^2 - 4 b c]/
                    2 /. {a -> - Gz c, 
                b -> Gz a  -  b B, c -> -B c, 
                d -> -Gp b  +  a B} /. sub, {t, 0, 500, 50}], 
        PlotStyle -> PointSize[0.02]]
]

Figure S9.1 - Epistasis

case1 = 
      Solve[{(Exp[alpha*x + beta*x^2] /. x -> 1) == 
              1.2, (Exp[alpha*x + beta*x^2] /. x -> 2) == 
              1.44}, {alpha, beta}] // Flatten

case2 = 
      Solve[{(Exp[alpha*x + beta*x^2] /. x -> 1) == 
              1.2, (Exp[alpha*x + beta*x^2] /. x -> 2) == 
              1.3}, {alpha, beta}] // Flatten

case3 = 
      Solve[{(Exp[alpha*x + beta*x^2] /. x -> 1) == 
              1.2, (Exp[alpha*x + beta*x^2] /. x -> 2) == 
              1.58}, {alpha, beta}] // Flatten

plot1 = 
      DisplayTogether[
ListPlot[{{0, 1}, {1, 1.2}, {2, 
              1.44}}, PlotStyle -> 
            PointSize[0.025]], 
Plot[
          Exp[alpha*x + beta*x^2] /. case1, {x, 0, 2}, 
          PlotStyle -> 
            Thickness[0.01]], 
ListPlot[{{0, 1}, {1, 
              1.2}, {2, 1.3}}, 
          PlotStyle -> PointSize[0.025]], 
Plot[
          Exp[alpha*x + beta*x^2] /. case2, {x, 0, 
            2}], 
ListPlot[{{0, 1}, {1, 1.2}, {2, 1.58}}, 
          PlotStyle -> PointSize[0.025]], 
Plot[
          Exp[alpha*x + beta*x^2] /. case3, {x, 0, 2}, 
          PlotStyle -> 
            Dashing[{0.01, 
                0.01}]], 
PlotRange -> {{-0.1, 
              2.1}, {0.9, 1.6}}
]

case1 = 
      Solve[{(Exp[alpha*x + beta*x^2] /. x -> 1) == 
              0.8, (Exp[alpha*x + beta*x^2] /. x -> 2) == 
              0.64}, {alpha, beta}] // Flatten

case2 = 
      Solve[{(Exp[alpha*x + beta*x^2] /. x -> 1) == 
              0.8, (Exp[alpha*x + beta*x^2] /. x -> 2) == 
              0.78}, {alpha, beta}] // Flatten

case3 = 
      Solve[{(Exp[alpha*x + beta*x^2] /. x -> 1) == 
              0.8, (Exp[alpha*x + beta*x^2] /. x -> 2) == 
              0.5}, {alpha, beta}] // Flatten

plot2 = 
      DisplayTogether[
ListPlot[{{0, 1}, {1, 0.8}, {2, 
              0.64}}, PlotStyle -> 
            PointSize[0.025]], 
Plot[
          Exp[alpha*x + beta*x^2] /. case1, {x, 0, 2}, 
          PlotStyle -> 
            Thickness[0.01]], 
ListPlot[{{0, 1}, {1, 
              0.8}, {2, 0.78}}, 
          PlotStyle -> PointSize[0.025]], 
Plot[
          Exp[alpha*x + beta*x^2] /. case2, {x, 0, 2}, 
          PlotStyle -> 
            Dashing[{0.01, 0.01}]], 
ListPlot[{{0, 1}, {1, 
              0.8}, {2, 0.5}}, 
          PlotStyle -> PointSize[0.025]], 
Plot[
          Exp[alpha*x + beta*x^2] /. case3, {x, 0, 
            2}], 
PlotRange -> {{-0.1, 2.1}, {0.0, 
              1.1}}
]

Display["graph1.eps", plot1, "EPS"]

Display["graph2.eps", plot2, "EPS"]