Generating the Figures in Chapter 13 of Otto and Day

Inputs (Enter first)

<< Graphics`Graphics`

<< Graphics`MultipleListPlot`

<<Statistics`DiscreteDistributions`

<<Statistics`ContinuousDistributions`

Off[General::"spell"]

Off[General::"spell1"]

Figure 13.1 - Family size distributions

After these packages have been loaded, Mathematica knows about the PoissonDistribution:

?PoissonDistribution

BarChart[Table[{PDF[PoissonDistribution[1.2],x],x},{x,0,6}]]

Figure 13.2 - Exponential Growth

SeedRandom[1219]

Table[Random[PoissonDistribution[1.2]], {x, 1, 10}]

Sum[%[[i]], {i, 1, 10}]

Table[Random[PoissonDistribution[1.2]], {x, 1, 10}]

Sum[%[[i]], {i, 1, 10}]

Random[PoissonDistribution[1.2 *10]]

Random[PoissonDistribution[1.2*21]]

1.2*21

Let's now develop a Mathematica program that simulates logistic growth over time
when the number of children that a parent has is distributed according to the Poisson.

The Poisson distribution is special in that the sum of n Poisson distributions with mean
μ is equal to a Poisson distribution with mean n μ.  We will use this fact in the following
way.  At time t, there are n[t] parents that each have μ kids on average.
If the number of kids per parent is Poisson distributed, then the number of kids produced
by the whole population will be Poisson distributed with mean n[t]*μ.

Clear[num, t, μ, n0, trial]
num[t_,μ_,n0_,trial_]:=num[t,μ,n0,trial]=If[num[t-1,μ,n0,trial]>0,
    Random[PoissonDistribution[μ num[t-1,μ,n0,trial]]],0]
num[0,μ_,n0_,trial_]:=n0

The "If" statement is needed to account for the possibility that the population might go extinct, in which case the
number of offspring remains zero (note that PoissonDistribution is not defined for μ = 0).
The variable "trial" that has been included so that we can run more than one trial with the same set of parameters.  To do this, just 
run the program once setting trial to one and then repeat setting trial to 2, 3, etc:

mu0 = 1.2;
num0 = 10;

DisplayTogether[
ListPlot[Table[{t,num[t,mu0,num0,1]},{t,0,30}],
	PlotJoined->True,PlotStyle->Thickness[0.02]],
ListPlot[Table[{t,num[t,mu0,num0,2]},{t,0,30}],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
ListPlot[Table[{t,num[t,mu0,num0,3]},{t,0,30}],
	PlotJoined->True],
ListPlot[Table[{t,num[t,mu0,num0,4]},{t,0,30}],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,num[t,mu0,num0,5]},{t,0,30}],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
	PlotRange->{{0,30},{0,2000}}
	]

DisplayTogether[
ListPlot[Table[{t,num[t,mu0,num0,1]},{t,0,30}],
	PlotJoined->True,PlotStyle->Thickness[0.02]],
ListPlot[Table[{t,num[t,mu0,num0,2]},{t,0,30}],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
ListPlot[Table[{t,num[t,mu0,num0,3]},{t,0,30}],
	PlotJoined->True],
ListPlot[Table[{t,num[t,mu0,num0,4]},{t,0,30}],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,num[t,mu0,num0,5]},{t,0,30}],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
	PlotRange->{{0,5},{0,50}}
	]

tab = Table[num[30, mu0, num0, j], {j, 1, 100}]

Count[tab, 0]

Quantile[tab, 0.025]

Quantile[tab, 0.975]

Mean[tab] // N

Sqrt[Variance[tab]] // N

Sqrt[Variance[tab]/100] // N

10 1.2^30

Figure 13.3 - Logistic Growth

SeedRandom[0706]

Clear[num]

10 (1 + 0.2 (1 - 10/100))

Random[PoissonDistribution[%]]

% (1 + 0.2 (1 - %/100))

Random[PoissonDistribution[%]]

Let's now develop a Mathematica program that simulates logistic growth over time
when the number of children that a parent has is distributed according to the Poisson.

The Poisson distribution is special in that the sum of n Poisson distributions with mean
μ is equal to a Poisson distribution with mean n μ.  We will use this fact in the following
way.  At time t, there are n[t] parents that each have μ kids on average.
If the number of kids per parent is Poisson distributed, then the number of kids produced
by the whole population will be Poisson distributed with mean n[t]*μ.

num[t_,r_,k_,n0_,trial_]:=num[t,r,k,n0,trial]=If[num[t-1,r,k,n0,trial]>0,
    Random[PoissonDistribution[
	num[t-1,r,k,n0,trial]+r*num[t-1,r,k,n0,trial]*
	(1-num[t-1,r,k,n0,trial]/k)]],0]
num[0,r_,k_,n0_,trial_]:=n0

The "If" statement is needed to account for the possibility that the population might go extinct, in which case the
number of offspring remains zero (note that PoissonDistribution is not defined for μ = 0).
The variable "trial" that has been included so that we can run more than one trial with the same set of parameters.  To do this, just 
run the program once setting trial to one and then repeat setting trial to 2, 3, etc:

DisplayTogether[
ListPlot[Table[{t,num[t,0.2,100,10,1]},{t,0,50}],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,num[t,0.2,100,10,2]},{t,0,50}],
	PlotJoined->True,PlotStyle->Thickness[0.02]],
ListPlot[Table[{t,num[t,0.2,100,10,3]},{t,0,50}],
	PlotJoined->True],
ListPlot[Table[{t,num[t,0.2,100,10,4]},{t,0,50}],
	PlotJoined->True,PlotStyle->Dashing[{0.025}]],
ListPlot[Table[{t,num[t,0.2,100,10,5]},{t,0,50}],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
	PlotRange->{{0,50},{0,140}}
	]

DisplayTogether[
ListPlot[Table[{t,num[t,0.2,100,10,1]},{t,0,100}],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,num[t,0.2,100,10,2]},{t,0,100}],
	PlotJoined->True,PlotStyle->Thickness[0.016]],
ListPlot[Table[{t,num[t,0.2,100,10,3]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,num[t,0.2,100,10,4]},{t,0,100}],
	PlotJoined->True,PlotStyle->Dashing[{0.025}]],
ListPlot[Table[{t,num[t,0.2,100,10,5]},{t,0,100}],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
	PlotRange->{{0,100},{0,150}}
	]

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

Figure 13.4 - Exponential Growth with bernoulli environmental stochasticity

SeedRandom[7261]

Clear[num, t, Rg, Rb, p, n0, trial]
num[t_,Rg_,Rb_,p_,n0_,trial_]:=num[t,Rg,Rb,p,n0,trial]=If[num[t-1,Rg,Rb,p,n0,trial]>0,
   yeartype=If[Random[]<p,Rg, Rb];
    Random[PoissonDistribution[yeartype num[t-1,Rg,Rb,p,n0,trial]]],0]
