Generating the Figures in Chapter 5 of Otto and Day

Inputs (Enter first)

<< Graphics`Graphics`

<< Graphics`MultipleListPlot`

<< Graphics`Arrow`

Off[General::"spell1"]

Off[General::"spell"]

Needs["Graphics`PlotField`"]

Figure 5.1 - Solving an implicit equation for a variable of interest

Plot[n - (1 + r) Exp[-a n] n  -  m /. r -> 1 /. 
          m -> 3 /. a -> Log[200/97]/100, {n, -5, 
        200}]

Plot[n - (1 + r) Exp[-a n] n  -  m /. r -> 1 /. 
          m -> 3 /. a -> Log[200/97]/100, {n, 98, 
        102}]

First, we can try to solve the equation analytically using Solve.  Mathematica fails to find a solution.  One hint, when replacing parameters in an equation, place all of the terms and replacements within parentheses.

Solve[(n - (1 + r) Exp[-a n] n  -  m /. 
                r -> 1 /. m -> 3 /. 
            a -> Log[200/97]/100) == 0, n]

Without the parentheses, Mathematica is confused by what exactly should equal 0 and interprets the equation as claiming that Log[200/97]/100==0:

Solve[n - (1 + r) Exp[-a n] n  -  m /. r -> 1 /. 
          m -> 3 /. a -> Log[200/97]/100 == 0, 
      n]

FindRoot numerically finds a solution to an equation.  There may, however, be more than one solution, and Mathematica uses starting conditions to identify solutions in a particular region.  Here, starting at n=0 and then trying n=10, Mathematica identifies a root at -2.8788:

FindRoot[(n - (1 + r) Exp[-a n] n  -  m /. 
                r -> 1 /. m -> 3 /. 
            a -> Log[200/97]/100) == 0, {n, 0, 10}]

Starting at n=0 and then trying n=1000, Mathematica identifies a root at 100:

FindRoot[(n - (1 + r) Exp[-a n] n  -  m /. 
                r -> 1 /. m -> 3 /. 
            a -> Log[200/97]/100) == 0, {n, 0, 1000}]

Figure 5.2 and 5.4 - Stability analysis

DisplayTogether[
Plot[
        Exp[-0.1 t], {t, 0, 
          50}], 
Plot[-Exp[-0.1 t], {t, 0, 
          50}], 
Ticks -> 
        False
]

DisplayTogether[
Plot[
        Exp[0.1 t], {t, 0, 
          50}], 
Plot[-Exp[0.1 t], {t, 0, 
          50}], 
Ticks -> 
        False
]

DisplayTogether[
Plot[
        Exp[-0.1 t], {t, 0, 
          50}], 
Plot[-Exp[-0.1 t], {t, 0, 
          50}], 
Ticks -> 
        False
]

Clear[num, lam]
\(num[lam_, t_] := \(num[lam, t] = 
		lam\ num[lam, t - 1]\)\)
num[lam_, 0] :=  1.0

ListPlot[Table[{t, num[-1.1, t]}, {t, 0, 20}], 
      PlotJoined -> True, 
Ticks -> False]

ListPlot[Table[{t, num[-0.9, t]}, {t, 0, 20}], 
      PlotJoined -> True, 
Ticks -> False, 
      PlotRange -> All]

Figure 5.3 - Diploid Selection (p[t+1] vs p[t])

Simplify[
      D[(p^2 WAA  +  p q WAa)/(p^2 WAA  +  2 p q WAa +  
                  q^2 Waa) /. q -> (1 - p), p] /. 
        p -> (WAa - Waa)/(2 WAa - WAA - Waa)]

\(DisplayTogether[
	
      Plot[\((p^2\ WAA\  + \ p\ q\ WAa)\)/\((p^2\ WAA\  + \ 2\ p\ q\ WAa + \ 
                q^2\ Waa)\) /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, \(-0.1\), 1.1}], 
	
      Plot[p\ , {p, \(-0.1\), 1.1}, 
        PlotStyle -> Dashing[{0.02}]]
]\)

{(WAa/Waa), (WAA (WAa - Waa)  + 
              Waa (WAa - WAA)) /(WAa^2 - WAA Waa) , 
        WAa/WAA} /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
        q -> (1 - p)}

\(DisplayTogether[
	
      Plot[WAa\/Waa*p /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, \(-0.01\), 0.01}], 
	
      Plot[p\ , {p, \(-0.01\), 0.01}, PlotStyle -> Dashing[{0.02}]], 
      Ticks -> {{}, {}}
]\)

