Generating the Figures in Chapter 8 of Otto and Day

Clear["Global`*"]

Figure 8.1 and 8.2

Clear["Global`*"]

f[s_,i_]:=θ-d s-β s i+γ i;
g[s_,i_]:=β s i-(d+v+γ) i;

Clear[v,d,γ,β,θ]

equil=Simplify[Solve[{f[s,i]==0,g[s,i]==0},{s,i}]]

endemic={s->((d+v+γ)/β),i->((-((d (d+v+γ))/β)+θ)/(d+v))};

m1=Simplify[{{D[f[s,i],s],D[f[s,i],i]},{D[g[s,i],s],D[g[s,i],i]}}/.endemic]

m2=Simplify[{{D[f[s,i],s],D[f[s,i],i]},{D[g[s,i],s],D[g[s,i],i]}}/.{s->(θ/d),i->0}]

v=.7;d=.1;γ=0.1;β=.005;θ=50;
equil
β θ/d-(d+v+γ)

Show[ContourPlot[{f[s,i]==0,g[s,i]==0},{s,0,550},{i,0,60},AspectRatio->1,PlotPoints->200,ContourStyle->Thickness[0.012]],ListPlot[{{165,38},{510,2}},PlotMarkers->{★,15}],PlotRange->{{0,550},{-2,60}}]

Clear[DifSI,s0,i0]
DifSI[s0_,i0_,tmin_,tmax_]:=DifSI[s0,i0,tmin,tmax]=NDSolve[{D[s[t],t]==θ-d s[t]-β s[t] i[t]+γ i[t],D[i[t],t]==β s[t] i[t]-(d+v+γ) i[t],s[tmin]==s0,i[tmin]==i0},{s[t],i[t]},{t,tmin,tmax}]

Using the linear approximation to the dynamics near the endemic equilibrium (note that the following describe the displacements from the equilibrium):

Clear[DifSI1,es0,ei0,es,is]
DifSI1[es0_,ei0_,tmin_,tmax_]:=DifSI1[es0,ei0,tmin,tmax]=NDSolve[{D[es[t],t]==m1[[1,1]] es[t]+m1[[1,2]] ei[t],D[ei[t],t]==m1[[2,1]] es[t] + m1[[2,2]] ei[t],es[tmin]==es0,ei[tmin]==ei0},{es[t],ei[t]},{t,tmin,tmax}]

Using the linear approximation to the dynamics near the disease-absent equilibrium (note that the following describe the displacements from the equilibrium):

Clear[DifSI2,es0,ei0,es,is]
DifSI2[es0_,ei0_,tmin_,tmax_]:=DifSI2[es0,ei0,tmin,tmax]=NDSolve[{D[es[t],t]==m2[[1,1]] es[t]+m2[[1,2]] ei[t],D[ei[t],t]==m2[[2,1]] es[t] + m2[[2,2]] ei[t],es[tmin]==es0,ei[tmin]==ei0},{es[t],ei[t]},{t,tmin,tmax}]

Starting at {165,38}:

Show[ContourPlot[{f[s,i]==0,g[s,i]==0},{s,0,550},{i,0,60},AspectRatio->1,PlotPoints->200,ContourStyle->Thickness[0.025]],ParametricPlot[Evaluate[{s[t],i[t]}/.DifSI[165,38,0,100]],{t,0,100}],ParametricPlot[Evaluate[{(s/.endemic)+es[t],(i/.endemic)+ei[t]}/.DifSI1[165-(s/.endemic),38-(i/.endemic),0,100]],{t,0,100},PlotStyle->Dashing[{0.02}]],ListPlot[{{165,38}},PlotMarkers->{★,15}],PlotRange->{{160,190},{32,47}},AxesOrigin->{160,32}]