num[0,Rg_,Rb_,p_,n0_,trial_]:=n0

The "If[Random[]<p,Rg, Rb]" statement determines whether the current year is a bad or a good one.

Rg0 = 1.5;
Rb0 = 0.5;
mu0 = 1.2;
num0 = 10;

p0 = 
      p /. (Solve[mu0  == p Rg0 + (1 - p) Rb0, p] // 
            Flatten)

The implied variance is:

p0 (Rg0 - mu0)^2 + (1 - p0) (Rb0 - mu0)^2

DisplayTogether[
ListPlot[Table[{t,num[t,Rg0,Rb0,p0,num0,1]},{t,0,30}],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,num[t,Rg0,Rb0,p0,num0,2]},{t,0,30}],
	PlotJoined->True],
ListPlot[Table[{t,num[t,Rg0,Rb0,p0,num0,3]},{t,0,30}],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
ListPlot[Table[{t,num[t,Rg0,Rb0,p0,num0,4]},{t,0,30}],
	PlotJoined->True],
ListPlot[Table[{t,num[t,Rg0,Rb0,p0,num0,5]},{t,0,30}],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
	PlotRange->{{0,30},{0,2000}}
	]

One line goes extinct at generation 8:

Table[{t, num[t, Rg0, Rb0, p0, num0, 4]}, {t, 0, 30}]

tab = Table[num[30, Rg0, Rb0, p0, num0, j], {j, 1, 100}]

Count[tab, 0]

Quantile[tab, 0.025]

Quantile[tab, 0.975]

Mean[tab] // N

Sqrt[Variance[tab]] // N

The coefficient of variation is:

Sqrt[Variance[tab]/Mean[tab]] // N

The standard error of the mean is:

Sqrt[Variance[tab]/100] // N

Rg0^p0 Rb0^(1 - p0)

10 (1.07883)^30

Figure 13.5 - Chaos vs. stochasticity?

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

DisplayTogether[
      ListPlot[Table[{t, num[t, 2.4, 100, 10, 1]}, {t, 0, 50}], 
        PlotJoined -> True], 
      PlotRange -> {{0, 50}, {0, 140}}]

ListPlot[Table[numD[100, 2.7, t], {t, 0, 50}], 
      PlotJoined -> True]

ListPlot[
      Table[{num[t - 1, 2.4, 100, 10, 
            1], (num[t, 2.4, 100, 10, 1] - 
              num[t - 1, 2.4, 100, 10, 1])}, {t, 1, 50}]]

?FourierTransform

ListPlot[
      Table[{numD[100, 2.7, 
            t - 1], (numD[100, 2.7, t] - numD[100, 2.7, t - 1])}, {t, 1, 
          50}]]

Figure 13.6 - Birth model

SeedRandom[1915]

Clear[n, t, b, n0, trial,time,sumtime]
time[n0_,b_,n_,trial_]:=time[n0,b,n,trial]=Random[ExponentialDistribution[(n0+n-1) b]]
sumtime[n0_,b_,n_,trial_]:=sumtime[n0,b,n,trial]=sumtime[n0,b,n-1,trial]+time[n0,b,n,trial]
sumtime[n0_,b_,0,trial_]:=0

b0 = 1/90;

time[1, b0, 1, 1]

time[2, b0, 1, 1]

time[3, b0, 1, 1]

b0 = 1/90;
num0 = 10;

In the following, we make a table that remains at num0 until the first birth, and then moves up to num0+1, etc.

DisplayTogether[
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,1],num0+i},{sumtime[num0,b0,i+1,1],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,2],num0+i},{sumtime[num0,b0,i+1,2],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,3],num0+i},{sumtime[num0,b0,i+1,3],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,4],num0+i},{sumtime[num0,b0,i+1,4],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,5],num0+i},{sumtime[num0,b0,i+1,5],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Thickness[0.02]],
PlotRange->{{0,50},{0,30}}
]

DisplayTogether[
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,1],num0+i},{sumtime[num0,b0,i+1,1],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,2],num0+i},{sumtime[num0,b0,i+1,2],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,3],num0+i},{sumtime[num0,b0,i+1,3],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,4],num0+i},{sumtime[num0,b0,i+1,4],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,5],num0+i},{sumtime[num0,b0,i+1,5],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Thickness[0.02]],
PlotRange->{{0,200},{0,160}}
]

DisplayTogether[
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,1],num0+i},{sumtime[num0,b0,i+1,1],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.01,0.02}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,2],num0+i},{sumtime[num0,b0,i+1,2],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,3],num0+i},{sumtime[num0,b0,i+1,3],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,4],num0+i},{sumtime[num0,b0,i+1,4],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,i,5],num0+i},{sumtime[num0,b0,i+1,5],
        num0+i}},{i,0,3000}],1],
	PlotJoined->True,PlotStyle->Thickness[0.02]],
PlotRange->{{0,400},{0,1500}}
]

Figure 13.7 - Birth-death model

SeedRandom[2101]

Clear[n,t,b,d,n0,trial,time,sumtime,pop]
pop[n0_,b_,d_,n_,trial_]:=pop[n0,b,d,n,trial]=
	If[pop[n0,b,d,n-1,trial]>0,pop[n0,b,d,n-1,trial]+If[ Random[]<b/(b+d),1,-1],0]
pop[n0_,b_,d_,0,trial_]:=n0
time[n0_,b_,d_,n_,trial_]:=time[n0,b,d,n,trial]=
	If[pop[n0,b,d,n-1,trial]>0,Random[ExponentialDistribution[pop[n0,b,d,n-1,trial] (b+d)]],1]
sumtime[n0_,b_,d_,n_,trial_]:=sumtime[n0,b,d,n,trial]=sumtime[n0,b,d,n-1,trial]+time[n0,b,d,n,trial]
sumtime[n0_,b_,d_,0,trial_]:=0

We hold the expected number of births minus deaths in a unit time constant by setting b0-d0 = 1/90:

b0 = 21/90;
d0 = 20/90;
num0 = 10;

DisplayTogether[
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,1],
  pop[num0,b0,d0,i,1]},{sumtime[num0,b0,d0,i+1,1],pop[num0,b0,d0,i,1]}},{i,0,12000,20}],1],
	PlotJoined->True],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,2],
  pop[num0,b0,d0,i,2]},{sumtime[num0,b0,d0,i+1,2],pop[num0,b0,d0,i,2]}},{i,0,12000,20}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,3],
  pop[num0,b0,d0,i,3]},{sumtime[num0,b0,d0,i+1,3],pop[num0,b0,d0,i,3]}},{i,0,12000,20}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.02}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,4],
  pop[num0,b0,d0,i,4]},{sumtime[num0,b0,d0,i+1,4],pop[num0,b0,d0,i,4]}},{i,0,12000,20}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,5],
  pop[num0,b0,d0,i,5]},{sumtime[num0,b0,d0,i+1,5],pop[num0,b0,d0,i,5]}},{i,0,12000,20}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
