Generating the Figures in Chapter 4 of Otto and Day

Inputs (Enter first)

<< Graphics`Graphics`

<< Graphics`MultipleListPlot`

<< Graphics`Arrow`

Off[General::"spell1"]

Needs["Graphics`PlotField`"]

<< Graphics`ImplicitPlot`

Figure 4.1 - Exponential Growth

Figure 4.1a

\(DisplayTogether[
	
      ListPlot[Table[R^t\ n0 /. {R -> 2, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      ListPlot[Table[R^t\ n0 /. {R -> 1.1, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      ListPlot[Table[R^t\ n0 /. {R -> 1.05, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      ListPlot[Table[R^t\ n0 /. {R -> 1.01, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      PlotRange -> {{0, 50.5}, {0, 10000}}
]\)

Export["Fig.eps", %, "EPS"];

\(DisplayTogether[
	
      ListPlot[Table[R^t\ n0 /. {R -> 0.5, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      ListPlot[Table[R^t\ n0 /. {R -> 0.9, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      ListPlot[Table[R^t\ n0 /. {R -> 0.99, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      ListPlot[Table[R^t\ n0 /. {R -> 0.999, \ n0 -> 1000}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.01]], 
	
      PlotRange -> {{0, 50.5}, {0, 1010}}
]\)

Figure 4.2 - Logistic Growth (Discrete Time)

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

Figure 4.2a

\(DisplayTogether[
	
      ListPlot[Table[{t, num[1000, 0.4, t]}, {t, 0, 50}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, num[1000, 0.7, t]}, {t, 0, 50}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, num[1000, 1.8, t]}, {t, 0, 50}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, num[1000, 2.1, t]}, {t, 0, 50}], 
        PlotJoined -> True], PlotRange -> {{0, 25}, {0, 1400}}
]\)

\(DisplayTogether[
	
      ListPlot[Table[{t, num[1000, 0.1, t]}, {t, 0, 50}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, num[1000, 1, t]}, {t, 0, 50}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, num[1000, 2.1, t]}, {t, 0, 50}], 
        PlotJoined -> True], PlotRange -> {{0, 50}, {0, 1400}}
]\)

Figure 4.2b

\(DisplayTogether[
	
      ListPlot[Table[num[1000, 2.5, t], {t, 0, 100}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 25}, {0, 1400}}
]\)

Figure 4.2c

\(DisplayTogether[
	
      ListPlot[Table[num[1000, 2.7, t], {t, 0, 100}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 100}, {0, 1400}}
]\)

Figure 4.2d

\(DisplayTogether[
	
      ListPlot[Table[num[1000, 3.01, t], {t, 0, 20}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 25}, {0, 1400}}
]\)

Figure 4.3 - Logistic Growth (Continuous Time)

Clear[DifLogistic, K, r, n0, n]
DifLogistic[K_, r_, n0_, tmin_, 
        tmax_] := DifLogistic[K, r, n0, tmin, tmax] = 
        NDSolve[{D[n[t], t] == r n[t] (1 - n[t] /K), n[tmin] == n0}, 
          n[t], {t, tmin, tmax}]

\(DisplayTogether[
	
      Plot[n[t] /. DifLogistic[1000, 1, 10, 0, 10], {t, 0, 10}], 
	
      Plot[n[t] /. DifLogistic[1000, 2.1, 10, 0, 10], {t, 0, 10}], 
	
      Plot[n[t] /. DifLogistic[1000, 2.7, 10, 0, 10], {t, 0, 10}], 
	
      Plot[n[t] /. DifLogistic[1000, 4, 10, 0, 10], {t, 0, 10}], 
\ \ \ 
      PlotRange -> {{0, 10}, {0, 1400}}
]\)

Figure 4.4 - Logistic Growth (Continuous vs Discrete Time)

\(DisplayTogether[
	
      ListPlot[Table[{t, num[1000, 0.05, t]}, {t, 0, 500}], 
        PlotStyle -> PointSize[0.007]], 
	
      ListPlot[Table[{t, num[1000, 0.1, t]}, {t, 0, 500}], 
        PlotStyle -> PointSize[0.007]], 
	
      Plot[n[t] /. DifLogistic[1000, 0.05, 10, 0, 500], {t, 0, 500}], 
	
      Plot[n[t] /. DifLogistic[1000, 0.1, 10, 0, 500], {t, 0, 
          500}], 
\ \ \ PlotRange -> {{0, 200}, {0, 1400}}
]\)

Figure 4.5 - Yeast consumer-resource model

LogGrowthC41a := {{0, 3.1*10^(5)}, {2, 4.8*10^(5)}, {4, 
            1.8*10^(6)}, {5, 3.6*10^(6)}, {8, 1.1*10^(7)}, {10, 
            3.2*10^(7)}, {12, 9.8*10^(7)}, {14, 3.2*10^(8)}, {16, 
            3.2*10^(8)}, {18, 3.0*10^(8)}, {20, 3.1*10^(8)}, {22, 
            3.0*10^(8)}, {24, 2.9*10^(8)}, {26, 3.5*10^(8)}, {28, 
            3.5*10^(8)}, {30, 3.5*10^(8)}, {32, 3.4*10^(8)}, {34, 
            3.8*10^(8)}, {36, 3.3*10^(8)}, {48, 3.9*10^(8)}, {72, 
            4.1*10^(8)}};

LogGrowthC41b := {{0, 3.0*10^(5)}, {2, 9.8*10^(5)}, {4, 
            1.4*10^(6)}, {5, 4.1*10^(6)}, {8, 1.6*10^(7)}, {10, 
            3.7*10^(7)}, {12, 8.2*10^(7)}, {14, 3.2*10^(8)}, {16, 
            3.0*10^(8)}, {18, 3.2*10^(8)}, {20, 3.5*10^(8)}, {22, 
            3.5*10^(8)}, {24, 3.6*10^(8)}, {26, 3.5*10^(8)}, {28, 
            3.4*10^(8)}, {30, 3.4*10^(8)}, {32, 3.3*10^(8)}, {34, 
            3.5*10^(8)}, {36, 3.6*10^(8)}, {48, 3.5*10^(8)}, {72, 
            3.9*10^(8)}};

logistic[K_,r_,N0_,t_]=K/(1+(K/N0-1)*Exp[-r*t]);

r41 = 0.545;
r42 = 0.548;

haploid0 = 305000;
diploid0 = 190000;

K41 = 3.7*10^8;

Clear[DifConsRes, theta, a, e, d, n10, n20, tmin, tmax]
DifConsRes[theta_, a_, e_, d_, n10_, n20_, tmin_, 
        tmax_] := DifConsRes[theta, a, e, d, n10, n20, tmin, tmax] = 
        NDSolve[{D[n1[t], t] == theta  -  a n1[t] n2[t],  
            D[n2[t], t] == e a n1[t] n2[t] -  d n2[t], 
            n1[tmin] == n10, n2[tmin] == n20}, {n1[t], n2[t]}, {t, tmin, 
            tmax}]

Plot1 = 
      MultipleListPlot[LogGrowthC41a, LogGrowthC41b, 
        PlotRange -> {{0, 73}, {-10^7, 4.2*10^8}}, 
        SymbolShape -> {PlotSymbol[Box], 
            MakeSymbol[{Line[{{3, 3}, {-3, -3}}], 
                Line[{{-3, 3}, {3, -3}}]}]}]

aequiv = r41/K41

(K41 - haploid0)

DisplayTogether[
Plot1,
Plot[n2[t]/.DifConsRes[0,aequiv,e,0,(K41-haploid0)/2,haploid0,0,100],{t,0,100},
  PlotStyle->{Thickness[0.01],Dashing[{0.02,0.04}]}],
Plot[logistic[K41,r41,haploid0,t],{t,0,100}]
]

Figure 4.6 - Haploid Selection (Continuous vs Discrete Time)

Clear[hap, s]
\(hap[s_, 
        t_] := \(hap[s, 
          t] = 
		hap[s, 
            t - 1] \((1 + 
                s)\)/\((hap[s, t - 1] \((1 + s)\) + \((1 - 
                    hap[s, t - 1])\))\)\)\)
hap[s_, 0] := 0.2

Clear[difhap, s, p0, p]
difhap[s_, p0_, tmin_, tmax_] := difhap[s, p0, tmin, tmax] = 
        NDSolve[{D[p[t], t] == s p[t] (1 - p[t] ), p[tmin] == p0}, 
          p[t], {t, tmin, tmax}]

\(DisplayTogether[
	
      ListPlot[Table[{t, hap[0.1, t]}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.007]], 
	
      Plot[p[t] /. difhap[0.1, 0.2, 0, 500], {t, 0, 100}], 
	
      ListPlot[Table[{t, hap[0.5, t]}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.007]], 
	
      Plot[p[t] /. difhap[0.5, 0.2, 0, 500], {t, 0, 100}], 
	
      ListPlot[Table[{t, hap[0.01, t]}, {t, 0, 100}], 
        PlotStyle -> PointSize[0.007]], 
	
      Plot[p[t] /. difhap[0.01, 0.2, 0, 500], {t, 0, 
          100}], 
\ \ \ PlotRange -> {{0, 100}, {0, 1.01}}
]\)

Figure 4.7 - Diploid Selection (Discrete Time)

Clear[dip, s]
\(dip[WAA_, WAa_, Waa_, dip0_, 
        t_] := \(dip[WAA, WAa, Waa, dip0, 
          t] = 
		\((dip[WAA, WAa, Waa, dip0, t - 1]^2  WAA + 
              dip[WAA, WAa, Waa, dip0, 
                  t - 1] \((1 - dip[WAA, WAa, Waa, dip0, t - 1])\) 
                WAa)\)/\((dip[WAA, WAa, Waa, dip0, t - 1]^2  WAA + 
              2  dip[WAA, WAa, Waa, dip0, 
                  t - 1] \((1 - 
                    dip[WAA, WAa, Waa, dip0, t - 1])\)\ WAa + \((1 - 
                      dip[WAA, WAa, Waa, dip0, t - 1])\)^2  
                Waa)\)\)\)
dip[WAA_, WAa_, Waa_, dip0_, 0] := dip0

Let's set the fitnesses to different orderings of 0.8, 0.95, and 1.  First consider the case where the heterozygote is most fit:

\(DisplayTogether[
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.99, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.9, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.7, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.5, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.3, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.1, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 1, 0.95, 0.01, t]}, {t, 0, 120}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 120}, {0, 1}}
]\)

\(DisplayTogether[
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8\ , 0.99, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8, 0.9, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8\ , 0.7, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8\ , 0.5, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8\ , 0.3, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8, 0.1, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[1, 0.95, 0.8\ , 0.01, t]}, {t, 0, 120}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 120}, {0, 1}}
]\)

\(DisplayTogether[
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1\ , 0.99, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1, 0.9, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1\ , 0.7, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1\ , 0.5, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1\ , 0.3, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1\ , 0.1, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.8, 0.95, 1\ , 0.01, t]}, {t, 0, 120}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 120}, {0, 1}}
]\)

\(DisplayTogether[
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.99, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.9, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.7, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.58, t]}, {t, 0, 120}], 
        PlotJoined -> True, PlotStyle -> Dashing[{0.02}]], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.56, t]}, {t, 0, 120}], 
        PlotJoined -> True, PlotStyle -> Dashing[{0.02}]], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.5, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.3, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.1, t]}, {t, 0, 120}], 
        PlotJoined -> True], 
	
      ListPlot[Table[{t, dip[0.95, 0.8, 1, 0.01, t]}, {t, 0, 120}], 
        PlotJoined -> 
          True], 
