restart:w[A]:= (w[AA]*p+w[Aa]*q); w[a]:= (w[Aa]*p+w[aa]*q); induced fitnesses of A and awmean:=(w[AA]*p^2+2*w[Aa]*p*q+w[aa]*q^2 ); average fitness of the populationq:=1-p;ps:=w[A]*p/wmean: qs:=w[a]*q/wmean: allele frequencies after selectionpm:=(1-mu)*ps+mu*qs: frequency of A after mutationDelta:=pm-p: change in p per generation w[AA]:=1-s: w[Aa]:=1-h : w[aa]:=1: genotype fitnessesDelta;factor(%);Eq:=numer(%);subs(s=0,Eq); symmetric underdominancefactor(%);sol:=solve(%=0,p);assume(h>0): simplify(taylor(sol[2],mu=0,2));simplify(taylor(sol[3],mu=0,2));subs(h=s/2,Eq); additive fitnessessol:=solve(%=0,p);assume(s>0): simplify(taylor(sol[1],mu=0,2));assume(s>0): simplify(taylor(sol[2],mu=0,2));subs(h=0,Eq); recessive allelef:=collect(%,p);factor(subs(mu=0,%));p:=sum(P[i]*mu^(epsilon*i),i=0..2);P[0]:=0:combine(expand(f));collect(%,mu,factor);subs(p=sqrt(mu)/sqrt(s)+P[2]*mu,f);combine(expand(%));mu^(3/2)/s^(1/2)+2*mu^(3/2)*s^(1/2)*P[2]; solve(%,P[2]);

Changes in allele frequency as modeled by binomial samplingrestart:with(stats):N:=50: p:=.4: T:=100: X:=array[0..T]:for i from 1 to T don:=0: for j from 1 to N do x:=random[uniform](1): if x<p then n:=n+1: fi od: p:=n/N: X[i]:=p; od:F:=seq([i,X[i]],i=0..T): plot({F}, style=point, labels=[ `generation number`,`freq`], title=`Genetic Drift`,color=[red,blue,green],titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);

restart:with(DEtools): General solution for the probability of fixationeq:=m(p)*diff(u(p),p)+v(p)*diff(u(p),p$2)/2;

NiM+SSNlcUc2IiwmKiYtSSJtR0YlNiNJInBHRiUiIiItSSVkaWZmR0kqcHJvdGVjdGVkR0YvNiQtSSJ1R0YlRipGK0YsRiwqJi1JInZHRiVGKkYsLUYuNiRGMS1JIiRHRi82JEYrIiIjRiwjRixGOw==

dsolve(eq=0,u(p));

NiMvLUkidUc2IjYjSSJwR0YmLCZJJF9DMUdGJiIiIiomLUkkSW50RzYkSSpwcm90ZWN0ZWRHRjBJKF9zeXNsaWJHRiY2JC1JJGV4cEdGLzYjLUYuNiQsJComLUkibUdGJkYnRistSSJ2R0YmRichIiIhIiNGKEYoRitJJF9DMkdGJkYrRis=

F:=exp(-2*int(m(p)/v(p),p));

NiM+SSJGRzYiLUkkZXhwRzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IywkLUkkaW50R0YoNiQqJi1JIm1HRiU2I0kicEdGJSIiIi1JInZHRiVGMyEiIkY0ISIj

U:=int(F,p=0..x)/int(F,p=0..1);

NiM+SSJVRzYiKiYtSSRpbnRHNiRJKnByb3RlY3RlZEdGKkkoX3N5c2xpYkdGJTYkLUkkZXhwR0YpNiMsJC1GKDYkKiYtSSJtR0YlNiNJInBHRiUiIiItSSJ2R0YlRjYhIiJGNyEiIy9GNzsiIiFJInhHRiVGOC1GKDYkRi0vRjc7Rj9GOEY7

Additive fitnessesm(p):=s*p*(1-p); v(p):=p*(1-p)/(2*N);

NiM+LUkidkc2IjYjSSJwR0YmLCQqKEYoIiIiLCZGK0YrRighIiJGK0kiTkdGJkYtI0YrIiIj

U;F:=subs(x=1/(2*N),%); assuming there is a single mutant initiallyplot( {subs(N=25,log10(F)),subs(N=50,log10(F)),subs(N=100,log10(F))},s=-0.1..0.1,y=-17..0,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12],thickness=2,labels=[ `s`,`log10(u)`], title=`Fixation probability`); Underdominant selection

restart:f:=p*q*(w[A]-w[a])/wmean;