PlotRange->{{0,200},{0,400}}
]

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

DisplayTogether[
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,1],
  pop[num0,b0,d0,i,1]},{sumtime[num0,b0,d0,i+1,1],pop[num0,b0,d0,i,1]}},{i,0,800}],1],
	PlotJoined->True],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,2],
  pop[num0,b0,d0,i,2]},{sumtime[num0,b0,d0,i+1,2],pop[num0,b0,d0,i,2]}},{i,0,500}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.04}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,3],
  pop[num0,b0,d0,i,3]},{sumtime[num0,b0,d0,i+1,3],pop[num0,b0,d0,i,3]}},{i,0,600}],1],
	PlotJoined->True,PlotStyle->Dashing[{0.02}]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,4],
  pop[num0,b0,d0,i,4]},{sumtime[num0,b0,d0,i+1,4],pop[num0,b0,d0,i,4]}},{i,0,200}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Flatten[Table[{{sumtime[num0,b0,d0,i,5],
  pop[num0,b0,d0,i,5]},{sumtime[num0,b0,d0,i+1,5],pop[num0,b0,d0,i,5]}},{i,0,600}],1],
	PlotJoined->True,PlotStyle->Thickness[0.01]],
PlotRange->{{0,50},{0,69}}
]

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

ListPlot[
      Flatten[Table[{{sumtime[num0, b0, d0, i, 4], 
              pop[num0, b0, d0, i, 4]}, {sumtime[num0, b0, d0, i + 1, 4], 
              pop[num0, b0, d0, i, 4]}}, {i, 0, 60, 1}], 1], 
      PlotJoined -> True, PlotStyle -> Thickness[0.01]]

Sum[(d/b)^i, {i, 0, Infinity}]

Simplify[
      1 - Sum[(d/b)^i, {i, 0, x - 1}]/
          Sum[(d/b)^i, {i, 0, Infinity}]]

% /. x -> {1, 10} /. d -> d0 /. b -> b0 // 
      N

tab = 
      Table[For[i = 0, i < 10001, i++, pop[num0, b0, d0, i, j]]; 
        pop[num0, b0, d0, 10000, j], {j, 1, 100}]

Count[tab, 0]

Quantile[tab, 0.025]

Quantile[tab, 0.975]

The standard error of the proportion of extinctions:

Sqrt[(Count[tab, 0]/100) ((100 - Count[tab, 0])/100)/100] // 
      N

Figure 13.8 - Wright-Fisher model

SeedRandom[1146]

?BinomialDistribution

Clear[p,num,t,trial,allele]
allele[p_,num_,t_,trial_]:=allele[p,num,t,trial]=
	Random[BinomialDistribution[num,allele[p,num,t-1,trial]]]/num
allele[p_,num_,0,trial_]:=p

num0 = 100;
p0 = 1/2;

allele[p0, num0, 1, 1]

DisplayTogether[
ListPlot[Table[{t,allele[p0,num0,t,1]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,t,2]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,t,3]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,t,4]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,t,5]},{t,0,100}],
	PlotJoined->True],
	PlotRange->{{0,100},{0,1}}
	]

Table[{t, allele[p0, num0, t, 1]}, {t, 0, 100}]

Table[{t, allele[p0, num0, t, 2]}, {t, 0, 100}]

Table[{t, allele[p0, num0, t, 5]}, {t, 0, 100}]

For[j = 1; lost = 0; fixed = 0; sumgen = 0; sumgen2 = 0, 
      j < 501,  j++,  
gen = 
        0; 
segregating = 1; 
While[
        segregating == 1, 
        gen = gen + 1; 
If[
          allele[p0, num0, gen, j] == 0, segregating = 0; lost++, 
          If[allele[p0, num0, gen, j] == 1, 
            segregating = 0; fixed++]]]; 
sumgen = 
        sumgen + gen; 
sumgen2 = 
        sumgen2 + gen^2
]

lost

fixed

meangen = sumgen/500 // N

vargen = (sumgen2 - sumgen^2/500)/(500 - 1) // N

The standard error is:

Sqrt[vargen/500]

242/500 // N

The probability of observing 242 losses or something more extreme is:

Sum[
        PDF[BinomialDistribution[500, 0.5], i], {i, 0, 
          242}] + 
Sum[
        PDF[BinomialDistribution[500, 0.5], i], {i, 258, 500}]

Figure 13.9 - Wright-Fisher model with selection

SeedRandom[1615]

?BinomialDistribution

Clear[p,num,WA,Wa,t,trial,allele]
allele[p_,num_,WA_,Wa_,t_,trial_]:=allele[p,num,WA,Wa,t,trial]=
	Random[BinomialDistribution[num,WA allele[p,num,WA,Wa,t-1,trial]/(WA allele[p,num,WA,Wa,t-1,trial]+Wa (1-allele[p,num,WA,Wa,t-1,trial]))]]/num
allele[p_,num_,WA_,Wa_,0,trial_]:=p

We hold the expected number of births minus deaths in a unit time constant by setting b0-d0 = 1/90:

num0 = 100;
p0 = 0.05;
WA0 = 1.1;
Wa0 = 1;

allele[p0, num0, WA0, Wa0, 1, 1]

The deterministic expectation is:

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.05

DisplayTogether[
ListPlot[Table[{t,hap[0.1,t]},{t,0,100}],PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,1]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,2]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,3]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,4]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,5]},{t,0,100}],
	PlotJoined->True],
	PlotRange->{{0,100},{0,1}}
	]

DisplayTogether[
      ListPlot[Table[{t, hap[0.1, t]}, {t, 0, 100}], PlotJoined -> True, 
        PlotStyle -> Thickness[0.01]], 
      ListPlot[Table[{t, allele[p0, num0, WA0, Wa0, t, 1]}, {t, 0, 100}], 
        PlotJoined -> True], 
      ListPlot[Table[{t, allele[p0, num0, WA0, Wa0, t, 2]}, {t, 0, 100}], 
        PlotJoined -> True], 
      ListPlot[Table[{t, allele[p0, num0, WA0, Wa0, t, 3]}, {t, 0, 100}], 
        PlotJoined -> True], 
      ListPlot[Table[{t, allele[p0, num0, WA0, Wa0, t, 4]}, {t, 0, 100}], 
        PlotJoined -> True], 
      ListPlot[Table[{t, allele[p0, num0, WA0, Wa0, t, 5]}, {t, 0, 100}], 
        PlotJoined -> True], 
      PlotRange -> {{0, 10}, {0, 0.21}}]