\ \ \ PlotRange -> {{0, 100}, {0, 1}}
]\)

Figure 4.8 - Lotka-Volterra Competition Model (Discrete Time)

Clear[pop1,pop2]
pop1[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,t_]:=pop1[r1,r2,K1,K2,a12,a21,n10,n20,t]=
	pop1[r1,r2,K1,K2,a12,a21,n10,n20,t-1]*(1+r1*
	(1-(pop1[r1,r2,K1,K2,a12,a21,n10,n20,t-1]+a12*pop2[r1,r2,K1,K2,a12,a21,n10,n20,t-1])/K1))
pop2[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,t_]:=pop2[r1,r2,K1,K2,a12,a21,n10,n20,t]=
	pop2[r1,r2,K1,K2,a12,a21,n10,n20,t-1]*(1+r2*
	(1-(pop2[r1,r2,K1,K2,a12,a21,n10,n20,t-1]+a21*pop1[r1,r2,K1,K2,a12,a21,n10,n20,t-1])/K2))
pop1[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,0]:=n10
pop2[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,0]:=n20

DisplayTogether[
	ListPlot[Table[{t,pop1[0.5,0.5,1000,1000,0,0.5,10,10,t]},{t,0,50}],
		PlotJoined->True,PlotRange->{Automatic,{0,1500}}],
	ListPlot[Table[{t,pop2[0.5,0.5,1000,1000,0,0.5,10,10,t]},{t,0,50}],
		PlotJoined->True,PlotRange->{Automatic,{0,1500}},PlotStyle->Dashing[{0.02}]]
		]

