restart:

F:=x->exp(-int(2*a(x)/b(x),x)); this is equation (2) in Box 3.1u:=p->int(F(x),x=0..p)/int(F(x),x=0..1); this is the probability of fixation (equation 1 in Box 3.1)t[low]:=2/(b(x)*F(x))*int(F(y),y=p..1); this is equation (2) in Box 3.2t[high]:=2/(b(x)*F(x))*int(F(y),y=x..1); this is equation (2) in Box 3.T:=int(t[low],x=0..p)+int(t[high],x=p..1); this is equation (1) in Box 3.2Phi(x):=2*int(F(x),x=0..1)/(b(x)*F(x));tau:=int(Phi(x)*u(x)*(1-u(x)),x=p..1)+(1-u(p))/u(p)*int(Phi(x)*u(x)^2,x=0..p);

neut:={a(x)=0,b(x)=x*(1-x)/(2*N),a(y)=0,b(y)=y*(1-y)/(2*N)}; mean and variance for a neutral alleleU_n:=simplify(subs(neut,u(p))); fixation probability of a neutral alleletau_n:=simplify(subs(neut,tau));(conditional) time until fixation of a neutral alleletau[n]:=-4*N*(1-p)/p*ln(1-p); Kimura and Ohta, 1969, Eq.14limit(tau[n],p=0); this is the waiting time until a neutral allele introduced in a single copy gets fixedsubs(p=1-1/(2*N),tau[n]); this is the waiting time until a neutral allele introduced in a single copy gets lost

ad:={a(x)=s*x*(1-x),b(x)=x*(1-x)/(2*N),a(y)=s*y*(1-y),b(y)=y*(1-y)/(2*N)}; mean and variance for an additive alleleU_a:=simplify(subs(ad,u(p))); fixation probability of an additive alleletau_a:=simplify(subs(ad,tau));(conditional) time until fixation of a neutral allelefactor(%);

T_n:=simplify(subs(neut,T)); average time until a neutral allele is lost or fixed NOT APPROPRIATE FORMULA: TWO ABSORBING STATES!!limit(x^x,x=0);

restart:plotsetup(window):plotsetup(default):plotsetup(cps,plotoutput=`./out.ps`,plotoptions=`portrait,noborder,height=5in,width=5in`);

f:=proc() if rand()/10^(12)>p then 1 else 0 fi end;with(plots):p:=.15: size:=50:densityplot(f,0..1,0..1,axes=NONE,grid=[size,size],title="p=0.15", titlefont=[TIMES,BOLDITALIC,12]);p:=.3: size:=50:densityplot(f,0..1,0..1,axes=NONE,grid=[size,size],title="p=0.30", titlefont=[TIMES,BOLDITALIC,12]);p:=.45: size:=50:densityplot(f,0..1,0..1,axes=NONE,grid=[size,size],title="p=0.45", titlefont=[TIMES,BOLDITALIC,12]);p:=.6: size:=50:densityplot(f,0..1,0..1,axes=NONE,grid=[size,size],title="p=0.60", titlefont=[TIMES,BOLDITALIC,12]);

restart:rho:=(N-K-1)!*(N-d)!/N!/(N-K-1-d)!; this is corrected formula: K->K-1simplify(rho);g:=z->exp(-z)*z^(z-1/2)*sqrt(2(Pi));g(N-K+1)*g(N-d+1)/g(N+1)/g(N-K-d+1);simplify(%);N:=100: r0:=subs(K=0,rho): r1:=subs(K=1,rho): r2:=subs(K=2,rho):r4:=subs(K=4,rho): r8:=subs(K=8,rho): r16:=subs(K=16,rho):plot({r0,r1,r2,r4,r8,r16},d=0..30,color=black,labels=["distance","correlation"],thickness=3,font=[TIMES,BOLDITALIC,14]);?colorsimplify(subs(K=0,rho));simplify(subs(K=1,rho));

restart:f:=nu*exp(-nu*t); this is the density of the waiting time till mutation given the gametic mutation rate is nu

assume(nu>0):int(t*f,t=0..infinity); average waiting time - OKint(t^2*f,t=0..infinity)-1/nu^2; variance of the waiting time - OK (exponential distribution)int(f,t);

nu:=mu*i; gametic mutation rate is gene mutation rate mu times the number of accessible neighbors ig:=exp(-n)*n^i/i!; distribution of the number of neighbors is Poisson with mean nF:=simplify(f*g); this is a contribution of term i into the distribution function of the time till mutation given the number of neighbors is randomsimplify(sum(F,i=2..infinity));assume(mu>0): assume(i>0):X:=int(t*F,t=0..infinity); this is the contribution of term i to mean waiting timet_mean:=simplify(sum(X,i=1..infinity)); average time, only genotypes with more than 2 neighbors contributeplot({n*exp(-n)*hypergeom([1, 1],[2, 2],n)},n=1..20); coefficient in front of 1/(mu*n)Y:=int(t^2*F,t=0..infinity);t2_mean:=simplify(sum(Y,i=1..infinity));t_var:=simplify(t2_mean-t_mean^2);plot({-n*exp(-n)*(-2*hypergeom([1, 1, 1],[2, 2, 2],n)+n*exp(-n)*hypergeom([1, 1],[2, 2],n)^2)},n=1..20);R:=simplify(t_var/t_mean^2);plot({R},n=1..20);fopen("dispersion.dat",WRITE):for i from 1 to 100 do N:=20*i/100; fprintf("dispersion.dat","%f %f \n",N,evalf(subs(n=N,R))): od:fclose(fd):