Table[{t, allele[p0, num0, WA0, Wa0, t, 2]}, {t, 0, 10}]

num0 = 10000;
p0 = 0.05;
WA0 = 1.1;
Wa0 = 1;

DisplayTogether[
ListPlot[Table[{t,hap[0.1,t]},{t,0,100}],PlotJoined->True,PlotStyle->Thickness[0.01]],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,1]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,2]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,3]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,4]},{t,0,100}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,5]},{t,0,100}],
	PlotJoined->True],
	PlotRange->{{0,100},{0,1}}
	]

Figure 13.10 - Moran model

SeedRandom[1146]

Clear[ran,pfit,p,num,WA,Wa,t,trial,allele]
allele[p_,num_,WA_,Wa_,t_,trial_]:=allele[p,num,WA,Wa,t,trial]=
	If[
	pfit=(WA allele[p,num,WA,Wa,t-1,trial])/(WA allele[p,num,WA,Wa,t-1,trial]+Wa 		
		(1-allele[p,num,WA,Wa,t-1,trial]));
	ran=Random[];
	ran<pfit*(1-allele[p,num,WA,Wa,t-1,trial]),
		allele[p,num,WA,Wa,t-1,trial]+1/num,
			If[ran>1-(1-pfit)*allele[p,num,WA,Wa,t-1,trial],allele[p,num,WA,Wa,t-1,trial]-1/num,
				allele[p,num,WA,Wa,t-1,trial]]]
allele[p_,num_,WA_,Wa_,0,trial_]:=p

We hold the expected number of births minus deaths in a unit time constant by setting b0-d0 = 1/90:

num0 = 100;
p0 = 1/2;
WA0 = 1;
Wa0 = 1;

DisplayTogether[
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,1]},{t,0,10000,20}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,2]},{t,0,10000,20}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,3]},{t,0,10000,20}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,4]},{t,0,10000,20}],
	PlotJoined->True],
ListPlot[Table[{t,allele[p0,num0,WA0,Wa0,t,5]},{t,0,10000,20}],
	PlotJoined->True],
	PlotRange->{{0,10000},{0,1}}
	]

DisplayTogether[
      ListPlot[Table[{t, allele[p0, num0, WA0, Wa0, t, 4]}, {t, 9900, 
            10000}], PlotJoined -> True], 
      PlotRange -> {{9975, 10000}, {0.50, 0.56}}]

For[tt = 1, 
      tt < 6, tt++, 
test = 
        0; 
While[
          allele[p0, num0, WA0, Wa0, test, tt] > 0 && 
            allele[p0, num0, WA0, Wa0, test, tt] < 
              1, test++]
Print[{test, 
            allele[p0, num0, WA0, Wa0, test, tt]}]]

For[j = 1; lost = 0; fixed = 0; sumgen = 0; sumgen2 = 0, 
      j < 501,  j++,  
gen = 
        0; 
segregating = 1; 
While[
        segregating == 1, 
        gen = gen + 1; 
If[
          allele[p0, num0, WA0, Wa0, gen, j] == 0, 
          segregating = 0; lost++, 
          If[allele[p0, num0, WA0, Wa0, gen, j] == 1, 
            segregating = 0; fixed++]]]; 
sumgen = 
        sumgen + gen; 
sumgen2 = 
        sumgen2 + gen^2
]

lost

fixed

sumgen

sumgen2

meangen = sumgen/500 // N

vargen = (sumgen2 - sumgen^2/500)/(500 - 1) // N

The mean in terms of generations is:

meangen/100

The variance in terms of generations is:

vargen/100^2

The standard error is:

Sqrt[%/500]

The probability of observing 242 losses or something more extreme is:

Sum[
        PDF[BinomialDistribution[500, 0.5], i], {i, 0, 
          242}] + 
Sum[
        PDF[BinomialDistribution[500, 0.5], i], {i, 258, 500}]

Figure 13.13 - Retinoblastoma

SeedRandom[1626]

tab = {0.497, 0.139, 0.118, 0.116, 0.006, 0.07, 0.046, 0.004, 
          0.004};

Sum[tab[[j]], {j, 1, 9}]

mtab = 
      m /. Table[Solve[2^m == tab[[j]] 20 10^6, m], {j, 1, 9}] // 
        Flatten

itab = IntegerPart[mtab]

ptab = 
      p /. Table[
            Solve[(1 - p) 2^itab[[j]] + p 2^(itab[[j]] + 1) == 
                tab[[j]] 20 10^6, p], {j, 1, 9}] // Flatten

Sum[
      41 - ((1 - ptab[[j]]) itab[[j]] + 
            ptab[[j]] (itab[[j]] + 1)), {j, 1, 9}]

Sum[41 - mtab[[j]], {j, 1, 9}]

totalS = 
      Sum[Sum[2^x, {x, 1, itab[[j]]}] + 
          ptab[[j]] 2^(itab[[j]] + 1), {j, 1, 9}]

muest = 1 - Exp[Log[1 - 0.95]/(2  totalS)]

nodisease = Simplify[(1 - mu)^num]

The probability that a heterozygous individual is unaffected:

nodisease /. mu -> muest /. num -> 2 totalS

Log[20 10^6]/Log[2. ]

The probability that an eye is unaffected (containing only half of the final number of cells):

nodisease /. mu -> muest /. num -> totalS

The probability of a unilateral case (see Knudson 1971):

2*(nodisease)*(1 - nodisease) /. mu -> muest /. 
      num -> totalS

The probability of a bilateral case (see Knudson 1971):

(1 - nodisease)^2 /. mu -> muest /. 
      num -> totalS

%/(% + %%)

((36 + 23*12)/2 + 23) // N

muest %

Because the mutation rate is so low, we run simulations by first asking whether any of the daughter cells in phases 1 and 2 carry a mutation (the answer is almost certainly not, in which case we need not simulate divisions during these phases exactly).

phases12[mu_, trial_] = 
        If[Random[] < (1 - mu)^(22 + 193), 0, 1];

In 100 trials, no mutations occurred during phases 1 or 2:

Max[Table[phases12[muest, ii], {ii, 1, 100}]]

The following routine simulates the number of unaffected (wildtype) cells, affected cells, and clusters after m binary cell divisions (ran reflects the number of newly mutated daughter cells) per eye:
[Note that m is the number of cell divisions that all cells undergo, while a proportion p undergo an extra division. i iterates through these m+1 cell divisions.]

