Generating the Figures in Chapter 3 of Otto and Day

Clear["Global`*"]

Figure 3.2 - Pheasant data from Lack (2001)

dates=Table[x,{x,1937,1942}]

pheasants={8,30,81,282,705,1325};

data=Transpose[{dates,pheasants}]

Show[ListPlot[Table[{1937+i,8 (R^i)/.R->3.5`},{i,0,5}],Joined->True,PlotStyle->Dashing[{0.04`}]],ListPlot[Table[{1937+i,8 (R^i)/.R->3},{i,0,5}],Joined->True],ListPlot[Table[{1937+i,8 (R^i)/.R->2.5`},{i,0,5}],Joined->True,PlotStyle->Dashing[{0.01`}]],ListPlot[data,PlotStyle->PointSize[0.02`]],PlotRange->{{1936.5`,1942.5`},{-20,1400}},AxesOrigin->{1936.5`,0}]

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

8*3^i*1.5/.i->50

Table[{1937+i,8*R^i/.R->3},{i,0,5}]

Figure 3.3 - Number of offspring per parent according to the logistic

Plot[(1+r )+(-r/K)n/.K->100/.r->1,{n,0,200},PlotRange->{Automatic,{0,2.1}}]

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

Figure 3.4 - Yeast data from Mable and Otto (2001)

The cell density per milliliter at various points in time for two replicate haploid strains of yeast {time, cell density}:

LogGrowthC41a:={
	{0,3.1*10^(5)},
	{2,4.8*10^(5)},
	{4,1.8*10^(6)},
	{6,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)},
	{6,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)}
	};

The cell density per milliliter at various points in time for two replicate diploid strains of yeast:

LogGrowthC42a:={
	{0,2.0*10^(5)},
	{2,4.5*10^(5)},
	{4,1.4*10^(6)},
	{6,3.5*10^(6)},
	{8,1.0*10^(7)},
	{10,2.7*10^(7)},
	{12,7.6*10^(7)},
	{14,1.6*10^(8)},
	{16,1.9*10^(8)},
	{18,1.9*10^(8)},
	{20,1.9*10^(8)},
	{22,2.0*10^(8)},
	{24,2.2*10^(8)},
	{26,1.8*10^(8)},
	{28,2.1*10^(8)},
	{30,2.2*10^(8)},
	{32,2.0*10^(8)},
	{34,2.0*10^(8)},
	{36,2.3*10^(8)},
	{48,2.2*10^(8)},
	{72,2.5*10^(8)}
	};

LogGrowthC42b:={
	{0,1.8*10^(5)},
	{2,6.8*10^(5)},
	{4,2.0*10^(6)},
	{6,3.3*10^(6)},
	{8,1.2*10^(7)},
	{10,2.8*10^(7)},
	{12,9.6*10^(7)},
	{14,2.0*10^(8)},
	{16,2.3*10^(8)},
	{18,2.1*10^(8)},
	{20,1.9*10^(8)},
	{22,1.9*10^(8)},
	{24,2.3*10^(8)},
	{26,2.1*10^(8)},
	{28,2.1*10^(8)},
	{30,2.1*10^(8)},
	{32,1.9*10^(8)},
	{34,2.1*10^(8)},
	{36,2.0*10^(8)},
	{48,2.9*10^(8)},
	{72,2.5*10^(8)}
	};

Plotting the haploid data:

axesrange={{0,73},{-10^7,4.2 10^8}};

Plot1=ListPlot[{LogGrowthC41a,LogGrowthC41b},PlotRange->axesrange,PlotMarkers->{{■,15},{"\[Cross]",15}}]

Plotting the diploid data:

Plot2=ListPlot[{LogGrowthC42a,LogGrowthC42b},PlotRange->axesrange,PlotMarkers->{{■,15},{"\[Cross]",15}}]

These plots make it clear that the carrying capacity hovers around 3.5×(10^8)in haploids and 2×(10^8)in diploids.  To be more precise, Mable and Otto (2001) set the carrying capacity to the average cell density observed after hour 32 (i.e., the 18th through 21st elements in the data arrays):

For haploids:

K41a=Sum[LogGrowthC41a[[i,2]],{i,18,21}]/4

K41b=Sum[LogGrowthC41b[[i,2]],{i,18,21}]/4

K41=(K41a+K41b)/2

For diploids:

K42a=Sum[LogGrowthC42a[[i,2]],{i,18,21}]/4

K42b=Sum[LogGrowthC42b[[i,2]],{i,18,21}]/4

K42=(K42a+K42b)/2

In Mable and Otto (2001), the growth rates were estimated by transforming the logistic equation into a line and then performing a linear regression (finding the line that minimized the sum of squared deviations of the data from the line).

In Chapter 6, we will show that:

	(Subscript[N, t]) = K/(1+(K/(Subscript[N, 0])-1)*Exp[–r*t])
	
satisfies the differential equation (3.5c) for the logistic model.  This equation gives us the value of N at time t (written here as (Subscript[N, t])).  It can be rearranged, as follows, to give a linear function of time, t:

Step 1:  (1+(K/(Subscript[N, 0])–1)*Exp[-r*t]) = K/(Subscript[N, t])  	{Divide both sides by Nt and multiply by (1+(K/(Subscript[N, 0])-1)*Exp[-r*t])}
Step 2:  (K/(Subscript[N, 0])–1)*Exp[–r*t] = (K/(Subscript[N, t])– 1) 	{Subtract 1 from both sides}
Step 3:  Exp[–r*t] = (K/(Subscript[N, t])– 1)/(K/(Subscript[N, 0])–1) 	{Divide both sides by (K/(Subscript[N, 0])-1)}
Step 4:  r*t  = –Log[(K/(Subscript[N, t]) – 1)/(K/(Subscript[N, 0])–1)] 	{Take the log of both sides and divide by –1}