This gives the change in species 1 and the change in species 2 as two surfaces over the x-y plane representing the number of species 1 and the number of species 2.  Only when both surfaces cross the x-y plane do both species remain constant in size.  Such 3-D plots are, however, hard to visualize and interpret.

\(DisplayTogether[
	
      Plot3D[\(\(n1*\((r1*\((1 - \((n1 + a12*n2)\)/K1)\))\) /. 
              K1 -> 1000\) /. r1 -> 0.5\) /. a12 -> 0, {n1, 0, 1500}, {n2, 0, 
          1500}], 
	
      Plot3D[\(\(n2*\((r2*\((1 - \((n2 + a21*n1)\)/K2)\))\) /. 
              K2 -> 1000\) /. r2 -> 0.5\) /. a21 -> 0.5, {n1, 0, 1500}, {n2, 
          0, 1500}], 
	
      PlotRange -> {Automatic, Automatic, {\(-100\), 100}}
]\)

Figure 4.9 - Haploid Selection (p[t+1] vs p[t])

Figure 4.8a

\(DisplayTogether[
	
      Plot[p\ WA\ /\((p\ WA\  + \ q\ Wa)\) /. 
		{WA -> 2, Wa -> 1, 
            q -> \((1 - p)\)}, {p, 0, 1}], 
	
      Plot[p\ , {p, 0, 1}, PlotStyle -> Dashing[{0.02}]]
]\)

Figure 4.8b

\(DisplayTogether[
	
      Plot[p\ WA\ /\((p\ WA\  + \ q\ Wa)\) /. 
		{WA -> 0.5, Wa -> 1, 
            q -> \((1 - p)\)}, {p, 0, 1}], 
	
      Plot[p\ , {p, 0, 1}, PlotStyle -> Dashing[{0.02}]]
]\)

WAvalue = 1;
Wavalue = 0.5;
initp = 0.5;

DisplayTogether[
      Plot[p WAvalue/(p WAvalue + (1 - p) Wavalue), {p, 0, 1}], 
      Plot[p, {p, 0, 1}, PlotStyle -> Dashing[{0.02}]], 
      ListPlot[pvalue = initp; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, 0}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, 
        PlotJoined -> True]]

WAvalue = 0.5;
Wavalue = 1.0;
initp = 0.5;