Clear[mu,m,p,trial,results,ran]
results[mu_,m_,p_,i_,trial_]:=results[mu,m,p,i,trial]=If[i<=m,
	{2 results[mu,m,p,i-1,trial][[1]]-
		(ran=Random[BinomialDistribution[2 results[mu,m,p,i-1,trial][[1]],mu]]),
	2 results[mu,m,p,i-1,trial][[2]]+ran,
	results[mu,m,p,i-1,trial][[3]]+ran},
	{results[mu,m,p,i-1,trial][[1]]+Round[p*results[mu,m,p,i-1,trial][[1]]]-
		(ran=Random[BinomialDistribution[2*Round[p* results[mu,m,p,i-1,trial][[1]]],mu]]),
	results[mu,m,p,i-1,trial][[2]]+Round[p*results[mu,m,p,i-1,trial][[2]]]+ran,
	results[mu,m,p,i-1,trial][[3]]+ran}]
results[mu_,m_,p_,0,trial_]={1,0,0};

results[muest, itab[[1]], ptab[[1]], itab[[1]] + 1, 
      1]

tab = 
      Table[Sum[
          results[muest, itab[[j]], ptab[[j]], itab[[j]] + 1, 
            ii], {j, 1, 9}], {ii, 1, 100}]

plot1 = 
      Histogram[Transpose[tab][[2]], 
        HistogramCategories -> Table[i, {i, 0, 50, 1}], 
        HistogramRange -> {0, 50}]

0.18*0.18

3.74 10^(-8) 179

Sort[Transpose[tab][[2]]]

Table[Count[%, i], {i, 0, 10}]

Mean[Transpose[tab][[2]]]

ListPlot[Table[{i, 100 PDF[PoissonDistribution[%], i]}, {i, 0, 50}], 
      PlotJoined -> True, PlotRange -> All]

Show[plot1, %]

plot2 = 
      Histogram[Transpose[tab][[3]], 
        HistogramCategories -> Table[i, {i, 0, 10, 1}], 
        HistogramRange -> {0, 10}]

Sort[Transpose[tab][[3]]]

Mean[Transpose[tab][[3]]]

ListPlot[
      Table[{i + 1/2, 100 PDF[PoissonDistribution[%], i]}, {i, 0, 10}], 
      PlotJoined -> True]

Show[plot2, %]

Figure 13.15 - Extinction-recolonization

SeedRandom[1216]

Clear[sources, savepop, pop, epop, extinct, i, j, k]

initp = 0.1;
ext = 0.16;
mig = 0.05;
gridsize = 10;
finalgen = 50000;
timearray = {0, 1, 2, 3, 4, 10, 20, 30, 40, 50, 100, 200, 300, 400, 
          500};

For[i = 1, 
      i < gridsize + 1, i++, 
For[j = 1, 
        j < gridsize + 1, j++, 
pop[{i, j}] = 
          If[Random[] < initp, 1, 0]; 
savepop[{i, j}, 0] = 
          pop[{i, j}]
]]

(* 
      sources describes the locations from which migrants can come *) 
 (* 
      Internal cells on the grid have 8 sources of migrants *) 
For[i = 2, 
      i < gridsize, i++, 
For[j = 2, 
        j < gridsize, j++, 
sources[{i, j}] = {{i - 1, 
              j - 1}, {i, j - 1}, {i + 1, j - 1}, {i - 1, j}, {i, 
              j - 1}, {i - 1, j + 1}, {i, j + 1}, {i + 1, 
              j + 1}}
]]

 (* 
      Edge cells on the grid have 5 sources of migrants *) 

    For[j = 2, 
      j < gridsize, j++, 
sources[{1, j}] = {{1, 
            j - 1}, {1, j + 1}, {2, j - 1}, {2, j}, {2, 
            j + 1}}; 
sources[{j, 1}] = {{j - 1, 
            1}, {j + 1, 1}, {j - 1, 2}, {j, 2}, {j + 1, 
            2}}; 
sources[{gridsize, j}] = {{gridsize, 
            j - 1}, {gridsize, j + 1}, {gridsize - 1, j - 1}, {gridsize - 1, 
            j}, {gridsize - 1, j + 1}}; 
sources[{j, 
            gridsize}] = {{j - 1, gridsize}, {j + 1, gridsize}, {j - 1, 
            gridsize - 1}, {j, gridsize - 1}, {j + 1, 
            gridsize - 
              1}};
]

 (* Corner cells on the grid have 3 sources of migrants *) 

    sources[{1, 1}] = {{1, 2}, {2, 1}, {2, 2}};

    sources[{1, gridsize}] = {{1, gridsize - 1}, {2, gridsize}, {2, 
            gridsize - 1}};

    sources[{gridsize, 1}] = {{gridsize - 1, 1}, {gridsize, 
            2}, {gridsize - 1, 2}};

    sources[{gridsize, gridsize}] = {{gridsize - 1, gridsize}, {gridsize, 
            gridsize - 1}, {gridsize - 1, gridsize - 1}};

For[gen = 1, 
      gen < finalgen + 
          1, gen++, 

 (*Migration happens first*) 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
mpop[{i, j}] = 
            If[pop[{i, j}] == 1, 
              1, 
pool = 
                Sum[pop[sources[{i, j}][[k]]], {k, 1, 
                    Length[sources[{i, j}]]}]; 
If[
                Random[] < pool*mig, 1, 
                0]]]]; 

 (*Extinction happens next*) 
extinct = 1; 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
pop[{i, j}] = 
            If[mpop[{i, j}] == 1, 
              If[Random[] < ext, 0, extinct = 0; 1], 0]; 
          If[MemberQ[timearray, gen], 
            savepop[{i, j}, gen] = 
              pop[{i, j}]]]]; 

 (*extinct tests if the entire population is extinct*) 
If[extinct == 1, Print[gen]; 
        gen = finalgen + 1]
]

