Generating the Figures in Chapter 15 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 15.2 - Probability distribution over time in Moran model

popsize=10;

STARTING IN THE MIDDLE:

ff[x_]:=If[And[x>.45,x<.55],10,0]

sol=NDSolve[{D[n[x,t],t]==(1/2) D[2 x (1-x)/(popsize)^2 n[x,t],x,x],n[x,0]==ff[x]},{n}, {x, 0, 1}, {t, 0,30}]

matrix={{1,0,0,0,0,0,0,0,0,0,0},{p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0,0,0,0},{0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0,0,0},{0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0,0},{0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0},{0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0},{0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0},{0,0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0},{0,0,0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0},{0,0,0,0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p)},{0,0,0,0,0,0,0,0,0,0,1}}/.p->(x/10);

matrix=Transpose[{{1,0,0,0,0,0,0,0,0,0,0},{(1/10) (1-(1/10)) 1,(((1-(1/10)))^2)+((1^2)/100),(1/10) (1-(1/10)) 1,0,0,0,0,0,0,0,0},{0,(1/10) (1-(2/10)) 2,(((1-(2/10)))^2)+((2^2)/100),(1/10) (1-(2/10)) 2,0,0,0,0,0,0,0},{0,0,(1/10) (1-(3/10)) 3,(((1-(3/10)))^2)+((3^2)/100),(1/10) (1-(3/10)) 3,0,0,0,0,0,0},{0,0,0,(1/10) (1-(4/10)) 4,(((1-(4/10)))^2)+((4^2)/100),(1/10) (1-(4/10)) 4,0,0,0,0,0},{0,0,0,0,(1/10) (1-(5/10)) 5,(((1-(5/10)))^2)+((5^2)/100),(1/10) (1-(5/10)) 5,0,0,0,0},{0,0,0,0,0,(1/10) (1-(6/10)) 6,(((1-(6/10)))^2)+((6^2)/100),(1/10) (1-(6/10)) 6,0,0,0},{0,0,0,0,0,0,(1/10) (1-(7/10)) 7,(((1-(7/10)))^2)+((7^2)/100),(1/10) (1-(7/10)) 7,0,0},{0,0,0,0,0,0,0,(1/10) (1-(8/10)) 8,(((1-(8/10)))^2)+((8^2)/100),(1/10) (1-(8/10)) 8,0},{0,0,0,0,0,0,0,0,(1/10) (1-(9/10)) 9,(((1-(9/10)))^2)+((9^2)/100),(1/10) (1-(9/10)) 9},{0,0,0,0,0,0,0,0,0,0,1}}];