DisplayTogether[
      Plot[p WAvalue/(p WAvalue + (1 - p) Wavalue), {p, 0, 1}], 
      Plot[p, {p, 0, 1}, PlotStyle -> Dashing[{0.02}]], 
      ListPlot[pvalue = initp; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, 0}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, PlotJoined -> True], 
      ListPlot[pvalue = newpvalue; 
        newpvalue = 
          pvalue WAvalue/(pvalue WAvalue + (1 - 
                        pvalue) Wavalue); {{pvalue, pvalue}, {pvalue, 
            newpvalue}, {newpvalue, newpvalue}}, 
        PlotJoined -> True]]

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

\(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}], 
	
      Plot[p\ , {p, 0, 1}, PlotStyle -> Dashing[{0.02}]]
]\)

Figure 4.11 - Haploid Selection (Change in p vs p[t])

Figure 4.10a

Plot[s  p q /. {s -> 0.01, q -> (1 - p)}, {p, 0, 1}]

Figure 4.10b

Plot[s  p q /. {s -> -0.01, q -> (1 - p)}, {p, 0, 1}]

Figure 4.12 - Logistic Growth (Change in n vs n[t])

\(DisplayTogether[
	
      Plot[r\ n\ \((1 - n/K)\) /. {r -> 2.5, K -> 1000}, {n, 0, 1200}], 
      Plot[r\ n\ \((1 - n/K)\) /. {r -> 1, K -> 1000}, {n, 0, 1200}]
]\)

Figure 4.13 - Phase-line Diagram

differenceeqn1 = s*p1*(1 - p1);
differenceeqn2 = 0;

DisplayTogether[
  Plot[0, {p1, 0, 1}],
\(PlotVectorField[{differenceeqn1, 
              differenceeqn2} /. s -> 0.1, {p1, 0.05, 0.95}, {p2, 0.0, 
            0.02}, \ 
		PlotPoints\  -> {20, 1}, 
          ScaleFunction\  -> \ \((0.01 &)\)]\),
PlotRange -> {{0, 1}, {-0.01, 0.01}},
AspectRatio->0.1,Ticks -> {Automatic, None}
]

Display["testA.eps", %, "EPS"]

DisplayTogether[
  Plot[0, {p1, 0, 1}],
\(PlotVectorField[{differenceeqn1, 
              differenceeqn2} /. s -> \(-0.1\), {p1, 0.05, 0.95}, {p2, 
            0.0, 0.02}, \ 
		PlotPoints\  -> {20, 1}, 
          ScaleFunction\  -> \ \((0.01 &)\)]\),
PlotRange -> {{0, 1}, {-0.01, 0.01}},
AspectRatio->0.1,Ticks -> {Automatic, None}
]

Display["testB.eps", %, "EPS"]

differenceeqn1 = -p1 + (((1 - p1) p1 WAa + (p1^2) WAA)/((((1 - p1))^2) Waa + 2 (1 - p1) p1 WAa + (p1^2) WAA));
differenceeqn2 = 0;

DisplayTogether[
  Plot[0, {p1, 0, 1}],
\(PlotVectorField[\(\({differenceeqn1, 
                  differenceeqn2} /. WAA -> 0.8\) /. WAa -> 1\) /. 
            Waa -> 0.95, {p1, 0.05, 0.95}, {p2, 0.0, 
            0.02}, \ 
		PlotPoints\  -> {20, 1}, 
          ScaleFunction\  -> \ \((0.01 &)\)]\),
PlotRange -> {{0, 1}, {-0.01, 0.01}},
AspectRatio->0.1,Ticks -> {Automatic, None}
]

Display["test1.eps", %, "EPS"]

DisplayTogether[
  Plot[0, {p1, 0, 1}],
\(PlotVectorField[\(\({differenceeqn1, 
                  differenceeqn2} /. WAA -> 1\) /. WAa -> 0.95\) /. 
            Waa -> 0.8, {p1, 0.05, 0.95}, {p2, 0.0, 
            0.02}, \ 
		PlotPoints\  -> {20, 1}, 
          ScaleFunction\  -> \ \((0.01 &)\)]\),
PlotRange -> {{0, 1}, {-0.01, 0.01}},
AspectRatio->0.1,Ticks -> {Automatic, None}
]

Display["test2.eps", %, "EPS"]

DisplayTogether[
  Plot[0, {p1, 0, 1}],
\(PlotVectorField[\(\({differenceeqn1, 
                  differenceeqn2} /. WAA -> 1\) /. WAa -> 0.8\) /. 
            Waa -> 0.95, {p1, 0.05, 0.95}, {p2, 0.0, 
            0.02}, \ 
		PlotPoints\  -> {20, 1}, 
          ScaleFunction\  -> \ \((0.01 &)\)]\),
PlotRange -> {{0, 1}, {-0.01, 0.01}},
AspectRatio->0.1,Ticks -> {Automatic, None}
]

Display["test3.eps", %, "EPS"]

Figure 4.14 - Phase-plane Diagram of the Lotka-Volterra Competition Model (Discrete Time)

Clear[pop1,pop2]
pop1[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,t_]:=pop1[r1,r2,K1,K2,a12,a21,n10,n20,t]=
	pop1[r1,r2,K1,K2,a12,a21,n10,n20,t-1]*(1+r1*
	(1-(pop1[r1,r2,K1,K2,a12,a21,n10,n20,t-1]+a12*pop2[r1,r2,K1,K2,a12,a21,n10,n20,t-1])/K1))