For[t = 0, t < 6, t++, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 10, t < 61, t = t + 10, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 100, t < 501, t = t + 100, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

Show[
      GraphicsArray[{{p[0], p[1], p[2], p[3], 
            p[4]}, 
{p[10], p[20], p[30], p[40], 
            p[50]}, 
{p[100], p[200], p[300], p[400], 
            p[500]}}]]

Sum[savepop[{i, j}, 500], {i, 10}, {j, 10}]

SeedRandom[1217]

Clear[sources, savepop, pop, epop, extinct, i, j, k]

initp = 1;
ext = 0.16;
mig = 0.05;
gridsize = 10;
finalgen = 50000;
timearray = {0, 1, 2, 3, 4, 10, 20, 30, 40, 50, 100, 200, 300, 400, 
          500};

For[i = 1, 
      i < gridsize + 1, i++, 
For[j = 1, 
        j < gridsize + 1, j++, 
pop[{i, j}] = 
          If[Random[] < initp, 1, 0]; 
savepop[{i, j}, 0] = 
          pop[{i, j}]
]]

(* 
      sources describes the locations from which migrants can come *) 
 (* 
      Internal cells on the grid have 8 sources of migrants *) 
For[i = 2, 
      i < gridsize, i++, 
For[j = 2, 
        j < gridsize, j++, 
sources[{i, j}] = {{i - 1, 
              j - 1}, {i, j - 1}, {i + 1, j - 1}, {i - 1, j}, {i, 
              j - 1}, {i - 1, j + 1}, {i, j + 1}, {i + 1, 
              j + 1}}
]]

 (* 
      Edge cells on the grid have 5 sources of migrants *) 

    For[j = 2, 
      j < gridsize, j++, 
sources[{1, j}] = {{1, 
            j - 1}, {1, j + 1}, {2, j - 1}, {2, j}, {2, 
            j + 1}}; 
sources[{j, 1}] = {{j - 1, 
            1}, {j + 1, 1}, {j - 1, 2}, {j, 2}, {j + 1, 
            2}}; 
sources[{gridsize, j}] = {{gridsize, 
            j - 1}, {gridsize, j + 1}, {gridsize - 1, j - 1}, {gridsize - 1, 
            j}, {gridsize - 1, j + 1}}; 
sources[{j, 
            gridsize}] = {{j - 1, gridsize}, {j + 1, gridsize}, {j - 1, 
            gridsize - 1}, {j, gridsize - 1}, {j + 1, 
            gridsize - 
              1}};
]

 (* Corner cells on the grid have 3 sources of migrants *) 

    sources[{1, 1}] = {{1, 2}, {2, 1}, {2, 2}};

    sources[{1, gridsize}] = {{1, gridsize - 1}, {2, gridsize}, {2, 
            gridsize - 1}};

    sources[{gridsize, 1}] = {{gridsize - 1, 1}, {gridsize, 
            2}, {gridsize - 1, 2}};

    sources[{gridsize, gridsize}] = {{gridsize - 1, gridsize}, {gridsize, 
            gridsize - 1}, {gridsize - 1, gridsize - 1}};

For[gen = 1, 
      gen < finalgen + 
          1, gen++, 

 (*Migration happens first*) 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
mpop[{i, j}] = 
            If[pop[{i, j}] == 1, 
              1, 
pool = 
                Sum[pop[sources[{i, j}][[k]]], {k, 1, 
                    Length[sources[{i, j}]]}]; 
If[
                Random[] < pool*mig, 1, 
                0]]]]; 

 (*Extinction happens next*) 
extinct = 1; 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
pop[{i, j}] = 
            If[mpop[{i, j}] == 1, 
              If[Random[] < ext, 0, extinct = 0; 1], 0]; 
          If[MemberQ[timearray, gen], 
            savepop[{i, j}, gen] = 
              pop[{i, j}]]]]; 

 (*extinct tests if the entire population is extinct*) 
If[extinct == 1, Print[gen]; 
        gen = finalgen + 1]
]

For[t = 0, t < 6, t++, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 10, t < 61, t = t + 10, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 100, t < 501, t = t + 100, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

Show[
      GraphicsArray[{{p[0], p[1], p[2], p[3], 
            p[4]}, 
{p[10], p[20], p[30], p[40], 
            p[50]}, 
{p[100], p[200], p[300], p[400], 
            p[500]}}]]

Sum[savepop[{i, j}, 500], {i, 10}, {j, 10}]

Slightly higher extinction rate (same seed):

SeedRandom[1218]

Clear[sources, savepop, pop, epop, extinct, i, j, k]

initp = 1;
ext = 0.18;
mig = 0.05;
gridsize = 10;
finalgen = 500;
timearray = {0, 1, 2, 3, 4, 10, 20, 30, 40, 50, 100, 200, 300, 400, 
          500};

For[i = 1, 
      i < gridsize + 1, i++, 
For[j = 1, 
        j < gridsize + 1, j++, 
pop[{i, j}] = 
          If[Random[] < initp, 1, 0]; 
savepop[{i, j}, 0] = 
          pop[{i, j}]
]]

(* 
      sources describes the locations from which migrants can come *) 
 (* 
      Internal cells on the grid have 8 sources of migrants *) 
For[i = 2, 
      i < gridsize, i++, 
For[j = 2, 
        j < gridsize, j++, 
sources[{i, j}] = {{i - 1, 
              j - 1}, {i, j - 1}, {i + 1, j - 1}, {i - 1, j}, {i, 
              j - 1}, {i - 1, j + 1}, {i, j + 1}, {i + 1, 
              j + 1}}
]]

 (* 
      Edge cells on the grid have 5 sources of migrants *) 

    For[j = 2, 
      j < gridsize, j++, 
sources[{1, j}] = {{1, 
            j - 1}, {1, j + 1}, {2, j - 1}, {2, j}, {2, 
            j + 1}}; 
sources[{j, 1}] = {{j - 1, 
            1}, {j + 1, 1}, {j - 1, 2}, {j, 2}, {j + 1, 
            2}}; 
sources[{gridsize, j}] = {{gridsize, 
            j - 1}, {gridsize, j + 1}, {gridsize - 1, j - 1}, {gridsize - 1, 
            j}, {gridsize - 1, j + 1}}; 
sources[{j, 
            gridsize}] = {{j - 1, gridsize}, {j + 1, gridsize}, {j - 1, 
            gridsize - 1}, {j, gridsize - 1}, {j + 1, 
            gridsize - 
              1}};
]

 (* Corner cells on the grid have 3 sources of migrants *) 

    sources[{1, 1}] = {{1, 2}, {2, 1}, {2, 2}};

    sources[{1, gridsize}] = {{1, gridsize - 1}, {2, gridsize}, {2, 
            gridsize - 1}};

    sources[{gridsize, 1}] = {{gridsize - 1, 1}, {gridsize, 
            2}, {gridsize - 1, 2}};

    sources[{gridsize, gridsize}] = {{gridsize - 1, gridsize}, {gridsize, 
            gridsize - 1}, {gridsize - 1, gridsize - 1}};

For[gen = 1, 
      gen < finalgen + 
          1, gen++, 

 (*Migration happens first*) 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
mpop[{i, j}] = 
            If[pop[{i, j}] == 1, 
              1, 
pool = 
                Sum[pop[sources[{i, j}][[k]]], {k, 1, 
                    Length[sources[{i, j}]]}]; 
If[
                Random[] < pool*mig, 1, 
                0]]]]; 

 (*Extinction happens next*) 