NiM+SSJmRzYiKipJInBHRiUiIiJJInFHRiVGKCwmJkkid0dGJTYjSSJBR0YlRigmRiw2I0kiYUdGJSEiIkYoSSZ3bWVhbkdGJUYy

w[A]:=w[AA]*p+w[Aa]*q;

NiM+Jkkid0c2IjYjSSJBR0YmLCYqJiZGJTYjSSNBQUdGJiIiIkkicEdGJkYuRi4qJiZGJTYjSSNBYUdGJkYuSSJxR0YmRi5GLg==

w[a]:=w[Aa]*p+w[aa]*q;

NiM+Jkkid0c2IjYjSSJhR0YmLCYqJiZGJTYjSSNBYUdGJiIiIkkicEdGJkYuRi4qJiZGJTYjSSNhYUdGJkYuSSJxR0YmRi5GLg==

wmean:=w[A]*p+w[a]*q;

NiM+SSZ3bWVhbkc2IiwmKiYsJiomJkkid0dGJTYjSSNBQUdGJSIiIkkicEdGJUYuRi4qJiZGKzYjSSNBYUdGJUYuSSJxR0YlRi5GLkYuRi9GLkYuKiYsJiomRjFGLkYvRi5GLiomJkYrNiNJI2FhR0YlRi5GNEYuRi5GLkY0Ri5GLg==

q:=1-p;

NiM+SSJxRzYiLCYiIiJGJ0kicEdGJSEiIg==

Delta[p]:=factor(f);

NiM+JkkmRGVsdGFHNiI2I0kicEdGJiwkKipGKCIiIiwmISIiRitGKEYrRissLComJkkid0dGJjYjSSNBQUdGJkYrRihGK0YtJkYxNiNJI0FhR0YmRi0qJkY0RitGKEYrIiIjJkYxNiNJI2FhR0YmRisqJkY5RitGKEYrRi1GKywuKiZGMEYrRihGOEYtRjchIiMqJkY0RitGKEY4RjhGOUYtRjxGOComRjlGK0YoRjhGLUYtRi0=

factor(subs({w[AA]=1,w[Aa]=1-s,w[aa]=1},Delta[p]));

NiMsJCosSSJwRzYiIiIiLCYhIiJGJ0YlRidGJ0kic0dGJkYnLCZGKUYnRiUiIiNGJywoKiZGJUYnRipGJyEiIyomRiVGLEYqRidGLEYnRidGKUYp

numer(%)/collect(denom(%),s,factor);

NiMsJCosSSJwRzYiIiIiLCYhIiJGJ0YlRidGJ0kic0dGJkYnLCZGKUYnRiUiIiNGJywmRidGJyooRipGJ0YoRidGJUYnRixGKUYp

if s<<1, the denominator (i.e. the average fitness of the population) is approximately 1. Unequal fitness of homozygotesfactor(subs({w[AA]=1,w[Aa]=1-s,w[aa]=1+h},Delta[p]));

NiMsJCoqSSJwRzYiIiIiLCYhIiJGJ0YlRidGJywqSSJzR0YmRikqJkYlRidGK0YnIiIjSSJoR0YmRikqJkYlRidGLkYnRidGJywuRiwhIiMqJkYlRi1GK0YnRi1GJ0YnRi5GJ0YvRjEqJkYlRi1GLkYnRidGKUYp

m(p):=s*p*(1-p)*(2*p-1+e*(p-1)); different peaks

NiM+LUkibUc2IjYjSSJwR0YmKipJInNHRiYiIiJGKEYrLCZGK0YrRighIiJGKywoRigiIiNGLUYrKiZJImVHRiZGKywmRihGK0YtRitGK0YrRis=

m(p):=s*p*(1-p)*(2*p-1); equal peaks

NiM+LUkibUc2IjYjSSJwR0YmKipJInNHRiYiIiJGKEYrLCZGK0YrRighIiJGKywmRigiIiNGLUYrRis=

U;

NiMsJComLCYtSSRlcmZHNiRJKnByb3RlY3RlZEdGKUkoX3N5c2xpYkc2IjYjKihJInNHRisiIiJJIk5HRitGLyomRi5GL0YwRi8jISIiIiIjRi8tRic2IyoqRi5GL0YwRi8sJkkieEdGK0Y0RjNGL0YvRjFGMkYvRi9GJkYzI0YvRjQ=

subs(x=1/(2*N),%);

NiMsJComLCYtSSRlcmZHNiRJKnByb3RlY3RlZEdGKUkoX3N5c2xpYkc2IjYjKihJInNHRisiIiJJIk5HRitGLyomRi5GL0YwRi8jISIiIiIjRi8tRic2IyoqRi5GL0YwRi8sJiokRjBGM0YvRjNGL0YvRjFGMkYvRi9GJkYzI0YvRjQ=