pop2[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,t_]:=pop2[r1,r2,K1,K2,a12,a21,n10,n20,t]=
	pop2[r1,r2,K1,K2,a12,a21,n10,n20,t-1]*(1+r2*
	(1-(pop2[r1,r2,K1,K2,a12,a21,n10,n20,t-1]+a21*pop1[r1,r2,K1,K2,a12,a21,n10,n20,t-1])/K2))
pop1[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,0]:=n10
pop2[r1_,r2_,K1_,K2_,a12_,a21_,n10_,n20_,0]:=n20

tab1 = 
        Table[{pop1[0.5, 0.5, 1000, 1000, 0, 0.5, 10, 10, t], 
            pop2[0.5, 0.5, 1000, 1000, 0, 0.5, 10, 10, t]}, {t, 0, 
            50}];
tab2 = 
        Table[{pop1[0.5, 0.5, 1000, 1000, 0.5, 0.5, 10, 10, t], 
            pop2[0.5, 0.5, 1000, 1000, 0.5, 0.5, 10, 10, t]}, {t, 0, 
            50}];
tab3 = 
        Table[{pop1[0.5, 0.5, 1000, 1000, 1, 0.5, 10, 10, t], 
            pop2[0.5, 0.5, 1000, 1000, 1, 0.5, 10, 10, t]}, {t, 0, 
            50}];
tab4 = 
        Table[{pop1[0.5, 0.5, 1000, 1000, 2, 0.5, 10, 10, t], 
            pop2[0.5, 0.5, 1000, 1000, 2, 0.5, 10, 10, t]}, {t, 0, 
            50}];

MultipleListPlot[tab1,tab2,tab3,tab4,PlotRange->{{0,1250},{0,1250}},
  SymbolShape -> {PlotSymbol[Triangle], PlotSymbol[Diamond,Filled->False], PlotSymbol[Star], MakeSymbol[
      {Line[{{2, 2}, {-2, -2}}],
       Line[{{-2, 2}, {2, -2}}]}
     ]},
    Epilog -> {Arrow[{130,40}, {280,190}]}
     ]

Figure 4.15 - Vector-Field Diagram of the Lotka-Volterra Competition Model (Discrete Time)

differenceeqn1 = n1*r1*(1 - (n1 + a12*n2)/K1);
differenceeqn2 = n2*r2*(1 - (n2 + a21*n1)/K2);

subparameters  =  {r1 -> 0.5, K1 -> 1000, r2 -> 0.5, K2 -> 1000, 
          a21 -> 0.5, a12 -> 0};

