Generating the Figures in Chapter 14 of Otto and Day

NOTE:  Some of the results involving random numbers differ slightly from those reported in the book because of updates in Mathematica between version 5.0 (used to write the book) and version 8.0 (here).

Clear["Global`*"]

Figure 14.2 - Five state model simulation

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

initvec={0,0,1,0,0};

SeedRandom[4323]

For[trial=1,trial<11,trial++,For[t = 1;vec[trial,0]=initvec,t< 101,t++,
newvec=(mat/.a->0.2/.b->0.35/.g->0.12).vec[trial,t-1];vec[trial,t]=RandomInteger[MultinomialDistribution[1,newvec]]
]
]

Show[ListPlot[Table[{t,{1,2,3,4,5}.vec[1,t]},{t,0,30}],Joined->True],ListPlot[Table[{t,{1,2,3,4,5}.vec[2,t]},{t,0,30}],Joined->True,PlotStyle->Dashing[{0.005,0.01}]],ListPlot[Table[{t,{1,2,3,4,5}.vec[3,t]},{t,0,30}],Joined->True,PlotStyle->Thickness[0.01]],ListPlot[Table[{t,{1,2,3,4,5}.vec[4,t]},{t,0,30}],Joined->True,PlotStyle->Dashing[{0.04}]],ListPlot[Table[{t,{1,2,3,4,5}.vec[1,t]},{t,0,30}],PlotStyle->PointSize[0.015]],ListPlot[Table[{t,{1,2,3,4,5}.vec[2,t]},{t,0,30}],PlotStyle->PointSize[0.015]],ListPlot[Table[{t,{1,2,3,4,5}.vec[3,t]},{t,0,30}],PlotStyle->PointSize[0.015]],ListPlot[Table[{t,{1,2,3,4,5}.vec[4,t]},{t,0,30}],PlotStyle->PointSize[0.015]],
PlotRange->{Automatic,{1,5}},
AxesOrigin->{0,1}]

Figure 14.3 - Five state model prob distribution

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

initvec={0,0,1,0,0};

MatrixPower[mat,2].initvec

BarChart[initvec/.a->0.2/.b->0.35/.g->0.12,PlotRange->{{0,5.5},{-0.08,1}},AspectRatio->0.4,ChartLabels->Table[i,{i,1,5}],Ticks->{False,Automatic}]

BarChart[MatrixPower[mat,10].initvec/.a->0.2/.b->0.35/.g->0.12,PlotRange->{{0,5.5},{-0.08,1}},AspectRatio->0.4,ChartLabels->Table[i,{i,1,5}],Ticks->{False,Automatic}]

BarChart[MatrixPower[mat,30].initvec/.a->0.2/.b->0.35/.g->0.12,PlotRange->{{0,5.5},{-0.08,1}},AspectRatio->0.4,ChartLabels->Table[i,{i,1,5}],Ticks->{False,Automatic}]

Figure 14.6 - Time to absorption in five state model

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

tvec={t1,t2,t3,t4,t5};

Here are two alternative ways to solve for the time to absorption:

Solve[{tvec==tvec.mat+{0,1,1,1,0},t1==0,t5==0},tvec]

or alternatively...

mat3=Transpose[Take[Transpose[Take[mat,{2,4}]],{2,4}]]

time={1,1,1}.Inverse[IdentityMatrix[3]-mat3]//Simplify

The probability of absorption in state 1 equals (from the left eigenvector):

fixvec1={1,((2 a+b)/(2 (a+b))),(1/2),(b/(2 (a+b))),0};

umat=Table[fixvec1[[i]]/fixvec1[[j]],{i,2,4},{j,2,4}]

the conditional time to absorption in state 1:

condtime1={1,1,1}.(umat*Inverse[IdentityMatrix[3]-mat3])//Factor

The probability of absorption in state 5 equals (from the left eigenvector):

fixvec5={0,(b/(2 (a+b))),(1/2),((2 a+b)/(2 (a+b))),1};

umat=Table[fixvec5[[i]]/fixvec5[[j]],{i,2,4},{j,2,4}]

the conditional time to absorption in state 5:

condtime5={1,1,1}.(umat*Inverse[IdentityMatrix[3]-mat3])//Factor

Factor[condtime1-time]

Factor[condtime5-time]

This correctly gives the overall time to absorption:

Take[fixvec1,{2,4}]*condtime1+Take[fixvec5,{2,4}]*condtime5//Factor

%-time

Starting from state 4, the following illustrates the unconditional absorption time (solid thin), the conditional time to absorption in state 1 (thick), and the conditional time to absorption in state 5 (dashed), as a function of gamma:

Show[
Plot[time[[3]]/.a->0.2/.b->0.2,{g,0,0.5},PlotRange->All],
Plot[condtime1[[3]]/.a->0.2/.b->0.2,{g,0,0.5},PlotStyle->Thickness[0.015],PlotRange->All],
Plot[condtime5[[3]]/.a->0.2/.b->0.2,{g,0,0.5},PlotStyle->Dashing[{0.02}],PlotRange->All],
PlotRange->{{0,0.5},{0,35}}
]