restart:p:=exp(-N)*N^i/i!;f:=simplify(sum(p,i=a..infinity));simplify(sum(p,i=0..a));GAMMA(2);N:=5: plot(log10(f),a=20..30);evalf(subs(a=25,f));f;2^100.;evalf(subs(a=25,f))*%;

restart:q:=exp(-(t/tau)^b);d:=L/2*(1-q);L:=24; d5:=subs({b=.97,tau=25},d): d2:=subs({b=.91,tau=88},d): d1:=subs({b=.82,tau=406},d):d09:=subs({b=.78,tau=547},d): d08:=subs({b=.71,tau=762},d): d07:=subs({b=.63,tau=1282},d):d06:=subs({b=.51,tau=2143},d): d055:=subs({b=.39,tau=2794},d):plot({d5,d1,d09,d08,d06,d07,d055},t=0..1000);d5;

restart:with(networks):H:=cube(5): draw(H);G:=delete({1,5,10,11,15,20,25,27,29,30},H): draw(G);components(G);with(linalg):max(evalf(eigenvals(adjacency(G))));evalf(eigenvals(adjacency(G)));NiskIStjMS1LIykhIzYkIiIiIiIhJCEiIkYoRiZGKUYmRilGJkYpQ:=graph({0,1,2,3,4,5},[[0,1],[0,2],[0,3],[0,4],[0,5]]):draw(Concentric([0],[1,2,3,4,5]),Q);Q:=graph({0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15},[[0,1],[0,2],[0,3],[0,4],[0,5],[1,6],[1,7],[2,8],[2,9],[3,10],[3,11],[4,12],[4,13],[5,14],[5,15]]):draw(Concentric([0],[1,2,3,4,5]),Q);K:=adjacency(Q); Y:=K+transpose(K);max(evalf(eigenvals(Y)));rowdim(Y);coldim(Y);evalf(eigenvects(Y));plotsetup(cps,plotoutput=`./galaxy.ps`,plotoptions=`portrait,noborder,height=5in,width=5in`);Q:=graph({0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20},[[0,1],[0,2],[0,3],[0,4],[0,5],[1,6],[1,7],[1,8],[2,9],[2,10],[2,11],[3,12],[3,13],[3,14],[4,15],[4,16],[4,17],[5,18],[5,19],[5,20]]):draw(Concentric([0],[1,2,3,4,5]),Q);with(linalg):max(evalf(eigenvals(adjacency(G))));K:=adjacency(Q); DOESNT DO IT RIGHT!?Y:=K+transpose(K);max(evalf(eigenvals(Y)));rowdim(Y);coldim(Y);evalf(eigenvects(Y));40/21.;

restart:

with(linalg):w:=matrix(4,4); fitness matrixx:=vector(4); vector of gamete frequencieswmean:=multiply( transpose(x),w,x); mean fitness of the population W:=multiply(w,x); vector of induced fitnesseseta:=[-1, 1, 1, -1];d:=x[1]*x[4]-x[2]*x[3]; linkage disequiibriumfor i from 1 to 4 do f[i]:=(W[i]*x[i]+r*eta[i]*w[2,3]*d)/wmean; od; gamete frequencies in the next generation

S:=matrix(4,4); define stability matrixfor i from 1 to 4 do for j from 1 to 4 do S[i,j]:=diff(f[i],x[j]); define elements of S od;od;

S1000:=subs({x[1]=1,x[2]=0,x[3]=0,x[4]=0},evalm(S)); stability matrix at equilibrium (1,0,0,0)S0100:=subs({x[1]=0,x[2]=1,x[3]=0,x[4]=0},evalm(S)); stability matrix at equilibrium (0,1,0,0)S0010:=subs({x[1]=0,x[2]=0,x[3]=1,x[4]=0},evalm(S)); stability matrix at equilibrium (0,0,1,0)S0001:=subs({x[1]=0,x[2]=0,x[3]=0,x[4]=1},evalm(S)); stability matrix at equilibrium (0,0,0,1)