extinct = 1; 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
pop[{i, j}] = 
            If[mpop[{i, j}] == 1, 
              If[Random[] < ext, 0, extinct = 0; 1], 0]; 
          If[MemberQ[timearray, gen], 
            savepop[{i, j}, gen] = 
              pop[{i, j}]]]]; 

 (*extinct tests if the entire population is extinct*) 
If[extinct == 1, Print[gen]; 
        gen = finalgen + 1]
]

For[t = 0, t < 6, t++, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 10, t < 61, t = t + 10, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 100, t < 501, t = t + 100, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

Show[
      GraphicsArray[{{p[0], p[1], p[2], p[3], 
            p[4]}, 
{p[10], p[20], p[30], p[40], 
            p[50]}, 
{p[100], p[200], p[300], p[400], 
            p[500]}}]]

Slightly lower migration rate:

SeedRandom[1219]

Clear[sources, savepop, pop, epop, extinct, i, j, k]

initp = 1;
ext = 0.16;
mig = 0.03;
gridsize = 10;
finalgen = 500;
timearray = {0, 1, 2, 3, 4, 10, 20, 30, 40, 50, 100, 200, 300, 400, 
          500};

For[i = 1, 
      i < gridsize + 1, i++, 
For[j = 1, 
        j < gridsize + 1, j++, 
pop[{i, j}] = 
          If[Random[] < initp, 1, 0]; 
savepop[{i, j}, 0] = 
          pop[{i, j}]
]]

(* 
      sources describes the locations from which migrants can come *) 
 (* 
      Internal cells on the grid have 8 sources of migrants *) 
For[i = 2, 
      i < gridsize, i++, 
For[j = 2, 
        j < gridsize, j++, 
sources[{i, j}] = {{i - 1, 
              j - 1}, {i, j - 1}, {i + 1, j - 1}, {i - 1, j}, {i, 
              j - 1}, {i - 1, j + 1}, {i, j + 1}, {i + 1, 
              j + 1}}
]]

 (* 
      Edge cells on the grid have 5 sources of migrants *) 

    For[j = 2, 
      j < gridsize, j++, 
sources[{1, j}] = {{1, 
            j - 1}, {1, j + 1}, {2, j - 1}, {2, j}, {2, 
            j + 1}}; 
sources[{j, 1}] = {{j - 1, 
            1}, {j + 1, 1}, {j - 1, 2}, {j, 2}, {j + 1, 
            2}}; 
sources[{gridsize, j}] = {{gridsize, 
            j - 1}, {gridsize, j + 1}, {gridsize - 1, j - 1}, {gridsize - 1, 
            j}, {gridsize - 1, j + 1}}; 
sources[{j, 
            gridsize}] = {{j - 1, gridsize}, {j + 1, gridsize}, {j - 1, 
            gridsize - 1}, {j, gridsize - 1}, {j + 1, 
            gridsize - 
              1}};
]

 (* Corner cells on the grid have 3 sources of migrants *) 

    sources[{1, 1}] = {{1, 2}, {2, 1}, {2, 2}};

    sources[{1, gridsize}] = {{1, gridsize - 1}, {2, gridsize}, {2, 
            gridsize - 1}};

    sources[{gridsize, 1}] = {{gridsize - 1, 1}, {gridsize, 
            2}, {gridsize - 1, 2}};

    sources[{gridsize, gridsize}] = {{gridsize - 1, gridsize}, {gridsize, 
            gridsize - 1}, {gridsize - 1, gridsize - 1}};

For[gen = 1, 
      gen < finalgen + 
          1, gen++, 

 (*Migration happens first*) 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
mpop[{i, j}] = 
            If[pop[{i, j}] == 1, 
              1, 
pool = 
                Sum[pop[sources[{i, j}][[k]]], {k, 1, 
                    Length[sources[{i, j}]]}]; 
If[
                Random[] < pool*mig, 1, 
                0]]]]; 

 (*Extinction happens next*) 
extinct = 1; 
For[i = 1, 
        i < 11, i++, 
For[j = 1, 
          j < 11, j++, 
pop[{i, j}] = 
            If[mpop[{i, j}] == 1, 
              If[Random[] < ext, 0, extinct = 0; 1], 0]; 
          If[MemberQ[timearray, gen], 
            savepop[{i, j}, gen] = 
              pop[{i, j}]]]]; 

 (*extinct tests if the entire population is extinct*) 
If[extinct == 1, Print[gen]; 
        gen = finalgen + 1]
]