Starting from state 4, the following illustrates the unconditional absorption time (solid thin), the conditional time to absorption in state 1 (thick), and the conditional time to absorption in state 5 (dashed), as a function of alpha:

Show[
Plot[time[[3]]/.g->0.25/.b->0.2,{a,0,0.5},PlotRange->All],
Plot[condtime1[[3]]/.g->0.25/.b->0.2,{a,0,0.5},PlotStyle->Thickness[0.015],PlotRange->All],
Plot[condtime5[[3]]/.g->0.25/.b->0.2,{a,0,0.5},PlotStyle->Dashing[{0.02}],PlotRange->All],
PlotRange->{{0,0.5},{0,35}}
]

Figure 14.7 - Fixation probability in Moran model

Simplify[Sum[(w)^i,{i,0,x-1}]/Sum[(w)^i,{i,0,N-1}]]

Show[
Plot[((-1+(w^x))/(-1+(w^N)))/.w->1/(1+s)/.x->1/.N->10,{s,-0.1,0.25}],
Plot[((-1+(w^x))/(-1+(w^N)))/.w->1/(1+s)/.x->1/.N->100,{s,-0.1,0.25}],
Plot[((-1+(w^x))/(-1+(w^N)))/.w->1/(1+s)/.x->1/.N->1000,{s,-0.1,0.25}],
PlotRange->{{-0.1,0.2},{0,0.2}},
AxesOrigin->{0,0}]

Here we check whether the two exact ways to calculate the fixation probability coincide:

num=100;

fix100=Table[{s,(Table[If[num==i+1,p'(1-p),0]/.p'->(1+s) p/(1+s p)/.p->i/num,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-Table[If[j==i,1-p'(1-p)-(1- p')p,If[j==i+1,p'(1-p),If[j==i-1,(1-p')p,0]]]/.p'->(1+s) p/(1+s p)/.p->i/num//N,{j,1,num-1},{i,1,num-1}]])[[1]]},{s,-0.1,0.1,0.2}]

Table[{s,((-1+(w^x))/(-1+(w^N)))/.w->1/(1+s)/.x->1/.N->num},{s,-0.1,0.1,0.2}]

Figure 14.8 - Time to absorption in the simple birth-death model (absorbing states and 0 and N)

Here we explore a birth-death model where we take N to equal the maximum population size, at which point the process is killed.  This might represent a threshold population size, above which we are less concerned about the viability of a population.

birth=b*i;
death=d*i;

The mean time until absorption for a branching process:

gamma[x_]=(d/b)^x;

Clear[t]

t[1,N_,b_,d_]:=t[1,N,b,d]=Sum[Sum[(d/b)^(i-j)/(b j),{j,1,i}],{i,0,N-1}]/Sum[(d/b)^(i),{i,0,N-1}];

t[x_,N_,b_,d_]:=t[x,N,b,d]=Sum[Sum[(d/b)^(i-j)/(b j),{j,1,i}],{i,x,N-1}]-t[1,N,b,d]Sum[(d/b)^(i),{i,x,N-1}]

num=100;

If birth and death rates are per day (e.g., a birth is expected each 500 days), then the following measures absorption times in years:

Show[ListPlot[Table[{x,(1/365) t[x,num,0.002,0.0019]},{x,0,num}],Joined->True,PlotStyle->Thickness[0.01]],ListPlot[Table[{x,N[(1/365) t[x,num,(2/1000),(38/10000)],30]},{x,0,num}],Joined->True],ListPlot[Table[{x,(1/365) t[x,num,0.004,0.0019]},{x,0,num}],Joined->True],ListPlot[Table[{x,(1/365) t[x,num,0.004,0.0038]},{x,0,num}],Joined->True,PlotStyle->Dashing[{0.01}]],PlotRange->{Automatic,{0,40}}]

With a smaller value of d:

Show[ListPlot[Table[{x,(1/365) t[x,num,0.002,0.0009]},{x,0,num}],Joined->True,PlotStyle->Thickness[0.01]],ListPlot[Table[{x,N[(1/365) t[x,num,(2/1000),(18/10000)],30]},{x,0,num}],Joined->True],ListPlot[Table[{x,(1/365) t[x,num,0.004,0.0009]},{x,0,num}],Joined->True],ListPlot[Table[{x,(1/365) t[x,num,0.004,0.0018]},{x,0,num}],Joined->True],PlotRange->{Automatic,{0,40}}]

Figure 14.9 - Mean number of cases

Show[
Plot[mu^t/.mu->0.9,{t,0,50}],
Plot[ mu^t+2Sqrt[var mu^(t-1) (1-mu^t)/(1-mu)]/.var->0.9/.mu->0.9,{t,0,50},PlotStyle->Dashing[{0.02}]],
Plot[ mu^t+2Sqrt[var mu^(t-1) (1-mu^t)/(1-mu)]/.var->4.5/.mu->0.9,{t,0,50},PlotStyle->Dashing[{0.04}]],
PlotRange->All
]