\(DisplayTogether[
	
      Plot[\((WAa - Waa)\)/\((2\ WAa - WAA - 
                  Waa)\) + \((WAA \((WAa - Waa)\)\  + 
                    Waa \((WAa - WAA)\))\)\ /\((WAa^2 - 
                    WAA\ Waa)\)\ \((p - \((WAa - Waa)\)/\((2\ WAa - Waa - 
                        WAA)\))\) /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, 2/3 - 0.01, 2/3 + 0.01}], 
	
      Plot[p\ , {p, 2/3 - 0.01, 2/3 + 0.01}, PlotStyle -> Dashing[{0.02}]], 
      AxesOrigin -> {2/3, 2/3}, Ticks -> {{}, {}}
]\)

\(DisplayTogether[
	
      Plot[1 - WAa\/WAA*\((1 - p)\) /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, 0.99, 1.01}], 
	
      Plot[p\ , {p, 0.99, 1.01}, PlotStyle -> Dashing[{0.02}]], 
      AxesOrigin -> {1, 1}, Ticks -> {{}, {}}
]\)

Figure 5.5 - Diploid Selection (dp/dt vs p)

The differential equation for the diploid model of selection is (see Problem 3.20):

(1 - 
          p) p (p (WAA - WAa)  + (1 - p) (WAa - 
                Waa))

\(DisplayTogether[
	
      Plot[\((1 - 
                p)\)\ p\ \((p \((WAA - WAa)\)\  + \((1 - p)\) \((WAa - 
                      Waa)\))\) /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, \(-0.1\), 
          1.1}], 
PlotRange -> {{\(-0.1\), 
            1.1}, {\(-0.12\), 0.12}}
]\)

Solve[ (p (WAA - WAa)  + (1 - p) (WAa - 
                  Waa)) == 0, p]

D[(1 - 
                p) p (p (WAA - WAa)  + (1 - p) (WAa - 
                      Waa)), p] /. 
        p -> ((WAa - Waa)/(2 WAa - Waa - WAA)) // Factor

{(WAa - 
            Waa), (-(((WAa - Waa) (WAa - WAA))/(2 WAa - 
                  Waa - WAA))), (WAa - WAA)} /. {WAA -> 0.6, 
        WAa -> 0.2, Waa -> 1}

\(DisplayTogether[
	
      Plot[\((WAa - Waa)\)*p /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, \(-0.01\), 0.01}], 
      PlotRange -> {{\(-0.01\), 0.01}, {\(-0.01\), 0.01}}, 

      Ticks -> {{}, {}}
]\)

\(DisplayTogether[
	
      Plot[\((\(-\(\(\((WAa - Waa)\)\ \((WAa - WAA)\)\)\/\(2\ WAa - Waa - 
                      WAA\)\)\))\)\ \((p - \((WAa - Waa)\)/\((2\ WAa - Waa - 
                      WAA)\))\) /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, 2/3 - 0.01, 2/3 + 0.01}], 

      PlotRange -> {{2/3 - 0.01, 2/3 + 0.01}, {\(-0.01\), 0.01}}, 

      AxesOrigin -> {2/3, 0}, Ticks -> {{}, {}}
]\)

\(DisplayTogether[
	
      Plot[\((WAa - WAA)\)*\((p - 1)\) /. {WAA -> 0.6, WAa -> 0.2, Waa -> 1, 
            q -> \((1 - p)\)}, {p, 0.99, 1.01}], 

      PlotRange -> {{1 - 0.01, 1 + 0.01}, {\(-0.01\), 0.01}}, 

      AxesOrigin -> {1, 0}, Ticks -> {{}, {}}
]\)

Figure 5.6 - Pathological behaviour of  Logistic Growth (Discrete Time)

Clear[num, r, K]
\(num[K_, r_, n0_, 
        t_] := \(num[K, r, n0, t] = 
		num[K, r, n0, t - 1] + 
          r\ num[K, r, n0, 
              t - 1]\ \((1 - num[K, r, n0, t - 1]\ /K)\)\)\)
num[K_, r_, n0_, 0] :=  n0

\(DisplayTogether[
	
      ListPlot[Table[{t, num[1000, \(-0.1\), 1025, t]}, {t, 0, 50}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, num[1000, \(-0.1\), 975, t]}, {t, 0, 50}], 
        PlotJoined -> True], 

      PlotRange -> {{0, 45}, {0, 2500}}
]\)

DisplayTogether[
Plot[(n + 
                  r n (1 - n /K))/n /. K -> 1000 /. 
          r -> -0.1, {n, 0, 
          1500}], 
Plot[(n + r n (1 - n /K))/
              n /. K -> 1000 /. r -> 0.1, {n, 0, 
          1500}]
]

Figure 5.8 - Meme model

\(DisplayTogether[
	
      Plot[f, {f, 0, 1}, PlotStyle -> Dashing[{0.02}]], 
	
      Plot[\(f\ \((1 - δ)\) + α\ f^2\ \((1 - 
                    f)\) /. δ -> 0.4\) /. α\  -> 
            2, {f, 0, 1}]
]\)