z=0;
Show[BarChart[popsize{0,0,0,0,0,1,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=1;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,0,0,0,1,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=3;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,0,0,0,1,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=10;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,0,0,0,1,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=30;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,0,0,0,1,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

STARTING AT p = 0.2:

ff[x_]:=If[And[x>.15,x<.25],10,0]

sol=NDSolve[{D[n[x,t],t]==(1/2) D[2 x (1-x)/(popsize)^2 n[x,t],x,x],n[x,0]==ff[x]},{n}, {x, 0, 1}, {t, 0,30}]

matrix={{1,0,0,0,0,0,0,0,0,0,0},{p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0,0,0,0},{0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0,0,0},{0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0,0},{0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0,0},{0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0,0},{0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0,0},{0,0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0,0},{0,0,0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p),0},{0,0,0,0,0,0,0,0,p (1-p),(p^2)+(((1-p))^2),p (1-p)},{0,0,0,0,0,0,0,0,0,0,1}}/.p->(x/10);

matrix=Transpose[{{1,0,0,0,0,0,0,0,0,0,0},{(1/10) (1-(1/10)) 1,(((1-(1/10)))^2)+((1^2)/100),(1/10) (1-(1/10)) 1,0,0,0,0,0,0,0,0},{0,(1/10) (1-(2/10)) 2,(((1-(2/10)))^2)+((2^2)/100),(1/10) (1-(2/10)) 2,0,0,0,0,0,0,0},{0,0,(1/10) (1-(3/10)) 3,(((1-(3/10)))^2)+((3^2)/100),(1/10) (1-(3/10)) 3,0,0,0,0,0,0},{0,0,0,(1/10) (1-(4/10)) 4,(((1-(4/10)))^2)+((4^2)/100),(1/10) (1-(4/10)) 4,0,0,0,0,0},{0,0,0,0,(1/10) (1-(5/10)) 5,(((1-(5/10)))^2)+((5^2)/100),(1/10) (1-(5/10)) 5,0,0,0,0},{0,0,0,0,0,(1/10) (1-(6/10)) 6,(((1-(6/10)))^2)+((6^2)/100),(1/10) (1-(6/10)) 6,0,0,0},{0,0,0,0,0,0,(1/10) (1-(7/10)) 7,(((1-(7/10)))^2)+((7^2)/100),(1/10) (1-(7/10)) 7,0,0},{0,0,0,0,0,0,0,(1/10) (1-(8/10)) 8,(((1-(8/10)))^2)+((8^2)/100),(1/10) (1-(8/10)) 8,0},{0,0,0,0,0,0,0,0,(1/10) (1-(9/10)) 9,(((1-(9/10)))^2)+((9^2)/100),(1/10) (1-(9/10)) 9},{0,0,0,0,0,0,0,0,0,0,1}}];

z=0;
Show[BarChart[popsize*{0,0,1,0,0,0,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=1;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,1,0,0,0,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=3;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,1,0,0,0,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=10;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,1,0,0,0,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

z=30;
Show[BarChart[popsize*MatrixPower[matrix,z].{0,0,1,0,0,0,0,0,0,0,0},ChartLabels->Table[i,{i,0,1,0.1}],Ticks->{False,Automatic},AxesOrigin->{1,0}],
Plot[{Evaluate[n[(x-1)/popsize,z]/.sol]},{x,1,popsize+1}]
]

Figure 15.3 - Brownian Motion

norm=PDF[NormalDistribution[x0+(t-t0)μ, Sqrt[(t-t0)σ^2 ]],x]

Simplify[-μ D[norm,x]+1/2 σ^2 D[norm,x,x]-D[norm,t]]

Simplify[μ D[norm,x0]+1/2 σ^2 D[norm,x0,x0]+D[norm,t0]]

Show[
Plot[norm/.σ->1/.μ->1/.t0->0/.x0->0/.t->1,{x,-10,20},PlotRange->All],
Plot[norm/.σ->1/.μ->1/.t0->0/.x0->0/.t->2,{x,-10,20},PlotRange->All],
Plot[norm/.σ->1/.μ->1/.t0->0/.x0->0/.t->3,{x,-10,20},PlotRange->All],
Plot[norm/.σ->1/.μ->1/.t0->0/.x0->0/.t->4,{x,-10,20},PlotRange->All],
Plot[norm/.σ->1/.μ->1/.t0->0/.x0->0/.t->5,{x,-10,20},PlotRange->All],
Plot[norm/.σ->1/.μ->1/.t0->0/.x0->0/.t->10,{x,-10,20},PlotRange->All],
PlotRange->{{-5,20},{0,0.41}}
]

Figure 15.4 - Allele frequency distribution in the Wright-Fisher model

We start by deriving the drift and diffusion parameters following Box 16.3.

The discrete change in allele frequency is:

prec=((1-ν1)(1+s) p+ν2(1-p))/((1+s) p+(1-p))

The expected change in number of A alleles:

Simplify[Sum[(j-i)*PDF[BinomialDistribution[N, pprime ],j],{j,0,N}]*N^(alpha)/N^(beta)/.i->x N^beta]/.(((-1+pprime))^(-N))->(((-1))^(-N))(((1-pprime))^(-N))

Factor[%/.i->x N^beta/.alpha->1/.beta->1/.x->p]

Factor[%/.pprime->prec/.s->ψ/N/.ν1->θ1/N/.ν2->θ2/N]

Simplify[Limit[%,N->Infinity]]

The expected squared change in number of A alleles:

Simplify[Simplify[Sum[(j-i)^2*PDF[BinomialDistribution[N, pprime],j],{j,0,N}]*N^(alpha)/N^(2 beta)/.i->x N^beta]/.(((-1+pprime))^(-N))->(((-1))^(-N))(((1-pprime))^(-N))/.(((-1+pprime))^N)->(((-1))^N)(((1-pprime))^N)/.(((-1+pprime))^(1+N))->(((-1))^(1+N))(((1-pprime))^(1+N))/.((((1/(1-pprime))))^N) ->(((1-pprime))^(-N)) ,{N>0}]

Factor[%/.i->x N^beta/.alpha->1/.beta->1/.x->p]

Simplify[%/.pprime->prec/.s->ψ/N/.ν1->θ1/N/.ν2->θ2/N]

Simplify[Limit[%,N->Infinity]]

The expected cubed change in number of A alleles:

Simplify[Simplify[Sum[(j-i)^3*PDF[BinomialDistribution[N, pprime],j],{j,0,N}]*N^(alpha)/N^(3 beta)/.i->x N^beta]/.(((-1+pprime))^(-N))->(((-1))^(-N))(((1-pprime))^(-N))/.(((-1+pprime))^N)->(((-1))^N)(((1-pprime))^N)/.(((-1+pprime))^(1+N))->(((-1))^(1+N))(((1-pprime))^(1+N))/.((((1/(1-pprime))))^N) ->(((1-pprime))^(-N)) ,{N>0}]

Factor[%/.i->x N^beta/.alpha->1/.beta->1/.x->p]/.pprime->rec/.s->ψ/N/.ν1->θ1/N/.ν2->θ2/N;

Simplify[Limit[%,N->Infinity]]

The Wright-Fisher diffusion process has a stationary distribution of:

Ax[x_]=Simplify[Integrate[2 ((1-x)θ2-x θ1+(1-x) x ψ)/((1-x) x),x]]

Simplify[Exp[Ax[x]]/((1-x) x)]

In the presence of selection, the solution involves numerical integration (represented by the hypergeometric functions):

Clear[stationary]
stationary[ψ_,θ1_,θ2_,x_]=Simplify[(E^(2 x ψ)) (((1-x))^(-1+2 θ1)) (x^(-1+2 θ2))/Integrate[(E^(2 x ψ)) (((1-x))^(-1+2 θ1)) (x^(-1+2 θ2)),{x,0,1}],{Re[θ1]>0&&Re[θ2]>0}]

In the absence of selection, this integrates to:

Simplify[(((1-x))^(-1+2 θ1)) (x^(-1+2 θ2))/Integrate[(((1-x))^(-1+2 θ1)) (x^(-1+2 θ2)),{x,0,1}],{Re[θ1]>0&&Re[θ2]>0}]

Using rules relating Gamma and sine functions, this can be more simply written as:

stationary[0,θ1_,θ2_,x_]=(((((1-x))^(-1+2 θ1)) (x^(-1+2 θ2)) Gamma[2 (θ1+θ2)])/(Gamma[2 θ1] Gamma[2 θ2]));

FullSimplify[%/%%]

In the absence of selection (equal forward and backward mutation rates):

Show[
Plot[stationary[0,0.1,0.1,x],{x,0,1},PlotRange->All],
Plot[stationary[0,0.5,0.5,x],{x,0,1},PlotRange->All],
Plot[stationary[0,1,1,x],{x,0,1},PlotRange->All],
Plot[stationary[0,5,5,x],{x,0,1},PlotRange->All],
PlotRange->{Automatic,{0,4.1}},AxesOrigin->{0,0}
]

With selection (equal forward and backward mutation rates):

Show[
Plot[stationary[-10,1,1,x],{x,0,1},PlotRange->All],
Plot[stationary[-2,1,1,x],{x,0,1},PlotRange->All],
Plot[stationary[0,1,1,x],{x,0,1},PlotRange->All],
Plot[stationary[2,1,1,x],{x,0,1},PlotRange->All],
Plot[stationary[10,1,1,x],{x,0,1},PlotRange->All]
]

Figure 15.5 - Probability of Fixation in the Wright-Fisher model

The Wright-Fisher diffusion process in the absence of mutation gives:

Ax[x_]=Simplify[Integrate[2 ((1-x) x ψ)/((1-x) x),x]]

Sx[x_]=Simplify[Integrate[Exp[-Ax[x]],x]]

The fixation probability allows us to answer the question:  what is the probability that the process will reach carrying capacity before going extinct:

Simplify[(Sx[x0]-Sx[0])/(Sx[N/N^beta]-Sx[0])]

Clear[fix,num,s]
fix[p0_]=Simplify[(Sx[x0]-Sx[0])/(Sx[N/N^beta]-Sx[0])/.x0->p0/.ψ->s N/.beta->1]

Limit[fix[p0],s->0]

The exact fixation probability (given by the left eigenvector associated with one):

num=10;

fix10=Table[{s,mat=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,0,num},{i,0,num}];
eigenvec=Table[e[i],{i,0,num}];
e[1]/.(Solve[{eigenvec==eigenvec.mat,e[0]==0,e[num]==1},eigenvec]//Flatten)},{s,-0.1,0.1,0.01}]

num=100;

(Greater precision was needed in the values of s used for these tiny fixation probabilties:)

fix100=Table[{N[s,5],mat=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,0,num},{i,0,num}];
eigenvec=Table[e[i],{i,0,num}];
e[1]/.(Solve[{eigenvec==eigenvec.mat,e[0]==0,e[num]==1},eigenvec]//Flatten)},{s,N[-1/10,32],N[1/10,32],N[1/100,32]}]

Show[
ListPlot[fix10,PlotStyle->PointSize[0.01]],
ListPlot[fix100,PlotStyle->PointSize[0.01]],
Plot[fix[p0]/.p0->1/N/.N->10,{s,-0.1,0.25}],
Plot[fix[p0]/.p0->1/N/.N->100,{s,-0.1,0.25}],
Plot[fix[p0]/.p0->1/N/.N->1000,{s,-0.1,0.25}],
PlotRange->{{-0.1,0.2},{-0.02,0.2}},
AxesOrigin->{0,0}]

Export["Fig15-5.eps",%];

Show[ListLogPlot[fix100,PlotStyle->PointSize[0.01]],LogPlot[fix[p0]/.p0->(1/N)/.N->100,{s,-0.1,0.25}]]

Figure 15.6 - Fixation time in the Wright-Fisher model (no selection)

For the neutral Wright-Fisher model, the waiting time until fixation is given by the diffusion approximation:

Simplify[2(Integrate[1,{y,p0,1}]+(1-p0)/p0 Integrate[y/(1-y),{y,0,p0}]),p0<1]

Normal[Series[%,{p0,0,0}]]

These fixation times are in units of N generations.

The mean time until fixation:

Plot[((2 (-1+p0) Log[1-p0])/p0),{p0,0,1}]

num=100;
s=0.;

The transition matrix excluding boundaries:

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

The fixation probabilities:

tabF=Table[(j/N)/(i/N)/.N->num,{j,1,num-1},{i,1,num-1}];

times100F=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*tab]

Show[
ListPlot[Table[{i/num,times100F[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.006]],Plot[((2 (-1+p0) Log[1-p0])/p0),{p0,0,1}],
PlotRange->All]

num=200;
s=0.;

The transition matrix excluding boundaries:

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

The fixation probabilities:

tabF=Table[(j/N)/(i/N)/.N->num,{j,1,num-1},{i,1,num-1}];

times200F=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*tab];

Show[
ListPlot[Table[{i/num,times200F[[i]]/num},{i,1,num-1}]],Plot[((2 (-1+p0) Log[1-p0])/p0),{p0,0,1}],
PlotRange->All]

The mean time until loss or fixation:

Plot[p0 ((2 (-1+p0) Log[1-p0])/p0)+ (1-p0)((2 (-p0) Log[p0])/(1-p0)),{p0,0,1}]

Figure 15.7 - Fixation time in the Wright-Fisher model (with selection)

The Wright-Fisher diffusion process in the absence of mutation gives:

Ax[x_]=Simplify[Integrate[2 ((1-x) x ψ)/((1-x) x),x]]

Sx[x_]=Simplify[Integrate[Exp[-Ax[x]],x]]

Clear[fix,num,s]
fix[p0_]=Simplify[(Sx[x0]-Sx[0])/(Sx[N/N^beta]-Sx[0])/.x0->p0/.beta->1]

To calculate the time to fixation:

Simplify[fix[y] (1-fix[y])/(Exp[-Ax[y]]  y (1-y))/.beta->1]

Integrate[%,y];
part1=(Limit[%,y->1])-(%/.y->p0);

Simplify[fix[y]^2/(Exp[-Ax[y]]  y (1-y))/.beta->1]

Integrate[%,y];
part2=(%/.y->p0)-(Limit[%,y->0]);

fixtime[p0_,ψ_]=2 (Sx[1]-Sx[0])(part1+(1-fix[p0])/fix[p0]part2)/.beta->1;

In the neutral case:

Simplify[Limit[fix[ y] (1-fix[y])/(Exp[-Ax[y]]  y(1-y))/.beta->1,ψ->0]]

Integrate[%,y];
part1=(%/.y->1)-(%/.y->p0);

Simplify[Limit[fix[y]^2/(Exp[-Ax[y]]y(1-y))/.beta->1,ψ->0]]

Integrate[%,y];
part2=(%/.y->p0)-(Limit[%,y->0]);

neutralfixtime=Simplify[Limit[2 (Sx[1]-Sx[0])(part1+(Sx[1]-Sx[p0])/(Sx[p0]-Sx[0])part2),ψ->0]/.beta->1]

The exact fixation time is given by:

num=100;
s=0.1/num;

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

fixE=(Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],num],{i,1,num-1}].Inverse[IdentityMatrix[num-1]-N[Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}],30]]);

tabF=Table[fixE[[j]]/fixE[[i]]/.ψ->N s/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fa=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*tab];

num=100;
s=1.0/num;

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

fixE=(Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],num],{i,1,num-1}].Inverse[IdentityMatrix[num-1]-N[Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}],30]]);