Show[ContourPlot[{f[s,i]==0,g[s,i]==0},{s,0,550},{i,0,60},AspectRatio->1,PlotPoints->200,ContourStyle->Thickness[0.025]],ParametricPlot[Evaluate[{s[t],i[t]}/.DifSI[505,1,0,100]],{t,0,100}],ParametricPlot[Evaluate[{(s/.endemic)+es[t],(i/.endemic)+ei[t]}/.DifSI1[505-(s/.endemic),1-(i/.endemic),0,100]],{t,0,100},PlotStyle->Dashing[{0.02}]],ListPlot[{{505,1}},PlotMarkers->{★,15}],PlotRange->{{460,519},{-0.2,4.8}},AxesOrigin->{460,0}]

Figure 8.3

Clear["Global`*"]

f[n1_,n2_]:=θ-a c n1 n2
g[n1_,n2_]:=ϵ a c n1 n2-d n2

Simplify[Solve[{f[n1,n2]==0,g[n1,n2]==0},{n1,n2}]]

m=Simplify[{{D[f[n1,n2],n1],D[f[n1,n2],n2]},{D[g[n1,n2],n1],D[g[n1,n2],n2]}}/.{n2->((ϵ θ)/d),n1->(d/(a c ϵ))}]

Det[m]

Simplify[Eigenvalues[m]]

θ=10;d=2;a=1;c=.1;ϵ=1;tmax=40;
NDSolve[{n1'[t]==θ-a c n1[t] n2[t],n2'[t]==ϵ a c n1[t] n2[t]-d n2[t],n1[0]==1,n2[0]==1},{n1[t],n2[t]},{t,0,tmax}][[1]];
Plot[{n1[t]/.%,n2[t]/.%},{t,0,tmax},PlotStyle->{Dashing[{0}],Dashing[{.02}]},PlotRange->{0,40}]
{-((a c ϵ θ+(Sqrt[a]) (Sqrt[c]) (Sqrt[ϵ]) (Sqrt[θ]) (Sqrt[-4 (d^2)+a c ϵ θ]))/(2 d)),((-a c ϵ θ+(Sqrt[a]) (Sqrt[c]) (Sqrt[ϵ]) (Sqrt[θ]) (Sqrt[-4 (d^2)+a c ϵ θ]))/(2 d))}
ϵ a c θ-4 (d^2)

4 Pi d/Sqrt[a c ϵ (4 d^2 - a c ϵ θ)]

result=NDSolve[{n1'[t]==θ-a c n1[t] n2[t],n2'[t]==ϵ a c n1[t] n2[t]-d n2[t],n1[0]==1,n2[0]==1},{n1[t],n2[t]},{t,0,tmax}][[1]]

plot1=ParametricPlot[Evaluate[{n1[t],n2[t]}/.result],{t,0,tmax},AspectRatio->Automatic,PlotRange->{{0,40},{0,40}},PlotStyle->Thickness[.01]]