For[t = 0, t < 6, t++, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 10, t < 61, t = t + 10, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

For[t = 100, t < 501, t = t + 100, 
      p[t] = Graphics[
          Raster[Table[1 - savepop[{i, j}, t], {i, 10}, {j, 10}]], 
          Frame -> True, FrameTicks -> False]]

Show[
      GraphicsArray[{{p[0], p[1], p[2], p[3], 
            p[4]}, 
{p[10], p[20], p[30], p[40], 
            p[50]}, 
{p[100], p[200], p[300], p[400], 
            p[500]}}]]

Figure 13.18 - The Coalescent

SeedRandom[4392]

First, we identify the branching times while there are n lineages:

Note that Mathematica defines the GeometricDistribution in a different manner, giving the probability that there are X trials before the successful one (not including the successful one).  Thus the range is from 0 to infinity rather than from 1 to infinity.  We thus add one to the random draw from the geometric to get the desired probability distribution (the generation at which the coalescence occurs).

?GeometricDistribution

Clear[n, t, trial,time,sumtime]
time[num_,n_,trial_]:=time[num,n,trial]=1+Random[GeometricDistribution[Binomial[n,2]/num]]
sumtime[num_,n_,trial_]:=sumtime[num,n,trial]=sumtime[num,n-1,trial]+time[num,n,trial]
sumtime[num_,1,trial_]:=0
sumtime[num_,0,trial_]:=0

n0 = 10;
num0 = 10000;

Table[sumtime[num0, n0, try], {try, 1, 4}]

Second, we identify the branching order from the base of the tree to the tip (a random number determines which lineage branches):

For[try = 1, 
      try ≤ 
        4, try++, 
branching[
          try] = {}; 
For[bb = 1, bb < 10, bb++, 
        branching[try] = 
          Append[branching[try], 
            Random[Integer, {1, bb}]]]
]

Third, we give heights for the branches, centering new branches around the parent branch from which they split off ("space" can be chosen between ~1 and ~2 for more-pleasing spacing between the branches):

space = 1.8;

Clear[height]
For[try = 1, 
      try ≤ 
        4, try++, 
height[try, 2] = {300, 
          450}; 
For[jj = 3, 
        jj < 11, jj++, 
temp = 
          height[try, 
            jj - 1]; 
temp[[branching[
                  try][[jj - 
                    1]]]] = {temp[[branching[
                      try][[jj - 1]]]] - 50*(3/jj)^space, 
            temp[[branching[try][[jj - 1]]]] + 
              50*(3/jj)^space}; 
height[try, jj] = 
          temp // Flatten
]
]

Fourth, we make lists of lines to draw the tree (tabH gives the horizontal lines, tabV gives the vertical lines):

Clear[tabH, tabV]
For[try = 1, 
      try ≤ 4, try++, 
tabH[try] = 
        Flatten[Table[
            Table[
{{sumtime[num0, ll - 1, try] - 
                    sumtime[num0, 10, try], height[try, 
                      ll][[kk]]}, {sumtime[num0, ll, try] - 
                    sumtime[num0, 10, try], height[try, 
                      ll][[kk]]}}, {kk, 1, ll}], {ll, 2, 10}], 
          1]; 
tabV[try] = 
        Table[
{{sumtime[num0, ll - 1, try] - 
                sumtime[num0, 10, try], height[try, 
                  ll][[branching[try][[ll - 1]]]]}, {sumtime[
                  num0, ll - 1, try] - 
                sumtime[num0, 10, try], height[try, 
                  ll][[branching[try][[ll - 1]] + 1]]}}, {ll, 
            2, 10}]
]

Finally, we draw each line and use MultipleListPlot to plot each of these lines:

Clear[line]
For[mm = 1; tabjoined = Join[tabH[1], tabV[1]], 
        mm ≤ 
          Length[Join[tabH[1], tabV[1]]], mm++, 
line[
            mm] = tabjoined[[mm]]];
minheight = Min[height[1, 10]] - 20; 
    maxheight = Max[height[1, 10]] + 10;

MultipleListPlot[line[1], line[2], line[3], line[4], line[5], line[6], 
      line[7], line[8], line[9], line[10], line[11], line[12], line[13], 
      line[14], line[15], line[16], line[17], line[18], line[19], line[20], 
      line[21], line[22], line[23], line[24], line[25], line[26], line[27], 
      line[28], line[29], line[30], line[31], line[32], line[33], line[34], 
      line[35], line[36], line[37], line[38], line[39], line[40], line[41], 
      line[42], line[43], line[44], line[45], line[46], line[47], line[48], 
      line[49], line[50], line[51], line[52], line[53], line[54], line[55], 
      line[56], line[57], line[58], line[59], line[60], line[61], line[62], 
      line[63], PlotRange -> {{0, -17000}, {minheight, maxheight}}, 
      PlotJoined -> True, AxesOrigin -> {0, minheight}, 
      Axes -> {True, True}, SymbolShape -> None, 
      PlotStyle -> GrayLevel[0], 
      Ticks -> {Automatic, None}]

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

Clear[line]
For[mm = 1; tabjoined = Join[tabH[2], tabV[2]], 
        mm ≤ 
          Length[Join[tabH[1], tabV[1]]], mm++, 
line[
            mm] = tabjoined[[mm]]];
minheight = Min[height[2, 10]] - 20; 
    maxheight = Max[height[2, 10]] + 10;

MultipleListPlot[line[1], line[2], line[3], line[4], line[5], line[6], 
      line[7], line[8], line[9], line[10], line[11], line[12], line[13], 
      line[14], line[15], line[16], line[17], line[18], line[19], line[20], 
      line[21], line[22], line[23], line[24], line[25], line[26], line[27], 
      line[28], line[29], line[30], line[31], line[32], line[33], line[34], 
      line[35], line[36], line[37], line[38], line[39], line[40], line[41], 
      line[42], line[43], line[44], line[45], line[46], line[47], line[48], 
      line[49], line[50], line[51], line[52], line[53], line[54], line[55], 
      line[56], line[57], line[58], line[59], line[60], line[61], line[62], 
      line[63], PlotRange -> {{0, -17000}, {minheight, maxheight}}, 
      PlotJoined -> True, AxesOrigin -> {0, minheight}, 
      Axes -> {True, True}, SymbolShape -> None, 
      PlotStyle -> GrayLevel[0], 
      Ticks -> {Automatic, None}]

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

Clear[line]
For[mm = 1; tabjoined = Join[tabH[3], tabV[3]], 
        mm ≤ 
          Length[Join[tabH[1], tabV[1]]], mm++, 
line[
            mm] = tabjoined[[mm]]];
minheight = Min[height[3, 10]] - 20; 
    maxheight = Max[height[3, 10]] + 10;

MultipleListPlot[line[1], line[2], line[3], line[4], line[5], line[6], 
      line[7], line[8], line[9], line[10], line[11], line[12], line[13], 
      line[14], line[15], line[16], line[17], line[18], line[19], line[20], 
      line[21], line[22], line[23], line[24], line[25], line[26], line[27], 
      line[28], line[29], line[30], line[31], line[32], line[33], line[34], 
      line[35], line[36], line[37], line[38], line[39], line[40], line[41], 
      line[42], line[43], line[44], line[45], line[46], line[47], line[48], 
      line[49], line[50], line[51], line[52], line[53], line[54], line[55], 
      line[56], line[57], line[58], line[59], line[60], line[61], line[62], 
      line[63], PlotRange -> {{0, -17000}, {minheight, maxheight}}, 
      PlotJoined -> True, AxesOrigin -> {0, minheight}, 
      Axes -> {True, True}, SymbolShape -> None, 
      PlotStyle -> GrayLevel[0], 
      Ticks -> {Automatic, None}]

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

Clear[line]
For[mm = 1; tabjoined = Join[tabH[4], tabV[4]], 
        mm ≤ 
          Length[Join[tabH[1], tabV[1]]], mm++, 
line[
            mm] = tabjoined[[mm]]];
minheight = Min[height[4, 10]] - 20; 
    maxheight = Max[height[4, 10]] + 10;

MultipleListPlot[line[1], line[2], line[3], line[4], line[5], line[6], 
      line[7], line[8], line[9], line[10], line[11], line[12], line[13], 
      line[14], line[15], line[16], line[17], line[18], line[19], line[20], 
      line[21], line[22], line[23], line[24], line[25], line[26], line[27], 
      line[28], line[29], line[30], line[31], line[32], line[33], line[34], 
      line[35], line[36], line[37], line[38], line[39], line[40], line[41], 
      line[42], line[43], line[44], line[45], line[46], line[47], line[48], 
      line[49], line[50], line[51], line[52], line[53], line[54], line[55], 
      line[56], line[57], line[58], line[59], line[60], line[61], line[62], 
      line[63], PlotRange -> {{0, -17000}, {minheight, maxheight}}, 
      PlotJoined -> True, AxesOrigin -> {0, minheight}, 
      Axes -> {True, True}, SymbolShape -> None, 
      PlotStyle -> GrayLevel[0], 
      Ticks -> {Automatic, None}]

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