tabF=Table[fixE[[j]]/fixE[[i]]/.ψ->N s/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fb=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*tab];

num=100;
s=10/num;

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

fixE=(Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],num],{i,1,num-1}].Inverse[IdentityMatrix[num-1]-N[Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}],30]]);

tabF=Table[fixE[[j]]/fixE[[i]]/.ψ->N s/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fc=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*N[tab,30]];

Show[
ListPlot[Table[{i/num,times100Fa[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Fb[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Fc[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],Plot[Re[fixtime[p0,0.1]],{p0,0,1}],
Plot[Re[fixtime[p0,1]],{p0,0,1}],
Plot[Re[fixtime[p0,10]],{p0,0,1}],
PlotRange->{Automatic,{0,2}}]

For negative selection:

num=100;
s=1/(1+0.1/num)-1;

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

fixE=(Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],num],{i,1,num-1}].Inverse[IdentityMatrix[num-1]-N[Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}],30]]);

tabF=Table[fixE[[j]]/fixE[[i]]/.ψ->N s/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fd=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*tab];

num=100;
s=1/(1+1.0/num)-1;

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

tabF=Table[fix[j/N]/fix[i/N]/.ψ->N s/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fe=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tabF*tab];

num=100;
s=s=1/(1+10./num)-1;

tab=Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}]/.Indeterminate->1;