X:=(erf(Y*(e-1))+erf(Y))/(2*erf(Y));

NiM+SSJYRzYiLCQqJiwmLUkkZXJmRzYkSSpwcm90ZWN0ZWRHRixJKF9zeXNsaWJHRiU2IyomSSJZR0YlIiIiLCZJImVHRiVGMSEiIkYxRjFGMS1GKjYjRjBGMUYxRjVGNCNGMSIiIw==

series(X,e=0,2);

NiMrJ0kiZUc2IioqSSNQaUdJKnByb3RlY3RlZEdGKCMhIiIiIiMtSSRleHBHNiRGKEkoX3N5c2xpYkdGJTYjLCQqJEkiWUdGJUYrRioiIiJGM0Y0LUkkZXJmR0YuNiNGM0YqRjQtSSJPR0YoNiNGNEYr

subs({Y=sqrt(S),e=1/N},exp(-Y^2)*Y/(sqrt(Pi)*erf(Y))*e);

NiMqLC1JJGV4cEc2JEkqcHJvdGVjdGVkR0YnSShfc3lzbGliRzYiNiMsJEkiU0dGKSEiIiIiIkYsI0YuIiIjSSNQaUdGJyNGLUYwLUkkZXJmR0YmNiMqJEYsRi9GLUkiTkdGKUYt

restart:with(DEtools): General solution for the probability of fixationm(p):=s*p*(1-p); v(p):=p*(1-p)/(2*N);

NiM+LUkibUc2IjYjSSJwR0YmKihJInNHRiYiIiJGKEYrLCZGK0YrRighIiJGKw==NiM+LUkidkc2IjYjSSJwR0YmLCQqKEYoIiIiLCZGK0YrRighIiJGK0kiTkdGJkYtI0YrIiIj

eq:=1+m(p)*diff(T(p),p)+v(p)*diff(T(p),p$2)/2;

NiM+SSNlcUc2IiwoIiIiRicqKkkic0dGJUYnSSJwR0YlRicsJkYnRidGKiEiIkYnLUklZGlmZkdJKnByb3RlY3RlZEdGLzYkLUkiVEdGJTYjRipGKkYnRicqKkYqRidGK0YnSSJOR0YlRiwtRi42JEYxLUkiJEdGLzYkRioiIiNGJyNGJyIiJQ==

dsolve(eq=0,T(p));

NiMvLUkiVEc2IjYjSSJwR0YmLCYtSSRJbnRHNiRJKnByb3RlY3RlZEdGLUkoX3N5c2xpYkdGJjYkLCgqKC1JJGV4cEdGLDYjLCQqKEkic0dGJiIiIkYoRjhJIk5HRiZGOCEiJUY4RjlGOC1JI0VpR0YsNiRGOEY1RjgiIiUqKEY5RjgtRjw2JEY4LCYqJkY3RjhGOUY4Rj5GNkY6RjgtRjM2IywkKihGN0Y4RjlGOCwmISIiRjhGKEY4RjhGOkY4RjoqJkYyRjhJJF9DMUdGJkY4RjhGKEY4SSRfQzJHRiZGOA==

restart;P0(p):=int(F(y),y=p..1)/int(F(y),y=0..1);

NiM+LUkjUDBHNiI2I0kicEdGJiomLUkkaW50RzYkSSpwcm90ZWN0ZWRHRi1JKF9zeXNsaWJHRiY2JC1JIkZHRiY2I0kieUdGJi9GMztGKCIiIkY2LUYrNiRGMC9GMzsiIiFGNiEiIg==

P1(p):=int(F(y),y=0..p)/int(F(y),y=0..1);

NiM+LUkjUDFHNiI2I0kicEdGJiomLUkkaW50RzYkSSpwcm90ZWN0ZWRHRi1JKF9zeXNsaWJHRiY2JC1JIkZHRiY2I0kieUdGJi9GMzsiIiFGKCIiIi1GKzYkRjAvRjM7RjZGNyEiIg==

F(y):=exp(-2*int(m(y)/v(y),y));

NiM+LUkiRkc2IjYjSSJ5R0YmLUkkZXhwRzYkSSpwcm90ZWN0ZWRHRixJKF9zeXNsaWJHRiY2IywkLUkkaW50R0YrNiQqJi1JIm1HRiZGJyIiIi1JInZHRiZGJyEiIkYoISIj

U:=int(F,p=0..x)/int(F,p=0..1);