v:=matrix(3,3):for i from 1 to 3 do for j from 1 to 3 do v[i,j]:=alpha[i]+beta[j]; additive fitness od;od;evalm(v);w:=matrix(4,4, [ [ v[1,1], v[1,2], v[2,1], v[2,2] ], [ v[1,2], v[1,3], v[2,2], v[2,3] ], [ v[2,1], v[2,2], v[3,1], v[3,2] ], [ v[2,2], v[2,3], v[3,2], v[3,3]] ]);simplify(W[1]-W[2]-W[3]+W[4]); epistatic term is zerod:='d';Eq:={ x[1]=p[1]*p[2]+d, x[2]=p[1]*(1-p[2])-d, x[3]=(1-p[1])*p[2]-d, x[4]=(1-p[1])*(1-p[2])+d}; expressing x_i in terms of p1,p2 and DWbar:=simplify(subs(Eq ,wmean)); average fitness doesn't depend on LDdeltap1:=f[1]+f[2]-(x[1]+x[2]);Delta1:=factor(subs(Eq,deltap1));factor(subs({d=0},numer(Delta1)));deltap2:=f[1]+f[3]-(x[1]+x[3]):Delta2:=factor(subs(Eq,deltap2));factor(subs({d=0},numer(Delta2)));

v:=matrix(3,3):for i from 1 to 3 do for j from 1 to 3 do v[i,j]:=alpha[i]*beta[j]; additive fitness od;od;evalm(v);w:=matrix(4,4, [ [ v[1,1], v[1,2], v[2,1], v[2,2] ], [ v[1,2], v[1,3], v[2,2], v[2,3] ], [ v[2,1], v[2,2], v[3,1], v[3,2] ], [ v[2,2], v[2,3], v[3,2], v[3,3]] ]);alpha[1]:=1-s: alpha[2]:=1: alpha[3]:=1-s: beta[1]:=1-s: beta[2]:=1: beta[3]:=1-s: special case Eq0:={x[1]=1/4,x[2]=1/4,x[3]=1/4,x[4]=1/4,d=0};for i from 1 to 4 do print(simplify( subs(Eq0,f[i]))); od; this shows that Eq0 is an equilibriumS0:=subs(Eq0,evalm(S));lambda:=eigenvals(S0);factor(1-lambda[1]);factor(1-lambda[3]); conditions for stability can be found from hereEq1:={x[1]=1/4+d,x[2]=1/4-d,x[3]=1/4-d,x[4]=1/4+d};subs(Eq1,f[1]-x[1]);factor(%); Su:=subs(Eq1, evalm(S));kappa:=eigenvals(Su);subs({D^2=(s^2-4*r)/(16*s^2), D^4=(s^2-4*r)^2/(16*s^2)^2},1-kappa[3]);factor(%);subs({D^2=(s^2-4*r)/(16*s^2), D^4=(s^2-4*r)^2/(16*s^2)^2},1-kappa[2]);factor(%);subs({D^2=(s^2-4*r)/(16*s^2), D^4=(s^2-4*r)^2/(16*s^2)^2},1-kappa[4]);factor(%);

v:=matrix(3,3,[[1-delta,1-beta,1-alpha],[1-gamma,1,1-gamma],[1-alpha,1-beta,1-delta]]);w:=matrix(4,4, [ [ v[1,1], v[1,2], v[2,1], v[2,2] ], [ v[1,2], v[1,3], v[2,2], v[2,3] ], [ v[2,1], v[2,2], v[3,1], v[3,2] ], [ v[2,2], v[2,3], v[3,2], v[3,3]] ]);d:='d';Eq0:={x[1]=1/4,x[2]=1/4,x[3]=1/4,x[4]=1/4,d=0};factor(subs(Eq0,f[1]+f[2]-x[1]-x[2]));S0:=subs(Eq0,evalm(S));lambda:=eigenvals(S0);Eq1:={x[1]=1/4+D,x[2]=1/4-D,x[3]=1/4-D,x[4]=1/4+D};subs(Eq1,f[1]-x[1]);factor(%); Su:=subs( Eq1, evalm(S));kappa:=eigenvals(Su);

Weak selection appoximation for a 3-locus 2-allele haploid population with parwise epistasis q1:=1-p1: q2:=1-p2: q3:=1-p3:f1:=p1*q1*(a1+2*(b12*p2+b13*p3)); # dynamic equationsf2:=p2*q2*(a2+2*(b21*p1+b23*p3));f3:=p3*q3*(a3+2*(b31*p1+b32*p2));poly:=solve( {a1+2*(b12*p2+b13*p3)=0, a2+2*(b21*p1+b23*p3)=0, a3+2*(b31*p1+b32*p2)=0}, {p1,p2,p3}); # solving for a polymorphic equilibriumwith(linalg): S:=matrix(3,3,[[ diff(f1,p1), diff(f1,p2), diff(f1,p3)],[diff(f2,p1), diff(f2,p2), diff(f2,p3)],[diff(f3,p1), diff(f3,p2), diff(f3,p3)]]); # stability matrixSpoly:=subs(poly,evalm(S)); # stability matrix at the polymorphic equilibriumsubs( {a1=a,a2=a,a3=a,b12=b,b21=b,b13=b,b31=b,b23=b,b32=b}, evalm(Spoly)); # symmetric caseeigenvals(%);subs( {a1=a,a2=a,a3=a,b12=b1*b2,b21=b1*b2,b13=b1*b3,b31=b1*b3,b23=b2*b3,b32=b2*b3}, evalm(Spoly)); # symmetric caseeigenvals(%);