fixE=(Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],num],{i,1,num-1}].Inverse[IdentityMatrix[num-1]-N[Table[PDF[BinomialDistribution[num,i/num*(1+s)/(1+s i/num )],j],{j,1,num-1},{i,1,num-1}],30]]);

tabF=Table[fixE[[j]]/fixE[[i]]/.ψ->N s/.N->num,{j,1,num-1},{i,1,num-1}];

times100Ff=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-N[tabF,50]*N[tab,50]];

Show[
ListPlot[Table[{i/num,times100Fd[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Fe[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Ff[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],Plot[Re[fixtime[p0,-0.1]],{p0,0,1}],
Plot[Re[fixtime[p0,-1]],{p0,0,1}],
Plot[Re[fixtime[p0,-10]],{p0,0,1}],
PlotRange->{Automatic,{0,2}}]

Altogether:

Show[
ListPlot[Table[{i/num,times100Fa[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Fb[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Fc[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],Plot[Re[fixtime[p0,0.1]],{p0,0,1}],
Plot[Re[fixtime[p0,1]],{p0,0,1}],
Plot[Re[fixtime[p0,10]],{p0,0,1}],
ListPlot[Table[{i/num,times100Fd[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Fe[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],ListPlot[Table[{i/num,times100Ff[[i]]/num},{i,1,num-1}],PlotStyle->PointSize[0.007]],Plot[Re[fixtime[p0,-0.1]],{p0,0,1}],
Plot[Re[fixtime[p0,-1]],{p0,0,1}],
Plot[Re[fixtime[p0,-10]],{p0,0,1}],
PlotRange->{Automatic,{0,2}}]

Figure 15.8 - Establishment and extinction times in the birth-death model

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 no longer concerned about the viability of a population.

First we show that the process can be modeled as a diffusion:

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

E[j-i] is:

exp1=Simplify[birth(+1) +death(-1)]

E[(j-i)^2] is:

exp2=Simplify[birth(+1)^2 +death(-1)^2]

E[(j-i)^3] is:

exp3=Simplify[birth(+1)^3 +death(-1)^3]

Scaled mean change:

mom1=Simplify[(exp1)*N^(alpha)/N^(beta)/.i->x N^beta]

Scaled squared change:

mom2=Simplify[(exp2)*N^(alpha)/N^(2 beta)/.i->x N^beta]

Scaled cubed change:

mom3=Simplify[(exp3)*N^(alpha)/N^(3 beta) /.i->x N^beta]

Simplify[{mom1,mom2,mom3}]

We choose alpha = beta on the basis of the variance:

Simplify[{mom1,mom2,mom3}/.alpha->beta]

If (b-d) is O(1), then there is no choice of beta that forces the third moment to be zero and the second moment to be finite.  Thus, we define rho as (b-d) N^beta.  Let's also define b+d = lambda as the per capita rate of events.

subbd=Solve[{b+d==lambda,(b-d)N^beta==rho},{b,d}]//Flatten

subbdR=Solve[{b+d==lambda,(b-d)N^beta==rho},{rho,lambda}]//Flatten

Simplify[{mom1,mom2,mom3}/.subbd/.alpha->beta]

A natural choice of beta = 1, and the moments are {rho x,lambda x,0}.

Simplify[{mom1,mom2,mom3}/.subbd/.alpha-> beta/.beta->1]

Limit[%,N->Infinity]

In this case, we have the birth-death diffusion process.

Ax[x_]=Simplify[Integrate[2 (rho x)/(lambda x),x]]

Sx[x_]=Simplify[Integrate[Exp[-Ax[x]],x]]

The fixation probability allows us to answer the question:  what is the probability that the process will reach carrying capacity before going extinct:

Simplify[(Sx[x0]-Sx[0])/(Sx[N/N^beta]-Sx[0])]

Simplify[(Sx[x0]-Sx[0])/(Sx[N/N^beta]-Sx[0])/.subbdR/.x0->n0/N^beta/.beta->1]

Which can be rewritten as:

fix[n0_]= (( -(E^((-2 (1-d/b) n0)/(1+d/b)))+1)/(-(E^((- 2 (1-d/b) N)/(1+d/b)))+1));

Limit[fix[n0],b->d]

For a branching process, the exact fixation probability is also known:

u[n0_]=Simplify[1 - Sum[(d/b)^i,{i,n0,N-1}]/Sum[(d/b)^i,{i,0,N-1}]]

The establishment probability is well approximated by the diffusion:

Show[
Plot[fix[n0]/.d->0.99 b/.N->1000,{n0,0,1000},PlotRange->All],
ListPlot[Table[{n0,u[n0]/.d->0.99 b/.N->1000},{n0,0,1000,10}]]
]

To calculate the waiting time until establishment (hits state N):

Simplify[fix[ y N] (1-fix[y N])/(Exp[-Ax[y]] lambda y)/.subbd/.beta->1]

Integrate[%,y];
part1=(Limit[%,y->1])-(%/.y->x0);

Simplify[fix[y N]^2/(Exp[-Ax[y]] lambda y)/.subbd/.beta->1]

Integrate[%,y];
part2=(%/.y->x0)-(Limit[%,y->0]);

fixtime=2 (Sx[1]-Sx[0])(part1+(1-fix[n0])/fix[n0]part2)/.subbdR/.x0->n0/N^beta/.beta->1;

To calculate the waiting time until extinction (hits state 0):

Simplify[(1-fix[y N])^2/(Exp[-Ax[y]] lambda y)/.subbd/.beta->1]

Integrate[%,y];
part1=(Limit[%,y->1])-(%/.y->x0);

Simplify[fix[ y N] (1-fix[ y N] )/(Exp[-Ax[y]] lambda y)/.subbd/.beta->1]

Integrate[%,y];
part2=(%/.y->x0)-(Limit[%,y->0]);

losstime=2 (Sx[1]-Sx[0])(fix[n0]/(1-fix[n0])part1+part2)/.subbdR/.x0->n0/N^beta/.beta->1;

In the neutral case:

Simplify[Limit[fix[ y N] (1-fix[y N])/(Exp[-Ax[y]] lambda y)/.subbd/.beta->1,rho->0]]

Integrate[%,y];
part1=(%/.y->1)-(%/.y->x0);

Simplify[Limit[fix[y N]^2/(Exp[-Ax[y]] lambda y)/.subbd/.beta->1,rho->0]]

Integrate[%,y];
part2=(%/.y->x0)-(Limit[%,y->0]);

neutralfixtime=Simplify[Limit[2 (Sx[N/N^beta]-Sx[0])(part1+(Sx[N/N^beta]-Sx[x0])/(Sx[x0]-Sx[0])part2),rho->0]/.subbdR/.x0->n0/N^beta/.beta->1]

The time to loss in the neutral case is:

Simplify[Limit[(1-fix[y N])^2 /(Exp[-Ax[y]] lambda y)/.subbd/.beta->1,rho->0]]

Integrate[%,y];
part1=(%/.y->1)-(%/.y->x0);

Simplify[Limit[(1-fix[y N])fix[y N] /(Exp[-Ax[y]] lambda y)/.subbd/.beta->1,rho->0]]

Integrate[%,y];
part2=(%/.y->x0)-(Limit[%,y->0]);

neutrallosstime=Simplify[Limit[2 (Sx[N/N^beta]-Sx[0])((Sx[x0]-Sx[0])/(Sx[N/N^beta]-Sx[x0])part1+part2),rho->0]/.subbdR/.x0->n0/N^beta/.beta->1]

These are consistent with the neutral time until absorption:  -((2 n0 Log[(n0/N)])/((b+d) N)):

Simplify[neutrallosstime (1-n0/N)+neutralfixtime(n0/N)]

The exact fixation time for a branching process is given in Ewen's (p. 71-72).

num=100;
bb=0.002;
dd=0.002;

tab=Table[If[j==i,1-bb i - dd i,If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num-1},{i,1,num-1}];

tabF=Table[(j/N)/(i/N)/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

times100F=Table[1,{i,1,num-1}].(tabF*Inverse[IdentityMatrix[num-1]-tab]);

tabL=Table[(1-j/N)/(1-i/N)/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

times100L=Table[1,{i,1,num-1}].(tabL*Inverse[IdentityMatrix[num-1]-tab]);

plot0F=Show[
ListPlot[Table[{i,times100F[[i]]/num},{i,1,num-1,2}]],Plot[Re[neutralfixtime/.N->100/.b->bb/.d->dd],{n0,0,num-1}],
PlotRange->All]

plot0L=Show[
ListPlot[Table[{i,times100L[[i]]/num},{i,1,num-1,2}]],Plot[Re[neutrallosstime/.N->100/.b->bb/.d->dd],{n0,0,num-1}],
PlotRange->All]

num=100;
bb=0.0021;
dd=0.0019;

tab=Table[If[j==i,1-bb i - dd i,If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num-1},{i,1,num-1}];

tabF=Table[u[j]/u[i]/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

tabL=Table[(1-u[j])/(1-u[i])/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fa=Table[1,{i,1,num-1}].(tabF*Inverse[IdentityMatrix[num-1]-tab]);

times100La=Table[1,{i,1,num-1}].(tabL*Inverse[IdentityMatrix[num-1]-tab]);

times100A=Table[1,{i,1,num-1}].Inverse[IdentityMatrix[num-1]-tab];

The average of the loss and fixation times is correctly the absorption time:

Table[u[i] times100Fa[[i]]+(1-u[i])times100La[[i]]/.b->bb/.d->dd/.N->num,{i,1,num-1}]/times100A

Note that the time to fixation conditional on fixation is the same if we switch the magnitudes of b and d.  The fixation probability is now lower because the birth rate is lower, but conditional on fixation the time is the same.

num=100;
dd=0.0021;
bb=0.0019;

tab=Table[If[j==i,1-bb i - dd i,If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num-1},{i,1,num-1}];

tabF=Table[u[j]/u[i]/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

tabL=Table[(1-u[j])/(1-u[i])/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fb=Table[1,{i,1,num-1}].(tabF*Inverse[IdentityMatrix[num-1]-tab]);

times100Lb=Table[1,{i,1,num-1}].(tabL*Inverse[IdentityMatrix[num-1]-tab]);

plot1F=Show[
ListPlot[Table[{i,times100Fa[[i]]/num},{i,1,num-1,2}]],ListPlot[Table[{i,times100Fb[[i]]/num},{i,1,num-1,2}]],Plot[Re[fixtime/.N->100/.b->bb/.d->dd],{n0,0,num-1}],
Plot[Re[fixtime/.N->100/.b->dd/.d->bb],{n0,0,num-1}],
PlotRange->All]

plot1L=Show[
ListPlot[Table[{i,times100La[[i]]/num},{i,1,num-1,2}]],ListPlot[Table[{i,times100Lb[[i]]/num},{i,1,num-1,2}]],Plot[Re[losstime/.N->100/.b->bb/.d->dd],{n0,0,num-1}],
Plot[Re[losstime/.N->100/.b->dd/.d->bb],{n0,0,num-1}],
PlotRange->All]

Increasing the difference:

num=100;
bb=0.0022;
dd=0.0018;

tab=Table[If[j==i,1-bb i - dd i,If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num-1},{i,1,num-1}];

tabF=Table[u[j]/u[i]/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

tabL=Table[(1-u[j])/(1-u[i])/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fc=Table[1,{i,1,num-1}].(tabF*Inverse[IdentityMatrix[num-1]-tab]);

times100Lc=Table[1,{i,1,num-1}].(tabL*Inverse[IdentityMatrix[num-1]-tab]);

num=100;
dd=0.0022;
bb=0.0018;

tab=Table[If[j==i,1-bb i - dd i,If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num-1},{i,1,num-1}];

tabF=Table[u[j]/u[i]/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

tabL=Table[(1-u[j])/(1-u[i])/.b->bb/.d->dd/.N->num,{j,1,num-1},{i,1,num-1}];

times100Fd=Table[1,{i,1,num-1}].(tabF*Inverse[IdentityMatrix[num-1]-tab]);

times100Ld=Table[1,{i,1,num-1}].(tabL*Inverse[IdentityMatrix[num-1]-tab]);

plot2F=Show[
ListPlot[Table[{i,times100Fc[[i]]/num},{i,1,num-1,2}]],ListPlot[Table[{i,times100Fd[[i]]/num},{i,1,num-1,2}]],Plot[Re[fixtime/.N->100/.b->bb/.d->dd],{n0,0,num-1}],
Plot[Re[fixtime/.N->100/.b->dd/.d->bb],{n0,0,num-1}],
PlotRange->All]

plot2L=Show[
ListPlot[Table[{i,times100Lc[[i]]/num},{i,1,num-1,2}]],ListPlot[Table[{i,times100Ld[[i]]/num},{i,1,num-1,2}]],Plot[Re[losstime/.N->100/.b->bb/.d->dd],{n0,0,num-1}],
Plot[Re[losstime/.N->100/.b->dd/.d->bb],{n0,0,num-1}],
PlotRange->All]

Show[plot0F,plot1F,plot2F]

Show[plot0L,plot1L,plot2L]

Note that the diffusion results apply to any population size assuming that (b+d) and (b-d)N are constant.

Solve[{0.004==(b+d),0.0002(100)== (b-d) N},{b,d}]//Flatten

Solve[{0.004==(b+d),0.0004(100)== (b-d) N},{b,d}]//Flatten

Show[Plot[Re[fixtime/.N->100/.b->0.0021/.d->0.0019/.n0->x*100],{x,0,1}],
Plot[Re[fixtime/.b->0.002+(0.01/N)/.d->0.002-(0.01/N)/.N->1000/.n0->x*1000],{x,0,1},PlotStyle->{Thick,Dashed}],
PlotRange->All
]

Figure 15.9 - Fixation time in the birth-death model with a hard carrying capacity of K

Here we explore the birth-death model where birth=b i and death=d i until a hard carrying capacity (K) is reached at which point the probability of birth is zero does meet the requirements for a diffusion approximation:

birth=b i;
death=d i;

E[j-i] is:

exp1=Simplify[birth(+1) +death(-1)]

E[(j-i)^2] is:

exp2=Simplify[birth(+1)^2 +death(-1)^2]

E[(j-i)^3] is:

exp3=Simplify[birth(+1)^3 +death(-1)^3]

Scaled mean change:

mom1=Simplify[(exp1)*K^(alpha)/K^(beta)/.p->i/K/.i->x K^beta]

Scaled squared change:

mom2=Simplify[(exp2)*K^(alpha)/K^(2 beta)/.p->i/K /.i->x K^beta]

Scaled cubed change:

mom3=Simplify[(exp3)*K^(alpha)/K^(3 beta)/.p->i/m /.i->x K^beta]

Expand[{mom1,mom2,mom3}/.subbd/.alpha->1/.beta->1/.gamma->1/.test->0]

Limit[%,K->Infinity]

Ax[x_]=Simplify[Integrate[2 ((rho x)/(lambda x)),x]/.rho->b-d/.lambda->b+d,{b>0,d>0}]

sigma[x_]=x (b+d)

Simplify[ Exp[Ax[y]]/sigma[y]]

Simplify[(E^(-Ax[z]))]

The diffusion approximation for the case of exponential growth (no density dependence):

Clear[timeE]
timeE[b_,d_,K_,n0_]:=timeE[b,d,K,n0]=NIntegrate[2 (E^(-((2 (b-d) z)/(b+d)))) ((E^((2 (b-d) y)/(b+d)))/((b+d) y)),{z,0,n0},{y,z,K}]

num=100;
bb=0.002;
dd=0.002;

tabE=Table[If[j==i,If[j==num,1 - dd i,1-bb i - dd i],If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num},{i,1,num}];

times100E=Table[1,{i,1,num}].(Inverse[IdentityMatrix[num]-N[tabE,30]]);

plotA=Show[ListPlot[Table[{i,((times100E[[i]])/num)},{i,1,num-1,2}]],ListPlot[Join[Table[{n0,((timeE[bb,dd,num,n0])/num)},{n0,0,5}],Table[{n0,((timeE[bb,dd,num,n0])/num)},{n0,10,num,5}]],Joined->True],PlotRange->{Automatic,{0,500}}]

Greater precision is needed when births exceed deaths:

num=100;
bb=22/10^4;
dd=18/10^4;

tabE=Table[If[j==i,If[j==num,1 - dd i,1-bb i - dd i],If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num},{i,1,num}];

times100E=Table[1,{i,1,num}].(Inverse[IdentityMatrix[num]-N[tabE,30]]);

times100E[[1]]/num

times100E[[100]]/num

plotB=Show[ListPlot[Table[{i,((times100E[[i]])/num)},{i,1,num-1,2}]],ListPlot[Join[Table[{n0,((timeE[bb,dd,num,n0])/num)},{n0,0,5}],Table[{n0,((timeE[bb,dd,num,n0])/num)},{n0,10,num,5}]],Joined->True],PlotRange->{Automatic,{0,1000000000}},AxesOrigin->{0,0}]

num=100;
bb=0.0018;
dd=0.0022;

tabE=Table[If[j==i,If[j==num,1 - dd i,1-bb i - dd i],If[j==i+1,bb i,If[j==i-1,dd i,0]]],{j,1,num},{i,1,num}];

times100E=Table[1,{i,1,num}].(Inverse[IdentityMatrix[num]-N[tabE,30]]);

plotC=Show[ListPlot[Table[{i,((times100E[[i]])/num)},{i,1,num-1,2}]],ListPlot[Join[Table[{n0,((timeE[bb,dd,num,n0])/num)},{n0,0,5}],Table[{n0,((timeE[bb,dd,num,n0])/num)},{n0,10,num,5}]],Joined->True],PlotRange->{Automatic,{0,500}}]

Show[plotA,plotB,plotC,
PlotRange->{Automatic,{0,600}}]

Figure 15.10 - Sperm migration

Show[
Plot[((-Exp[((-2 μ x0)/σ)]+1)/(-Exp[((-2 μ b)/σ)]+1))/.{b->1,x0->.01,σ->1},{μ,0,10}],
Plot[((-Exp[((-2 μ x0)/σ)]+1)/(-Exp[((-2 μ b)/σ)]+1))/.{b->1,x0->.01,σ->3},{μ,0,10}],
PlotRange->{{0,10},{0,0.06}}
]

Export["fig15-10.eps",%];

Figure 15SM3.2 - Bezykin model

Clear[s,m]

Letting m = σ2, Bezykin's model is:

p''[x]==(2 s/m)(1-p) p (1-2 p)

(Actually, this differs from Bezykin by a factor of 2, because he defined migration slightly differently than is standard in the Brownian motion model.)

Solving with Mathematica (takes a while and gives an uninterpretable result):

DSolve[p''[x]==(2 s/m)(1-p[x]) p[x] (1-2 p[x]),p[x],x]

Using the condition that p[0] = 1/2 (fails):

DSolve[{p''[x]==(2 s/m)(1-p[x]) p[x] (1-2 p[x]),p[0]==1/2},p[x],x]

Solving by guessing the form that dp/dx = c p[x] (1-p[x]):

Factor[D[c(1-p[x]) p[x],x]/.p'[x]->c p[x] (1-p[x])]

Solving by substitution:  First define u[x] = p[x] (1-p[x]) and note that du/dx = (1-2p[x]) dp/dx.  Consequently, (d^2)p/d(x^2)becomes (2 s/m) u[x] du/dx /(dp/dx).  Then let A = dp/dx to get that (d^2)p/d(x^2)= dA/dx = (2 s/m) u[x] du/dx /A.  Rearranging, A dA = (2 s/m) u du.  Integrating both sides and noting that dp/dx must equal zero when u[x] =0, we have the solution (A^2)= (2 s/m) (u^2)so that dp/dx = (Sqrt[2 s/m]) p[x] (1-p[x]):

Show[
Plot[Exp[Sqrt[2 s/m] x]/(1+Exp[Sqrt[2 s/m] x])/.s->0.5/.m->0.1,{x,-15,15}],
Plot[Exp[Sqrt[2 s/m] x]/(1+Exp[Sqrt[2 s/m] x])/.s->0.05/.m->0.1,{x,-15,15}],
Plot[Exp[Sqrt[2 s/m] x]/(1+Exp[Sqrt[2 s/m] x])/.s->0.005/.m->0.1,{x,-15,15}],
AspectRatio->0.25
]

The form of the solution makes it clear that we should be able to get the same results by increasing m as we do by decreasing s:

Show[
Plot[Exp[Sqrt[2 s/m] x]/(1+Exp[Sqrt[2 s/m] x])/.s->0.05/.m->0.01,{x,-10,10}],
Plot[Exp[Sqrt[2 s/m] x]/(1+Exp[Sqrt[2 s/m] x])/.s->0.05/.m->0.1,{x,-10,10}],
Plot[Exp[Sqrt[2 s/m] x]/(1+Exp[Sqrt[2 s/m] x])/.s->0.05/.m->1,{x,-10,10}],
AspectRatio->0.25
]