The data are plugged into the right hand-side at various time points, t, and the best fitting line with slope r is found.

Only time points between 4 hours and 12 hours were used (earlier time points were subject to measurement error while later time points were close to the carrying capacity, which can cause the above transformation to generate complex numbers).

Estimating r for haploids (Transpose[LogGrowthC41a][[2]] turns the rows into columns for the data in LogGrowthC41a and then chooses the second row giving all of the cell densities):

-Log[(K41a/Transpose[LogGrowthC41a][[2]]-1)/(K41a/LogGrowthC41a[[1,2]]-1)]

Putting back in the times:

Transpose[{Transpose[LogGrowthC41a][[1]],%}]

haploida=Table[%[[i]],{i,3,7}]

Using Mathematica's regression tools:

hapamodel=LinearModelFit[haploida,{1,t},t]
hapamodel["ParameterTable"]

hapamodel["BestFit"]
hapra=%[[2]]/(t)

The second replication haploid data:

-Log[(K41b/Transpose[LogGrowthC41b][[2]]-1)/(K41b/LogGrowthC41b[[1,2]]-1)]

Putting back in the times:

Transpose[{Transpose[LogGrowthC41b][[1]],%}]

haploidb=Table[%[[i]],{i,3,7}]

hapbmodel=LinearModelFit[haploidb,{1,t},t]
hapbmodel["ParameterTable"]

hapbmodel["BestFit"]
haprb=%[[2]]/(t)

The average haploid growth rate:

r41=(hapra+haprb)/2

Estimating r for diploids:

-Log[(K42a/Transpose[LogGrowthC42a][[2]]-1)/(K42a/LogGrowthC42a[[1,2]]-1)]

Putting back in the times:

Transpose[{Transpose[LogGrowthC42a][[1]],%}]

diploida=Table[%[[i]],{i,3,7}]

dipamodel=LinearModelFit[diploida,{1,t},t]
dipamodel["ParameterTable"]

dipamodel["BestFit"]
dipra=%[[2]]/(t)

The second replication haploid data:

-Log[(K42b/Transpose[LogGrowthC42b][[2]]-1)/(K42b/LogGrowthC42b[[1,2]]-1)]

Putting back in the times:

Transpose[{Transpose[LogGrowthC42b][[1]],%}]

diploidb=Table[%[[i]],{i,3,7}]

dipbmodel=LinearModelFit[diploidb,{1,t},t]
dipbmodel["ParameterTable"]

dipbmodel["BestFit"]
diprb=%[[2]]/(t)

The average diploid growth rate:

r42=(dipra+diprb)/2

The average initial population sizes for the haploids (n10) and diploids (n20):

n10=(LogGrowthC41a[[1,2]]+LogGrowthC41b[[1,2]])/2
n20=(LogGrowthC42a[[1,2]]+LogGrowthC42b[[1,2]])/2

The solution to the logistic model is given in Chapter 6, see equation (6.14b):

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

Show[
Plot1,
Plot[logistic[K41,r41,n10,t],{t,0,100},PlotRange->axesrange]
]

Show[
Plot2,
Plot[logistic[K42,r42,n20,t],{t,0,100},PlotRange->axesrange]
]

Alternatively, we can use Mathematica's NonlinearFit function to fit the logistic in a way that minimizes the sum of squared deviations from the logistic curve.

This procedure generates similar estimates of r of 0.547855 and 0.553945 for haploids (using the first data point as the initial number of cells:

functhapa=NonlinearModelFit[LogGrowthC41a,(K41a/(1+((K41a/(LogGrowthC41a[[1,2]]))-1) (E^(-r t)))),{r},t]["BestFit"]

functhapb=NonlinearModelFit[LogGrowthC41b,(K41b/(1+((K41b/(LogGrowthC41b[[1,2]]))-1) (E^(-r t)))),{r},t]["BestFit"]

Show[
Plot1,
Plot[functhapa,{t,0,100},PlotRange->axesrange],
Plot[functhapb,{t,0,100},PlotRange->axesrange]
]

...and similar estimates of r of 0.538108 and 0.583806 for diploids:

functdipa=NonlinearModelFit[LogGrowthC42a,(K42a/(1+((K42a/(LogGrowthC42a[[1,2]]))-1) (E^(-r t)))),{r},t]["BestFit"]

functdipb=NonlinearModelFit[LogGrowthC42b,(K42b/(1+((K42b/(LogGrowthC42b[[1,2]]))-1) (E^(-r t)))),{r},t]["BestFit"]

Show[
Plot2,
Plot[functdipa,{t,0,100},PlotRange->axesrange],
Plot[functdipb,{t,0,100},PlotRange->axesrange]
]

Figure 3.6 - Allele frequency change in the haploid model of selection

Show[Plot[s p (1-p)/.s->0.1,{p,0,1},PlotStyle->Thickness[0.01]],Plot[s p (1-p)/(1+s p)/.s->0.1,{p,0,1},PlotStyle->Dashing[{0.04,0.04}]],PlotRange->{{0,1},{0,0.03}}]

Figure 3.10 - Comparing contact rates

Show[
Plot[c s/.c->0.001,{s,0,10000},PlotStyle->Thickness[0.01]],
Plot[c s/n/.c->6/.n->s/0.9,{s,0,10000}],
Plot[c s/(b+n)/.c->6/.b->500/.n->s/0.9,{s,0,10000},PlotStyle->Dashing[{0.01}],
PlotRange->{{0,10000},{0,10}}]
]

Export["test.eps",%,"EPS"]