DisplayTogether[
	ListPlot[{{1000, 500}}, PlotStyle -> PointSize[0.03]],
\(PlotVectorField[{differenceeqn1, 
              differenceeqn2} /. subparameters, {n1, 50, 1200}, {n2, 50, 
            1200}, \ 
		PlotPoints\  -> 10, 
          ScaleFunction\  -> \ \((1  # &)\), ScaleFactor\  -> None]\),
PlotRange -> {{0, 1250}, {0, 1250}},
AspectRatio->1/GoldenRatio
]

1200 + 0.5 1200 (1 - (1200 + 0.5 50)/1000)

n1 + differenceeqn1 /. subparameters /. n1 -> 50

n2 + differenceeqn2 /. subparameters /. n1 -> 50 /. 
      n2 -> 1200

Figure 4.16 - Null-clines of the Lotka-Volterra Competition Model (Discrete Time)
(needs definitions from Figure 4.12)

differenceeqn1 = n1*r1*(1 - (n1 + a12*n2)/K1);
differenceeqn2 = n2*r2*(1 - (n2 + a21*n1)/K2);

subparameters  =  {r1 -> 0.5, K1 -> 1000, r2 -> 0.5, K2 -> 1000, 
          a21 -> 0.5, a12 -> 0};

Solve[differenceeqn1 == 0, n1]

Solve[differenceeqn1 == 0, n2]

Solve[differenceeqn2 == 0, n1]

Solve[differenceeqn2 == 0, n2]

In this example with a12=0, the null-cline for n1 = K1.  This corresponds to a spike at 1000, which can be drawn by connecting two dots, one at {1000,0} and one at {1000,1100} using ListPlot.

\(DisplayTogether[
	
      Plot[\(K2 - a21\ n1 /. K2 -> 1000\) /. a21 -> 0.5\ , {n1, 0, 
          1250}], 
	
      ListPlot[{{1000, 0}, {1000, 1250}}, PlotJoined -> True, 
        PlotStyle -> Dashing[{0.01, 0.01}]], 
	
      MultipleListPlot[tab1, 
        SymbolShape\  -> \ {PlotSymbol[
              Triangle]}], 
\ \ \ \ \ \ \ ListPlot[{{1000, 
            500}}, PlotStyle -> PointSize[0.03]], 
	
      PlotRange -> {{0, 1250}, {0, 
            1250}}, 
\ \ \ \ Epilog\  -> \ {Arrow[{130, 40}, \ {280, 
              190}]}
]\)

More generally, we can use ContourPlot:

\(DisplayTogether[

      ContourPlot[
        differenceeqn1 /. subparameters, {n1, 1, 1250}, {n2, 1, 
          1250}, 
		Contours -> {0}, ContourShading -> False, 
        ContourStyle -> Dashing[{0.01}]], 

      ContourPlot[
        differenceeqn2 /. subparameters, {n1, 1, 1250}, {n2, 1, 
          1250}, 
		Contours -> {0}, ContourShading -> False], 

      MultipleListPlot[
        tab1, 
\ \ SymbolShape\  -> \ {PlotSymbol[Triangle]}\ \ \ \ \ ], 
	
      PlotRange -> {{0, 1200}, {0, 
            1200}}, 
\ \ \ \ Epilog\  -> \ {Arrow[{130, 40}, \ {280, 
              190}]}, 
AspectRatio -> 
        1/GoldenRatio
]\)

Figure 4.17 - Phase-plane Diagram of the Lotka-Volterra Predator-Prey Model (Continuous Time)

Here, we plot the predator-prey model in continuous time using the differential equations:

Prey:  		dn1/dt = r n1[t] - c n1[t] n2[t]
Predator:	dn2/dt =  e c n1[t] n2[t] - d  n2[t]

where we define c to equal the product of the contact rate and the probability that the predator successfully consumes the prey (= a c in the text).

In the following, we define a function, "system" that gives the numerical solution to these two differential equations for parameters of interest:

Clear[system,n1,n2]
system[r_,c_,e_,d_,n10_,n20_,tmin_,tmax_] := 
	system[r,c,e,d,n10,n20,tmin,tmax] =
	NDSolve[{n1'[t] == n1[t]*(r - c*n2[t]),
				n2'[t] == n2[t]*(e*c*n1[t] - d),
				n1[0]  == n10,
				n2[0]  == n20},
			{n1[t], n2[t]}, {t, tmin, tmax}]

DisplayTogether[
ParametricPlot[Evaluate[{n1[t], n2[t]} /. 
		system[0.5,0.01,0.02,0.1,500,45,0,100] ],
	{t, 0, 30}],
ParametricPlot[Evaluate[{n1[t], n2[t]} /. 
		system[0.5,0.01,0.02,0.1,500,40,0,100] ],
	{t, 0, 30}],
ParametricPlot[Evaluate[{n1[t], n2[t]} /. 
		system[0.5,0.01,0.02,0.1,500,35,0,100] ],
	{t, 0, 30}],
ParametricPlot[Evaluate[{n1[t], n2[t]} /. 
		system[0.5,0.01,0.02,0.1,500,30,0,100] ],
	{t, 0, 30}],
    PlotRange -> {{0, 1250}, {0, 100}}
]

differentialeqn1= n1*(r - c*n2);
differentialeqn2=n2*(e*c*n1 - d);

subparameters  =  {r -> 0.5, c -> 0.01, e -> 0.02, 
          d -> 0.1 };

The equilibrium requires that the number of species 1 remains constant as does the number of species 2.  These requirements are met only at the following two points:

DisplayTogether[
ParametricPlot[
        Evaluate[{n1[t], n2[t]} /. 
            system[0.5, 0.01, 0.02, 0.1, 500, 45, 0, 100]], {t, 0, 
          30}], 
ParametricPlot[
        Evaluate[{n1[t], n2[t]} /. 
            system[0.5, 0.01, 0.02, 0.1, 500, 40, 0, 100]], {t, 0, 30}], 
      ParametricPlot[
        Evaluate[{n1[t], n2[t]} /. 
            system[0.5, 0.01, 0.02, 0.1, 500, 35, 0, 100]], {t, 0, 30}], 
      ParametricPlot[
        Evaluate[{n1[t], n2[t]} /. 
            system[0.5, 0.01, 0.02, 0.1, 500, 30, 0, 100]], {t, 0, 30}], 
      Plot[r/c /. subparameters, {n1, 0, 1250}, 
        PlotStyle -> Dashing[{0.01, 0.01}]], 
      ListPlot[{{500, 0}, {500, 150}}, PlotJoined -> True], 
      ListPlot[{{500, 50}}, PlotStyle -> PointSize[0.03]], 
      PlotRange -> {{0, 1250}, {0, 100}}]

Export["Fig.eps", %, "EPS"];

Figure 4.18 - Vector-Field Diagram of the Lotka-Volterra Predator-Prey Model (Continuous Time)

differentialeqn1= n1*(r - c*n2);
differentialeqn2=n2*(e*c*n1 - d);

subparameters  =  {r -> 0.5, c -> 0.01, e -> 0.02, 
          d -> 0.1 };

The equilibrium requires that the number of species 1 remains constant as does the number of species 2.  These requirements are met when:

Solve[differentialeqn1  ==  0, {n1, n2}]

Solve[differentialeqn2  ==  0, {n1, n2}]

Solve[{differentialeqn1  ==  0, differentialeqn2 ==  0}, {n1, 
        n2}]

This equilibrium was used to help choose parameters that gave approximately 500 prey and 50 predators at equilibrium:

{n1 -> (d/(c e)), n2 -> (r/c)} /. subparameters

The null-cline for n2 occurs at r/c.  This corresponds to a line at 50.The null-cline for n1 occurs at d/(c e).  This corresponds to a spike at 500, which can be drawn by connecting two dots, one at {500,0} and one at {500,105} using ListPlot.

\(DisplayTogether[
	
      Plot[r/c /. subparameters\ , {n1, 0, 1250}, 
        PlotStyle -> Dashing[{0.01, 0.01}]], 
	
      ListPlot[{{500, 0}, {500, 150}}, PlotJoined -> True], 

      ListPlot[{{500, 50}}, PlotStyle -> PointSize[0.03]], 
      PlotVectorField[{differentialeqn1, differentialeqn2} /. 
          subparameters, {n1, 50, 1100}, {n2, 5, 90}, PlotPoints -> 15, 
        ScaleFunction -> None, 
        ScaleFactor -> 
          None], 
PlotRange -> {{0, 1200}, {0, 
            100}}, 
AspectRatio -> 
        1/GoldenRatio]\)

\(DisplayTogether[
	
      Plot[r/c /. subparameters\ , {n1, 0, 1250}, 
        PlotStyle -> Dashing[{0.01, 0.01}]], 
	
      ListPlot[{{500, 0}, {500, 150}}, PlotJoined -> True], 

      ListPlot[{{500, 50}}, PlotStyle -> PointSize[0.03]], 
      PlotVectorField[{differentialeqn1, differentialeqn2} /. 
          subparameters, {n1, 50, 1100}, {n2, 5, 90}, PlotPoints -> 15, 
        ScaleFunction -> None, 
        ScaleFactor -> 
          None], 
PlotRange -> {{0, 1200}, {0, 
            100}}, 
AspectRatio -> 
        1/GoldenRatio]\)

Alternatively, we can use ContourPlot:

DisplayTogether[
ContourPlot[(differentialeqn1 /. 
            subparameters), {n1, 0, 1100}, {n2, 0, 1100}, 
        Contours -> {0}, ContourShading -> False, 
        PlotPoints -> 100], 
      ContourPlot[(differentialeqn2 /. subparameters), {n1, 0, 
          1100}, {n2, 0, 1100}, Contours -> {0}, 
        ContourShading -> False, ContourStyle -> Dashing[{0.01}], 
        PlotPoints -> 100], 
ListPlot[{{500, 50}}, 
        PlotStyle -> PointSize[0.03]], 
      PlotVectorField[{differentialeqn1, differentialeqn2} /. 
          subparameters, {n1, 50, 1200}, {n2, 5, 95}, PlotPoints -> 15, 
        ScaleFunction -> None, 
        ScaleFactor -> 
          None], 
PlotRange -> {{0, 1250}, {0, 
            105}}, 
AspectRatio -> 
        1/GoldenRatio]

Figures 4.19 - Phillips' model

Change the value of beta here to explore its effect on the figure:

newbeta = 0.00027;

\(\(VarSet = {Γ -> 1.36, μ -> 
            1.36*10^\((\(-3\))\), τ -> 0.2, β -> 
            newbeta, 
		\ p -> 0.1, α -> 
            3.6*10^\((\(-2\))\), σ -> 2, δ -> 0.33, πP -> 
            100};\)\)

{β - ((α δ μ σ + δ (μ^2) σ)/(α Γ πP τ + (1 - p) Γ μ πP τ))} /. VarSet

Clear[Solution]
\(Solution[t_] = 

      NDSolve[{
	
            D[r[t], t] == Γ*τ - μ*r[t] - β*r[t]*
                  v[t], 
	
            D[l[t], t] == 
              p*β*r[t]*v[t] - μ*l[t] - α*l[t], 
	
            D[e[t], t] == \((1 - p)\)*β*r[t]*v[t] + α*
                  l[t] - δ*e[t], 
	
            D[v[t], t] == πP*e[t] - σ*v[t], 
	
	
            r[0] == 1000*τ, 
	l[0] == 0, 
	e[0] == 0, 
	
            v[0] == 4*10^\((\(-7\))\)}
 /. VarSet, 
{r[t], l[t], e[t], 
          v[t]}, {t, 0, 220}]\)

plot1 = 
      DisplayTogether[
LogPlot[(v[t] /. 
                Solution[t])[[1]], {t, 0, 120}, 
          PlotRange -> {Automatic, {0.1, 10000}}, 
          PlotStyle -> {Dashing[{0}]}], 
LogPlot[(e[
                  t] /. Solution[t])[[1]], {t, 0, 120}, 
          PlotStyle -> {Dashing[{ .01}]}], 
LogPlot[(l[t] /. Solution[t])[[1]], {t, 0, 120}, 
          PlotStyle -> {Dashing[{ .05}]}]]

Display["F1.EPS", %, "EPS"]

Problem 4.20 - Diploid model

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

Problem 4.21 - SI model

differenceeqn1 = 
        b - d n1  -  a c n1 n2  +  rho n2;
differenceeqn2 = 
        a c  n1 n2 - rho  n2  -  delta n2;

Solve[differenceeqn1 == 0, n1]

Solve[differenceeqn1 == 0, n2]

Solve[differenceeqn2 == 0, n1]

Solve[differenceeqn2 == 0, n2]

Solve[{differenceeqn1 == 0, differenceeqn2 == 0}, {n1, n2}]

equil = {((delta + rho)/(a c)), ((b - ((d (delta + rho))/(a c)))/delta)};

subparameters  =  {b -> 4,  d -> 0.01,  delta -> 0.02, 
          rho -> 0.2, a -> 0.1, c -> 0.01};

equil /. subparameters

Here the arrows are drawn with magnitude of 8:

\(DisplayTogether[

      PlotVectorField[{differenceeqn1, differenceeqn2} /. subparameters, {n1, 
          185, 275}, {n2, 5, 120}, \ PlotPoints\  -> 18, 

        ScaleFunction\  -> \ \((5 &)\), 
		\ \ ScaleFactor\  -> None], 
	
      PlotRange -> {{180, 280}, {0, 125}}, 
      Frame -> True, 
AspectRatio -> 
        1/GoldenRatio
]\)

Figure 4.1.1-  Sensitivity to Initial Conditions

Clear[num, r, K, n0]
\(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

Figure 4.2c

\(DisplayTogether[
	
      ListPlot[Table[num[20000, 2.7, 1, t], {t, 0, 500}], 
        PlotJoined -> True], 
	
      ListPlot[Table[num[20000, 2.7, 10000, t], {t, 0, 500}], 
        PlotJoined -> True, PlotStyle -> Dashing[{0.01}]], 
	
      ListPlot[Table[num[20000, 2.7, 10001, t], {t, 0, 500}], 
        PlotJoined -> True, 
        PlotStyle -> 
          Dashing[{0.04}]], 
\ \ \ PlotRange -> {{20, 40}, {5000, 
            Automatic}}, 
AxesOrigin -> {20, 
          5000}
]\)

\(DisplayTogether[
	
      ListPlot[Table[num[20000, 2.7, 1, t], {t, 0, 500}], 
        PlotJoined -> True], 
	
      ListPlot[Table[num[20000, 2.7, 10000, t], {t, 0, 500}], 
        PlotJoined -> True, PlotStyle -> Dashing[{0.01}]], 
	
      ListPlot[Table[num[20000, 2.7, 10001, t], {t, 0, 500}], 
        PlotJoined -> True, 
        PlotStyle -> 
          Dashing[{0.04}]], 
\ \ \ PlotRange -> {{485, 500}, {5000, 
            Automatic}}, 
AxesOrigin -> {485, 
          5000}
]\)

Figure 4.1.2 - Bifurcation diagram

The following makes a big long list ("thelist") of r values over the final "takelast" iterations of a simulation that runs for a specified number of steps ("timesteps").

Clear[bifurcationplot]
\(bifurcationplot[{rbegin_, rend_, rstep_}, timesteps_, 
          takelast_] := 
        Module[{n, \ iter, r, 
            thelist = {}}, 
n[K_, r_, 
              t_] := \(n[K, r, t] = 
		n[K, r, t - 1] + 
                r\ n[K, r, 
                    t - 1]\ \((1 - 
                      n[K, r, t - 1]\ /K)\)\); 
n[K_, r_, 
              0] := \ 10; 
iter[r_] := 
            Table[n[1000.0, r, t], {t, 0, timesteps}]; 
For[
            r = rbegin, r < rend, r += rstep, 
            thelist = 
              Flatten[{thelist, \(({r, #1} &)\) /@ 
                    Take[iter[r], \(-takelast\)]}, 
                1]]; 
ListPlot[thelist, 
            PlotStyle -> PointSize[0.002], 
            PlotRange -> {Automatic, {0, 2000}}]]\)

Timing[bifurcationplot[{1, 3,  0.003},  200,  20]]

Timing[bifurcationplot[{1, 3,  0.001},  200,  20]]

Timing[bifurcationplot[{1, 3,  0.003},  300,  30]]

Figure 4.2.1 - 4.2.3 - cob-webbing

rvalue = 0.7;
Kvalue = 1000;
initn = 400;

n + rvalue n (1 - n/Kvalue) /. n -> {400, 568}

\(DisplayTogether[
ListPlot[nvalue = initn; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            0}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
Plot[
        n + rvalue\ n\ \((1 - n/Kvalue)\), {n, 0, 
          1200}], 
Plot[n, {n, 0, 1200}, 
        PlotStyle -> Dashing[{0.02}]]
]\)

Make the point that a negative slope at K implies overshooting (cob-webbing generates boxes around equilibrium).

rvalue = 1.8;
Kvalue = 1000;
initn = 825;

n + rvalue n (1 - n/Kvalue) /. n -> 825

D[n + rvalue n (1 - n/Kvalue), n] /. n -> 1000

\(DisplayTogether[
ListPlot[nvalue = initn; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            0}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
Plot[
        n + rvalue\ n\ \((1 - n/Kvalue)\), {n, 0, 
          1200}], 
Plot[n, {n, 0, 1200}, 
        PlotStyle -> Dashing[{0.02}]]
]\)

Cob-webbing with an extremely negative slope around K implies that the oscillations increase (the boxes around the equilibrium are getting larger).

rvalue = 2.5;
Kvalue = 1000;
initn = 950;

n + rvalue n (1 - n/Kvalue) /. n -> 825

D[n + rvalue n (1 - n/Kvalue), n] /. n -> 1000

\(DisplayTogether[
ListPlot[nvalue = initn; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            0}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
ListPlot[
        nvalue = newnvalue; 
        newnvalue = 
          nvalue + rvalue\ nvalue\ \((1 - nvalue/Kvalue)\); {{nvalue, 
            nvalue}, {nvalue, newnvalue}, {newnvalue, newnvalue}}, 
        PlotJoined -> True], 
Plot[
        n + rvalue\ n\ \((1 - n/Kvalue)\), {n, 0, 
          1200}], 
Plot[n, {n, 0, 1200}, 
        PlotStyle -> Dashing[{0.02}]]
]\)