plot2=(VectorPlot[{θ-a c x y,ϵ a c x y-d y},{x,0,40},{y,0,40},Axes->True,VectorScale->{Medium,0.5,#5&},AspectRatio->1,PlotRange->{{0,40},{0,40}}])

plot3=ContourPlot[{θ-a c x y==0,ϵ a c x y-d y==0},{x,0,40},{y,0,40},PlotPoints->300]

Show[plot1,plot2,plot3]

Export["s.EPS",%];

θ=10;d=.3;a=1;c=.1;ϵ=1;tmax=40;
NDSolve[{n1'[t]==θ-a c n1[t] n2[t],n2'[t]==ϵ a c n1[t] n2[t]-d n2[t],n1[0]==1,n2[0]==1},{n1[t],n2[t]},{t,0,tmax}][[1]];
Plot[{n1[t]/.%,n2[t]/.%},{t,0,tmax},PlotStyle->{Dashing[{0}],Dashing[{.02}]},PlotRange->{0,40}]
{-((a c ϵ θ+(Sqrt[a]) (Sqrt[c]) (Sqrt[ϵ]) (Sqrt[θ]) (Sqrt[-4 (d^2)+a c ϵ θ]))/(2 d)),((-a c ϵ θ+(Sqrt[a]) (Sqrt[c]) (Sqrt[ϵ]) (Sqrt[θ]) (Sqrt[-4 (d^2)+a c ϵ θ]))/(2 d))}
ϵ a c θ-4 (d^2)

result=NDSolve[{n1'[t]==θ-a c n1[t] n2[t],n2'[t]==ϵ a c n1[t] n2[t]-d n2[t],n1[0]==1,n2[0]==1},{n1[t],n2[t]},{t,0,tmax}][[1]]

plot4=ParametricPlot[Evaluate[{n1[t],n2[t]}/.result],{t,0,tmax},AspectRatio->Automatic,PlotRange->{{0,40},{0,40}},PlotStyle->Thickness[.01]]

plot5=(VectorPlot[{θ-a c x y,ϵ a c x y-d y},{x,0,40},{y,0,40},Axes->True,VectorScale->{Medium,0.5,#5&},AspectRatio->1,PlotRange->{{0,40},{0,40}}])

plot6=ContourPlot[{θ-a c x y==0,ϵ a c x y-d y==0},{x,0,40},{y,0,40},PlotPoints->300]

Show[plot4,plot5,plot6]

Export["ss.EPS",%];

{p->(θ/μ),n->(μ/s)}
{-((s θ+(Sqrt[s]) (Sqrt[θ]) (Sqrt[s θ-4 (μ^2)]))/(2 μ)),((-s θ+(Sqrt[s]) (Sqrt[θ]) (Sqrt[s θ-4 (μ^2)]))/(2 μ))}

Figure 8.4

Clear["Global`*"]

critβ=δ μ σ(α+μ)/(Γ πP τ (α+(1-p) μ))

Here is where we plug in the parameter values from Phillips' (1996) Table 1.

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

Show[
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->fac*0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->fac*3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10},PlotStyle->Dashing[{0.01}]],
PlotRange->{Automatic,{0.00003,0.00004}},AxesOrigin->{-1,0.00003}]

Show[
LogLinearPlot[critβ/.{Γ->1.36,μ->fac*1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->fac*0.33,πP->100},{fac,0.1,10}],
PlotRange->{Automatic,{0,0.0001}},AxesOrigin->{-1,0}]

Show[
LogLinearPlot[critβ/.{Γ->fac*1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->fac*100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->fac*0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
PlotRange->{Automatic,{0,0.0001}},AxesOrigin->{-1,0}]

Show[
LogLinearPlot[critβ/.{Γ->fac*1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->fac*1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10},PlotStyle->Dashing[{0.01}]],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->fac*0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->fac*0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->fac*3.6*10^(-2),σ->2,δ->0.33,πP->100},{fac,0.1,10},PlotStyle->Dashing[{0.01}]],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->fac*2,δ->0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->fac*0.33,πP->100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->2,δ->0.33,πP->fac*100},{fac,0.1,10}],
LogLinearPlot[critβ/.{Γ->1.36,μ->1.36*10^(-3),τ->0.2,
		 p->0.1,α->3.6*10^(-2),σ->fac*2,δ->0.33,πP->100},{fac,0.1,10}],
PlotRange->{Automatic,{0,0.00035}},AxesOrigin->{-1,0}]

Figure 8.5

Clear["Global`*"]

r=0.5;
Plot[{(1+r) Exp[-.01 n],(1+r) Exp[-.1 n]},{n,0,100},PlotRange->{0,2}]

Export["graph.EPS",%];

Figure 8.6 - Density-dependent natural selection

Clear["Global`*"]

w[n_,r_,α_]:=(1+r) Exp[-α n]

r=29;
(1-Log[1+r])//N

n[1]=10;p[1]:=.9;
αAA=.01;αAa=.015;αaa=.02;
r=4;
{1+r,Exp[2]}//N
{αAA,αAa}//N
tmax=100;
For[t=1,t≤tmax,t++,
n[t+1]=n[t] (((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(1-p[t])^2 w[n[t],r,αaa]);p[t+1]=p[t]*((p[t] w[n[t],r,αAA]+(1-p[t]) w[n[t],r,αAa])/(((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(((1-p[t]))^2) w[n[t],r,αaa]))]

plot1=Show[
ListPlot[Table[n[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010}],ListPlot[Table[1000 p[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010},PlotStyle->Dashing[{0.015}]],
Frame->True]

Export["s1.EPS",%];

n[1]=10;p[1]:=.9;
αAA=.01;αAa=.015;αaa=.02;
r=7;
{1+r,Exp[2]}//N
{αAA,αAa}//N
tmax=100;
For[t=1,t≤tmax,t++,
n[t+1]=n[t] (((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(1-p[t])^2 w[n[t],r,αaa]);p[t+1]=p[t]*((p[t] w[n[t],r,αAA]+(1-p[t]) w[n[t],r,αAa])/(((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(((1-p[t]))^2) w[n[t],r,αaa]))]

plot2=Show[
ListPlot[Table[n[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010}],ListPlot[Table[1000 p[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010},PlotStyle->Dashing[{0.015}]],
Frame->True]

Export["s2.EPS",%];

n[1]=10;p[1]:=.99;
αAA=.0155;αAa=.015;αaa=.02;
r=4;
{1+r,Exp[2]}//N
{αAA,αAa}//N
tmax=100;
For[t=1,t≤tmax,t++,
n[t+1]=n[t] (((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(1-p[t])^2 w[n[t],r,αaa]);p[t+1]=p[t]*((p[t] w[n[t],r,αAA]+(1-p[t]) w[n[t],r,αAa])/(((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(((1-p[t]))^2) w[n[t],r,αaa]))]

plot3=Show[
ListPlot[Table[n[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010}],ListPlot[Table[1000 p[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010},PlotStyle->Dashing[{0.015}]],
Frame->True]

Export["s3.EPS",%];

n[1]=10;p[1]:=.99;
αAA=.0155;αAa=.015;αaa=.02;
r=7;
{1+r,Exp[2]}//N
{αAA,αAa}//N
tmax=100;
For[t=1,t≤tmax,t++,
n[t+1]=n[t] (((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(1-p[t])^2 w[n[t],r,αaa]);p[t+1]=p[t]*((p[t] w[n[t],r,αAA]+(1-p[t]) w[n[t],r,αAa])/(((p[t])^2) w[n[t],r,αAA]+2 p[t] (1-p[t]) w[n[t],r,αAa]+(((1-p[t]))^2) w[n[t],r,αaa]))]

plot4=Show[ListPlot[Table[n[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010}],
ListPlot[Table[1000 p[i],{i,1,tmax}],Joined->True,PlotRange->{0,1010},PlotStyle->Dashing[{0.015}]],
Frame->True]

Export["s4.EPS",%];

Figure 8.7 -  Two-locus

Clear["Global`*"]

wbar=Sum[x[i] x[j]w[i,j],{i,1,4},{j,1,4}] ;
dis=x[1] x[4] (w[1,4]+w[4,1])/(2 wbar)-x[2] x[3] (w[2,3]+w[3,2])/(2 wbar);
eqn1=x[1] Sum[x[j] (w[1,j]+w[j,1])/(2 wbar),{j,1,4}] - r dis;
eqn2=x[2] Sum[x[j] (w[2,j]+w[j,2])/(2 wbar),{j,1,4}] + r dis;
eqn3=x[3] Sum[x[j] (w[3,j]+w[j,3])/(2 wbar),{j,1,4}] + r dis;
eqn4=x[4] Sum[x[j] (w[4,j]+w[j,4])/(2 wbar),{j,1,4}] - r dis;

mat={{D[eqn1,x[1]],D[eqn1,x[2]],D[eqn1,x[3]],D[eqn1,x[4]]},
{D[eqn2,x[1]],D[eqn2,x[2]],D[eqn2,x[3]],D[eqn2,x[4]]},
{D[eqn3,x[1]],D[eqn3,x[2]],D[eqn3,x[3]],D[eqn3,x[4]]},
{D[eqn4,x[1]],D[eqn4,x[2]],D[eqn4,x[3]],D[eqn4,x[4]]}}/.{x[1]->1,x[2]->0,x[3]->0,x[4]->0}//Simplify;
MatrixForm[mat]

Histogram3D[Join[Table[{0,0},{i,1,10}],Table[{1,0},{i,1,8}],Table[{0,1},{i,1,8}],Table[{1,1},{i,1,12}]],{{-0.1,0.1,0.3,0.5,0.7,0.9,1.1}},Ticks->{None,None,{0,2,4,6,8,10,12}},BoxRatios->{1,1,1}]

Export["graph.EPS",%];

Figure 8.12 - Evolution of haploidy/diploidy

Clear["Global`*"]

x1=frequency(A C1)
x2=frequency(a C1)
x3=frequency(A C2)
x4=frequency(a C2)

where x1+x2+x3+x4=1.

Here we describe the recursions for the more general model, where dd gives the probability that an organism is diploid at the time of selection and 1-dd is the probability that it is haploid.  In the book, we assume that dd11 = 0, dd12 = 1, and dd22 = 1.

eqn1=(1-mu)*(x1*x1+x1*x2*(1-h*s*dd11)+x1*x3+(1-r)*x1*x4*(1-h*s*dd12)+r*x2*x3*(1-h*s*dd12));
eqn2=((mu)*(x1*x1+x1*x2*(1-h*s*dd11)+x1*x3+(1-r)*x1*x4*(1-h*s*dd12)+r*x2*x3*(1-h*s*dd12))
+(x1*x2*(1-s+s*dd11-h*s*dd11)+r*x1*x4*(1-s+s*dd12-h*s*dd12)+x2*x2*(1-s)+(1-r)*x2*x3*(1-s+s*dd12-h*s*dd12)+x2*x4*(1-s)));
eqn3=(1-mu)*(x1*x3+x1*x4*r*(1-h*s*dd12)+x2*x3*(1-r)*(1-h*s*dd12)+x3*x3+x3*x4*(1-h*s*dd22));
eqn4=((mu)*((x1*x3+x1*x4*r*(1-h*s*dd12)+x2*x3*(1-r)*(1-h*s*dd12)+x3*x3+x3*x4*(1-h*s*dd22)))
+(x1*x4*(1-r)*(1-h*s*dd12-s+s*dd12)+x2*x3*(1-h*s*dd12-s+dd12*s)*r+x2*x4*(1-s)+x3*x4*(1-s+s*dd22-h*s*dd22)+x4*x4*(1-s)));

total=eqn1+eqn2+eqn3+eqn4;

The recursions are thus given by eqnX/total.

Factor[eqn1-x1*total/.x3->0/.x4->0/.x1->1-x2/.dd11->0]

Focusing on the case where the initial population is haploid allows us to simplify the original equilibrium to x2 = mu/s.

subx={x1->x1hat+eps1*sm,x2->x2hat+eps2*sm,x3->eps3*sm,x4->eps4*sm};

mat={{Coefficient[Series[eqn1/total/.subx,{sm,0,1}],eps1],
Coefficient[Series[eqn1/total/.subx,{sm,0,1}],eps2],
Coefficient[Series[eqn1/total/.subx,{sm,0,1}],eps3],
Coefficient[Series[eqn1/total/.subx,{sm,0,1}],eps4]},
{Coefficient[Series[eqn2/total/.subx,{sm,0,1}],eps1],
Coefficient[Series[eqn2/total/.subx,{sm,0,1}],eps2],
Coefficient[Series[eqn2/total/.subx,{sm,0,1}],eps3],
Coefficient[Series[eqn2/total/.subx,{sm,0,1}],eps4]},
{Coefficient[Series[eqn3/total/.subx,{sm,0,1}],eps1],
Coefficient[Series[eqn3/total/.subx,{sm,0,1}],eps2],
Coefficient[Series[eqn3/total/.subx,{sm,0,1}],eps3],
Coefficient[Series[eqn3/total/.subx,{sm,0,1}],eps4]},
{Coefficient[Series[eqn4/total/.subx,{sm,0,1}],eps1],
Coefficient[Series[eqn4/total/.subx,{sm,0,1}],eps2],
Coefficient[Series[eqn4/total/.subx,{sm,0,1}],eps3],
Coefficient[Series[eqn4/total/.subx,{sm,0,1}],eps4]}}/.sm->1;

mat=Simplify[mat (x1hat+x2hat)/.dd11->0]/.(x1hat+x2hat)->1;

MatrixForm[mat]

To simplify matters,

MatrixForm[Factor[mat/.dd12->1/.x1hat->1-x2hat/.x2hat->mu/s]]

mat1=Transpose[Drop[Transpose[Drop[Factor[mat/.x1hat->1-x2hat/.x2hat->mu/s],-2]],-2]];

MatrixForm[mat1]

Eigenvalues[mat1]

mat2=Transpose[Drop[Transpose[Drop[Factor[mat/.x1hat->1-x2hat/.x2hat->mu/s],2]],2]];
MatrixForm[mat2]

When dd12 = 1, this can be rewritten as:

MatrixForm[mat2/.dd12->1]

mat2b={{1-((mu (r+h s-h r s))/s),((r (s-mu) (1-h s))/s)},
{((mu (r + s -h r s)-mu^2(r + h s -h r s))/((1-mu) s)),1-(((s-mu) (h s+(1-mu) r (1-h s)))/((1-mu) s))}};

MatrixForm[mat2b]

MatrixForm[Factor[mat2b/.mu->0]]

Simplify[Eigenvalues[%]]

Factor[Det[mat2-lam*IdentityMatrix[2]]]

char=Factor[Det[mat2-lam*IdentityMatrix[2]]]

Factor[Normal[Series[char/.lam->oone+omu*mu,{mu,0,0}]]]

Thus, one eigenvalue is near one (oone = 1).  To order mu, this eigenvalue is:

Normal[Series[char/.lam->1+omu*mu,{mu,0,1}]]

The following indicates the critical value of the eigenvalue (where lam = 1+omu*mu) to leading order in mu:

criterion=Simplify[omu/.Solve[%==0,omu]][[1]]

For the case considered in the book, where the modifier causes a fully diploid life cycle in heterozygotes:

Simplify[criterion/.dd12->1]

which can be rewritten as (( r (1-2 h)(1-h s)- (h^2) s)/(h s+r (1 - h s)))

Show[
Plot[mu(( r (1-2 h)(1-h s)- (h^2) s)/(h s+r (1 - h s)))/.mu->10^(-6)/.s->0.1/.r->0.5,{h,0,1}],
Plot[mu(( r (1-2 h)(1-h s)- (h^2) s)/(h s+r (1 - h s)))/.mu->10^(-6)/.s->0.1/.r->0.05,{h,0,1}],
Plot[mu(( r (1-2 h)(1-h s)- (h^2) s)/(h s+r (1 - h s)))/.mu->10^(-6)/.s->0.1/.r->0.005,{h,0,1}],
Plot[Eigenvalues[mat2b][[2]]-1/.mu->10^(-6)/.s->0.1/.r->0.5,{h,0,1}],
Plot[Eigenvalues[mat2b][[2]]-1/.mu->10^(-6)/.s->0.1/.r->0.05,{h,0,1}],
Plot[Eigenvalues[mat2b][[2]]-1/.mu->10^(-6)/.s->0.1/.r->0.005,{h,0,1}]
]

Export["graph.EPS",%];

The criterion for invasion of diploidy (with dd12=1):

hsoln=h/.Solve[(r (1-2 h)(1-h s)- (h^2) s)==0,h][[1]];
Show[
Plot[hsoln/.r->0.01,{s,0,1},PlotRange->All],
Plot[hsoln/.r->0.1,{s,0,1},PlotRange->All],
Plot[hsoln/.r->0.4999,{s,0,1},PlotRange->All],
PlotRange->{{0,1},{0,1}},
AxesOrigin->{0,0}
]

It's harder for a little bit of diploidy to invade:

Show[
Plot[Simplify[h/.Solve[(criterion/.r->0.01/.dd12->0.01)==0,h][[2]]],{s,0,1},PlotRange->All,PlotStyle->Dashed],
Plot[Simplify[h/.Solve[(criterion/.r->0.1/.dd12->0.01)==0,h][[2]]],{s,0,1},PlotRange->All,PlotStyle->Dashed],
Plot[Simplify[h/.Solve[(criterion/.r->0.499/.dd12->0.01)==0,h][[2]]],{s,0,1},PlotRange->All,PlotStyle->Dashed],
PlotRange->{{0,1},{0,1}},
AxesOrigin->{0,0}
]

Figure 8.1.1 - Taylor series

Show[
Plot[Exp[x]Sin[y]/.y->0.025,{x,0,0.101}],
Plot[Exp[x]Sin[y]/.y->0.05,{x,0,0.101}],
Plot[Exp[x]Sin[y]/.y->0.075,{x,0,0.101}],
Plot[y/.y->0.025,{x,0,0.101},PlotStyle->Dashing[{0.01}]],
Plot[y/.y->0.05,{x,0,0.101},PlotStyle->Dashing[{0.01}]],
Plot[y/.y->0.075,{x,0,0.101},PlotStyle->Dashing[{0.01}]],
PlotRange->{{0,0.1},{0,0.1}},
AxesOrigin->{0,0}
]

Show[
Plot[Exp[x]Sin[y]/.x->0.025,{y,0,0.101}],
Plot[Exp[x]Sin[y]/.x->0.05,{y,0,0.101}],
Plot[Exp[x]Sin[y]/.x->0.075,{y,0,0.101}],
Plot[y/.x->0.025,{y,0,0.101},PlotStyle->Dashing[{0.01}]],
Plot[y/.x->0.05,{y,0,0.101},PlotStyle->Dashing[{0.01}]],
Plot[y/.x->0.075,{y,0,0.101},PlotStyle->Dashing[{0.01}]],
PlotRange->{{0,0.1},{0,0.1}},
AxesOrigin->{0,0}
]

Plot3D[Exp[x]Sin[y],{x,0,0.101},{y,0,0.101},PlotRange->{Automatic,Automatic,{0,0.11}}]

Export["test.eps",%];

Sup Mat Figure 8.1.1 - Parabolas

parabola [lam_,b_,c_]:= lam^2 - b lam +c

Show[
Plot[parabola[lam,1,3/4],{lam,-4,4},PlotStyle->Thickness[0.011]],
Plot[parabola[lam,2,6/4],{lam,-4,4}],
PlotRange->{{-2.5,2.5},{-1,2}}]

Export["graph.EPS",%];

Solve[parabola[lam,-1,3/4]==0,lam]
Abs[lam/.%]//N

Solve[parabola[lam,2,6/4]==0,lam]
Abs[lam/.%]//N

Show[
Plot[parabola[lam,-1,-1/2],{lam,-4,4}],
Plot[parabola[lam,-1,1/8],{lam,-4,4},PlotStyle->Thickness[0.011]],
Plot[parabola[lam,1,-1/2],{lam,-4,4}],
Plot[parabola[lam,1,1/8],{lam,-4,4},PlotStyle->Thickness[0.011]],
PlotRange->{{-2.5,2.5},{-1,2}}
]

Export["graph.EPS",%];

Solve[parabola[lam,-1,-1/2]==0,lam]//N
Solve[parabola[lam,-1,1/8]==0,lam]//N
Solve[parabola[lam,1,-1/2]==0,lam]//N
Solve[parabola[lam,1,1/8]==0,lam]//N

Show[
Plot[parabola[lam,-4,7/2],{lam,-4,4},PlotStyle->Dashing[{0.04}]],
Plot[parabola[lam,4,7/2],{lam,-4,4}],
PlotRange->{{-3.5,3.5},{-1,2}}
]

Export["graph.EPS",%];

Solve[parabola[lam,-4,7/2]==0,lam]//N
Solve[parabola[lam,4,7/2]==0,lam]//N

The following is a graphical "representation" of complex roots (the reflected parabola indicates roots with the same absolute magnitude and real part):

Show[
Plot[parabola[lam,2/3,3/4],{lam,-4,4},PlotStyle->Thickness[0.011]],
Plot[(2(-((b^2)/4)+c)/.b->2/3/.c->3/4)-parabola[lam,2/3,3/4],{lam,-4,4},PlotStyle->Dashing[{0.01}]],
PlotRange->{{-3,3},{-1,3}}
]

Solve[2(-((b^2)/4)+c)-parabola[lam,b,c]==0,lam]

Solve[parabola[lam,b,c]==0,lam]

Solve[parabola[lam,b,c]==0,lam]

Abs[lam/.%]

Solve[D[parabola[lam,b,c],lam]==0,lam]

parabola[lam,b,c]/.%

Sup Mat Figure 8.2.1 - Lotka-Volterra conditions

Equilibrium with both species present:

Sub1=Simplify[Solve[{(1+r1*(1-(N1+a12*N2)/K1))==1,
	(1+r2*(1-(N2+a21*N1)/K2))==1},{N1,N2}]][[1]]

Sub2=Simplify[Solve[{(1+r1*(1-(N1+a12*N2)/K1))==1,
	(1+r2*(1-(N2+a21*N1)/K2))==1},{K1,K2}]][[1]]

Sub3=Solve[{N1==(-K1 + a12*K2)/(-1 + a12*a21), N2 ==(a21*K1 - K2)/(-1 + a12*a21)},{a12,a21}]

Throughout the following, it turns out that the conditions depend strongly on
N1/K1 r1 and N2/K2 r2, rather than r1 and r2 themselves.  Let er1 = N1/K1 r1
and er2 = N2/K2 r2.

Sub4= Solve[{er1 == N1/K1 r1,
	er2 == N2/K2 r2},{r1,r2}][[1]]

Sub4b= Solve[({er1 == N1/K1 r1,
	er2 == N2/K2 r2}/.Sub1),{er1,er2}][[1]]

The following plots the region in {r1,r2} space that allows the equilibrium to be stable:

Show[
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.01)==0,r2][[1]],{r1,0,2.02},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.25)==0,r2][[1]],{r1,0,2.6},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.5)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.75)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotRange->All],
PlotRange->{{0,4.1},{0,4.1}},
AspectRatio->1
]

Export["graph.EPS",%];

Note that these restrictions are all based on the requirement that P>0 at lambda = -1 (i.e., that 4-2*er1-2*er2+er1*er2-a12*a21*er1*er2 be positive).  The additional requirement that the derivative be negative at lambda = -1 (-4+er1+er2 < 0) is more permissive in these examples:

Show[
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.01)==0,r2][[1]],{r1,0,2.02},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.5)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.01)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.5)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
PlotRange->{{0,4.1},{0,4.1}}
]

The above picture when K2 = 5 K1:

Show[
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->5 K/.a12->a/.a21->a/.a->0.2)==0,r2][[1]],{r1,0,2.5},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.5)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.7)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.2)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.5)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.7)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->5K/.a12->a/.a21->a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
PlotRange->{{0,4.1},{0,4.1}}
]

The above picture when a21 = a12/2:

Show[
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2  a/.a->0.2)==0,r2][[1]],{r1,0,2.5},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.5)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.7)==0,r2][[1]],{r1,0,3.999},PlotRange->All],
Plot[r2/.Solve[(4-2*er1-2*er2+er1*er2-a12*a21*er1*er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.2)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.5)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.7)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
Plot[r2/.Solve[(-4+er1+er2/.Sub4b/.K1->K/.K2->K/.a12->a/.a21->1/2 a/.a->0.99)==0,r2][[1]],{r1,0,5},PlotStyle->Dashing[{0.01}],PlotRange->All],
PlotRange->{{0,4.1},{0,4.1}}
]