MATHEMATICAL EVOLUTIONARY THEORY - MATH 583plotsetup(default):plotsetup(window):

k:=5: t_max:=100: number of alleles and maximum timef[1]:=.2: f[2]:=.2: f[3]:=.2: f[4]:=.2: f[5]:=.2: initial allele frequenciesw[1]:=.92: w[2]:=.94: w[3]:=.96: w[4]:=.98: w[5]:=1: fitnesseswmean:=sum( 'f[i]*w[i]^t', 'i'=1..k): Wmean:=sum( 'f[i]*w[i]^(t+1)', 'i'=1..k):wbar:= Wmean/wmean: mean fitness of the population in generation tfor i from 1 to k do p[i]:=f[i]*w[i]^t/wmean od: compute allele frequencies at generation tplot( {p[1],p[2],p[3],p[4],p[5],wbar}, t=0..t_max, title=`Allele frequencies and mean fitness\n in 1-locus k-allele haploid population with constant fitnesses`, labels=[generation_number,p_i],titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]); plot

A system of two linear ODE: x' = ax + by, y' = cx + dy.Its solutions depend on the eigenvalues of the stability matrix [ a b] S = [ c d].S := linalg[matrix](2,2,[a,b,c,d]); define a 2x2 matrix S

NiM+SSJTRzYiLUknbWF0cml4RzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IzckNyRJImFHRiVJImJHRiU3JEkiY0dGJUkiZEdGJQ==

eq := [ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)]; with(DEtools):

NiM+SSNlcUc2IjckLy1JJWRpZmZHSSpwcm90ZWN0ZWRHRio2JC1JInhHRiU2I0kidEdGJUYvLCYqJkkiYUdGJSIiIkYsRjNGMyomSSJiR0YlRjMtSSJ5R0YlRi5GM0YzLy1GKTYkRjZGLywmKiZJImNHRiVGM0YsRjNGMyomSSJkR0YlRjNGNkYzRjM=

Several cases to consider.Both eigenvalues real and positive.a := 1: b := 0: c := 1: d := 1.5: evalf(Eigenvals(S));

NiMtSSd2ZWN0b3JHNiI2IzckJCIjOiEiIiQiIiIiIiE=

phaseportrait([ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)],[x(t),y(t)],-2..1,{[0,.5,.5],[0,.5,-.5],[0,-.5,-.5],[0,-.5,.5]},stepsize=.2,x=-1..1,y=-1..1,title=`Unstable node\n (eigenvalues real and positive)`,arrows=THICK,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);Both eigenvalues real and negativea := -.9: b := 0: c := 1: d := -1.5: evalf(Eigenvals(S));

NiMtSSd2ZWN0b3JHNiI2IzckJCEjOiEiIiQhIipGKg==

phaseportrait([ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)],[x(t),y(t)],-2..1,{[0,.5,.5],[0,.5,-.5],[0,-.5,-.5],[0,-.5,.5]},stepsize=.2,x=-1..1,y=-1..1,title=`Stable node\n (eigenvlues real and negative)`,arrows=THICK,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);One positive eigenvalue and one negative eigenvaluea := 1: b := 0: c := 1: d := -1.5: evalf(Eigenvals(S));

NiMtSSd2ZWN0b3JHNiI2IzckJCEjOiEiIiQiIiIiIiE=

phaseportrait([ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)],[x(t),y(t)],-2..1,{[0,.5,.5],[0,.5,-.5],[0,-.5,-.5],[0,-.5,.5]}, stepsize=.2,x=-1..1,y=-1..1, title=`Saddle\n (one negative, one positive eigenvalue)`, arrows=THICK, titlefont=[TIMES,BOLD,12], font=[TIMES,BOLD,12]);Complex eigenvalues, Re(l)>0a := .1: b := 1: c := -1: d := .5: evalf(Eigenvals(S));

NiMtSSd2ZWN0b3JHNiI2IzckXiQkIisrKysrSSEjNSQiK3IqZXp6KkYrXiRGKSQhK3IqZXp6KkYr

phaseportrait([ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)],[x(t),y(t)],-2..1,{[0,.5,.5],[0,.5,-.5],[0,-.5,-.5],[0,-.5,.5]},stepsize=.2,x=-1..1,y=-1..1,title=`Unstable focus\n (complex eigenvalues, Re(l)>0)`,arrows=THICK,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);Complex eigenvalues, Re(l)<0a :=- .5: b := 1: c := -1: d := -.5: evalf(Eigenvals(S));

NiMtSSd2ZWN0b3JHNiI2IzckXiQkIisrKysrSSEjNSQiK3IqZXp6KkYrXiRGKSQhK3IqZXp6KkYr

phaseportrait([ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)],[x(t),y(t)],-2..1,{[0,.5,.5],[0,.5,-.5],[0,-.5,-.5],[0,-.5,.5]},stepsize=.2,x=-1..1,y=-1..1,title=`Stable focus\n (complex eigenvalues, Re(l)<0)`,arrows=THICK,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);Complex eigenvalues, Re(l)=0a :=0: b := 1: c := -1: d := .0: evalf(Eigenvals(S)); phaseportrait([ diff(x(t),t)=a*x(t)+b*y(t), diff(y(t),t) = c*x(t)+d*y(t)],[x(t),y(t)],-2..1,{[0,.5,.5],[0,.5,-.5],[0,-.5,-.5],[0,-.5,.5]},stepsize=.2,x=-1..1,y=-1..1,title=`center (complex, Re(l)=0)`,arrows=THICK,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);?plot

The following system will be used as an examplerestart: f1:=x*(a-b*x-c*y); specify f1(x,y)

NiM+SSNmMUc2IiomSSJ4R0YlIiIiLChJImFHRiVGKComSSJiR0YlRihGJ0YoISIiKiZJImNHRiVGKEkieUdGJUYoRi1GKA==

f2:=y*(d*x-e-f*y); specify f2(x,y)

NiM+SSNmMkc2IiomSSJ5R0YlIiIiLCgqJkkiZEdGJUYoSSJ4R0YlRihGKEkiZUdGJSEiIiomSSJmR0YlRihGJ0YoRi5GKA==

with(linalg): load libraries

Warning, the protected names norm and trace have been redefined and unprotected

sols:=solve({f1,f2},{x,y});

NiM+SSVzb2xzRzYiNiY8JC9JInhHRiUiIiEvSSJ5R0YlRio8JEYoL0YsLCQqJkkiZUdGJSIiIkkiZkdGJSEiIkY0PCRGKy9GKSomSSJhR0YlRjJJImJHRiVGNDwkL0YpKiYsJiomSSJjR0YlRjJGMUYyRjIqJkYzRjJGOEYyRjJGMiwmKiZGP0YySSJkR0YlRjJGMiomRjNGMkY5RjJGMkY0L0YsLCQqJiwmKiZGOEYyRkNGMkY0KiZGOUYyRjFGMkYyRjJGQUY0RjQ=

S:=array([[diff(f1,x),diff(f1,y)],[diff(f2,x),diff(f2,y)]]); define stability matrix

NiM+SSJTRzYiLUknbWF0cml4RzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IzckNyQsKEkiYUdGJSIiIiomSSJiR0YlRjBJInhHRiVGMCEiIyomSSJjR0YlRjBJInlHRiVGMCEiIiwkKiZGM0YwRjZGMEY4NyQqJkY3RjBJImRHRiVGMCwoKiZGPUYwRjNGMEYwSSJlR0YlRjgqJkkiZkdGJUYwRjdGMEY0

S0:=subs( [x=0,y=0],evalm(S)); substitute equilibrium values (0,0) in S

NiM+SSNTMEc2Ii1JJ21hdHJpeEc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JDckSSJhR0YlIiIhNyRGLywkSSJlR0YlISIi

lambda0:=eigenvals(evalm(S0)); find eigenvalues

NiM+SShsYW1iZGEwRzYiNiRJImFHRiUsJEkiZUdGJSEiIg==

S1:=subs( [x=a/b,y=0], evalm(S)); substitute equilibrium values (a/b,0) in S

NiM+SSNTMUc2Ii1JJ21hdHJpeEc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JDckLCRJImFHRiUhIiIsJCooRi8iIiJJImJHRiVGMEkiY0dGJUYzRjA3JCIiISwmKihJImRHRiVGM0YvRjNGNEYwRjNJImVHRiVGMA==

lambda1:=eigenvals( evalm(S1)); find eigevnalues

NiM+SShsYW1iZGExRzYiNiQsJEkiYUdGJSEiIiwkKiYsJiomRigiIiJJImRHRiVGLkYpKiZJImJHRiVGLkkiZUdGJUYuRi5GLkYxRilGKQ==

nontriv:=solve( {a-b*x-c*y=0,d*x-e-f*y=0},{x,y}); find nontrivial equilibrium

NiM+SShub250cml2RzYiPCQvSSJ4R0YlKiYsJiomSSJjR0YlIiIiSSJlR0YlRi1GLSomSSJmR0YlRi1JImFHRiVGLUYtRi0sJiomRixGLUkiZEdGJUYtRi0qJkYwRi1JImJHRiVGLUYtISIiL0kieUdGJSwkKiYsJiomRjFGLUY0Ri1GNyomRjZGLUYuRi1GLUYtRjJGN0Y3

S2:=subs( nontriv, evalm(S)); substitute into S

NiM+SSNTMkc2Ii1JJ21hdHJpeEc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JDckLChJImFHRiUiIiIqKEkiYkdGJUYwLCYqJkkiY0dGJUYwSSJlR0YlRjBGMComSSJmR0YlRjBGL0YwRjBGMCwmKiZGNUYwSSJkR0YlRjBGMComRjhGMEYyRjBGMCEiIiEiIyooRjVGMCwmKiZGL0YwRjtGMEY9KiZGMkYwRjZGMEYwRjBGOUY9RjAsJCooRjNGMEY5Rj1GNUYwRj03JCwkKihGQEYwRjlGPUY7RjBGPSwoKihGO0YwRjNGMEY5Rj1GMEY2Rj0qKEY4RjBGQEYwRjlGPSIiIw==

lambda2:=eigenvals( evalm(S2)); find eigenvalues

NiM+SShsYW1iZGEyRzYiNiQsJComLCYqJkkiY0dGJSIiIkkiZEdGJUYsRiwqJkkiZkdGJUYsSSJiR0YlRixGLCEiIiwsKihGL0YsRjBGLEkiYUdGJUYsRjEqKEYvRixGMEYsSSJlR0YlRixGLCooRi9GLEYtRixGNEYsRjEqKEYwRixGK0YsRjZGLEYxKiQsPCooRi8iIiNGMEY8RjRGPEYsKipGNEYsRi9GPEYwRjxGNkYsRjwqKkY0RjxGL0Y8RjBGLEYtRiwhIiMqLEYvRixGMEY8RjRGLEYrRixGNkYsRjwqKEYvRjxGMEY8RjZGPEYsKixGL0Y8RjBGLEY2RixGLUYsRjRGLEY/KipGMEY8RitGLEY2RjxGL0YsRjwqKEYvRjxGLUY8RjRGPEYsKi5GL0YsRi1GLEY0RixGMEYsRitGLEY2RixGPCooRjBGPEYrRjxGNkY8RiwqKkYrRixGLUY8Ri9GLEY0RjwhIiUqKkYrRjxGLUY8RjZGLEY0RixGSCoqRitGPEYtRixGNkY8RjBGLCIiJSNGLEY8RixGLEZMLCQqJkYpRjEsLEYzRjFGNUYsRjdGMUY4RjFGOUYxRixGTA==

det(S2); find determinant

NiMsJCooLCYqJkkiY0c2IiIiIkkiZUdGKEYpRikqJkkiZkdGKEYpSSJhR0YoRilGKUYpLCYqJkYtRilJImRHRihGKSEiIiomSSJiR0YoRilGKkYpRilGKSwmKiZGJ0YpRjBGKUYpKiZGLEYpRjNGKUYpRjFGMQ==

factor( trace(S2)); find trace

NiMqJiwqKihJImZHNiIiIiJJImJHRidGKEkiYUdGJ0YoISIiKihGJkYoRilGKEkiZUdGJ0YoRigqKEYmRihJImRHRidGKEYqRihGKyooRilGKEkiY0dGJ0YoRi1GKEYrRigsJiomRjFGKEYvRihGKComRiZGKEYpRihGKEYr

g:=x-> x * (wAA*x+wAa*(1-x)) / (wAA*x^2 +2*wAa*x*(1-x)+waa*(1-x)^2 ): dynamic equationT:=1000: Z:=array[0..T]: Z[0]:= .00001: wAA:=1: wAa:=.95: waa:=.9: parameters and initial conditionfor i from 1 to T do y:=g(Z[i-1]): Z[i]:=y od: iterate dynamic equationplot( [seq([i,Z[i]],i=0..T)], style=point, labels=[ `generation number `, `x` ], title=`1-locus 2-allele model`,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]); plot graphrestart: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==

factor(f);

NiMsJCoqSSJwRzYiIiIiLCYhIiJGJ0YlRidGJywsKiYmSSJ3R0YmNiNJI0FBR0YmRidGJUYnRikmRi02I0kjQWFHRiZGKSomRjBGJ0YlRiciIiMmRi02I0kjYWFHRiZGJyomRjVGJ0YlRidGKUYnLC4qJkYsRidGJUY0RilGMyEiIyomRjBGJ0YlRjRGNEY1RilGOEY0KiZGNUYnRiVGNEYpRilGKQ==

sols:=solve(f,p);

NiM+SSVzb2xzRzYiNiUiIiEiIiIqJiwmJkkid0dGJTYjSSNBYUdGJSEiIiZGLDYjSSNhYUdGJUYoRigsKCZGLDYjSSNBQUdGJUYoRishIiNGMEYoRi8=

S:=diff(f,p);

NiM+SSJTRzYiLCoqKCwmIiIiRilJInBHRiUhIiJGKSwqKiYmSSJ3R0YlNiNJI0FBR0YlRilGKkYpRikqJiZGLzYjSSNBYUdGJUYpRihGKUYpKiZGM0YpRipGKUYrKiYmRi82I0kjYWFHRiVGKUYoRilGK0YpLCYqJiwmRi1GKUYyRilGKUYqRilGKSomLCZGNkYpRjdGKUYpRihGKUYpRitGKSooRipGKUYsRilGO0YrRisqKkYqRilGKEYpLChGLkYpRjMhIiNGOEYpRilGO0YrRikqLEYqRilGKEYpRixGKUY7RkMsLiomLCZGLkYpRjNGK0YpRipGKUYpRi1GKUYyRikqJiwmRjNGKUY4RitGKUYoRilGKUY2RitGN0YrRilGKw==

subs(p=sols[1],S);

NiMqJiwmJkkid0c2IjYjSSNBYUdGJyIiIiZGJjYjSSNhYUdGJyEiIkYqRitGLg==

subs(p=sols[2],S);

NiMsJComLCYmSSJ3RzYiNiNJI0FBR0YoIiIiJkYnNiNJI0FhR0YoISIiRitGJkYvRi8=

subs(p=sols[3],S);

NiMsKiooLCYiIiJGJiomLCYmSSJ3RzYiNiNJI0FhR0YrISIiJkYqNiNJI2FhR0YrRiZGJiwoJkYqNiNJI0FBR0YrRiZGKSEiI0YvRiZGLkYuRiYsKiooRjNGJkYoRiZGMkYuRiYqJkYpRiZGJUYmRiYqKEYpRiZGKEYmRjJGLkYuKiZGL0YmRiVGJkYuRiYsJiooLCZGOEYmRjlGJkYmRihGJkYyRi5GJiomLCZGOkYmRjtGJkYmRiVGJkYmRi5GJioqRihGJkYyRi5GN0YmRjxGLkYuKihGKEYmRiVGJkY8Ri5GJiouRihGJkYyRi5GJUYmRjdGJkY8RjYsLiooLCZGM0YmRilGLkYmRihGJkYyRi5GJkY4RiZGOUYmKiYsJkYpRiZGL0YuRiZGJUYmRiZGOkYuRjtGLkYmRi4=

simplify(%);

NiMsJCooLCYmSSJ3RzYiNiNJI0FhR0YoISIiJkYnNiNJI2FhR0YoIiIiRi8sJiZGJzYjSSNBQUdGKEYrRiZGL0YvLCYqJkYsRi9GMUYvRi8qJEYmIiIjRitGK0Yr

restart:Two fitness matrices in one-locus three-allele systems illustrating that heterosis is neither a mathematically necessary nor sufficient condition for a stable feasible equilibrium with more than two alleles (Lewontin et al., 1978, 88, 149-170).with(linalg):z:=1-x-y:w:=array(1..3,1..3,[[.66,.75,.89],[.75,.28,.77],[.89,.77,.61]]);W1:=w[1,1]*x+w[1,2]*y+w[1,3]*z: induced fitness of allele 1W2:=w[2,1]*x+w[2,2]*y+w[2,3]*z: induced fitness of allele 2W3:=w[3,1]*x+w[3,2]*y+w[3,3]*z: induced fitness of allele 3solve({W2-W1=0,W3-W1=0},{x,y}); no feasible polymorphic solutionw:=array(1..3,1..3,[[.2358,.8457,.7482],[.8457,.1837,.3927],[.7482,.3927,.3964]]);W1:=w[1,1]*x+w[1,2]*y+w[1,3]*z: induced fitness of allele1 W2:=w[2,1]*x+w[2,2]*y+w[2,3]*z: induced fitness of allele 2W3:=w[3,1]*x+w[3,2]*y+w[3,3]*z: induced fitness of allele 3solve({W2-W1=0,W3-W1=0},{x,y}); polymorphic equilibrium is feasibleeigenvals(w); and is stable

restart:g:=x-> 1 / (1 + epsilon * x^2 ); 1-locus 2-allele model

NiM+SSJnRzYiZio2I0kieEdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlKiQsJiIiIkYuKiZJKGVwc2lsb25HRiVGLjkkIiIjRi4hIiJGJUYlRiU=

T:=200: Z:=array[0..T]: Z[0]:= .4: epsilon:=4.1: parameters and initial conditionfor i from 1 to T do y:=g(Z[i-1]): Z[i]:=y od: iterate dynamic equationplot( [seq([i,Z[i]],i=0..T)], style=point, labels=[ `generation number `, `x` ], title=`1-locus 2-allele model`,titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]); plot

LSUlUExPVEc2KC0lJ0NVUlZFU0c2JTdldzckJCIiIUYrJCIzQysrKysrKytTISM9NyQkIiIiRiskIjMlKSoqKioqSE1aJ1FnRi43JCQiIiNGKyQiMykqKioqKioqelYqeSslRi43JCQiIiRGKyQiMyEqKioqKip6ajUjSGdGLjckJCIiJUYrJCIzJSoqKioqKipRJjNhLCVGLjckJCIiJkYrJCIzLisrKyJmUS0tJ0YuNyQkIiInRiskIjM6KysrRHpjQVNGLjckJCIiKEYrJCIzTSsrKypRKnA2Z0YuNyQkIiIpRiskIjNAKysrOGRSSFNGLjckJCIiKkYrJCIzIikqKioqKnA2a04rJ0YuNyQkIiM1RiskIjMiKioqKioqNDE4Zi4lRi43JCQiIzZGKyQiM20qKioqKkhVMWUqZkYuNyQkIiM3RiskIjMmKioqKioqSExSQC8lRi43JCQiIzhGKyQiM0srKytTQVMpKWZGLjckJCIjOUYrJCIzQCsrKyFHIzRbU0YuNyQkIiM6RiskIjNaKysrTCZIOClmRi43JCQiIztGKyQiMyIpKioqKio0QSl5YFNGLjckJCIjPEYrJCIzNisrKyU0b1goZkYuNyQkIiM9RiskIjMmKSoqKioqNDxVI2ZTRi43JCQiIz5GKyQiM0QrKytRJCo0b2ZGLjckJCIjP0YrJCIzJioqKioqKkgqellrU0YuNyQkIiNARiskIjNWKysrXmghPidmRi43JCQiI0FGKyQiMzErKyslXHklcFNGLjckJCIjQkYrJCIzKCkqKioqKlJ1c2YmZkYuNyQkIiNDRiskIjMtKysrRmJHdVNGLjckJCIjREYrJCIzLSsrK0FYR11mRi43JCQiI0VGKyQiMyMpKioqKip6MysqeVNGLjckJCIjRkYrJCIzUSsrK256I1slZkYuNyQkIiNHRiskIjM3KysrJFFLTDMlRi43JCQiI0hGKyQiMzgrKytdMGZSZkYuNyQkIiNJRiskIjMxKysrMz5mKDMlRi43JCQiI0pGKyQiM1srKytEMWNNZkYuNyQkIiNLRiskIjNBKysrOHZvIjQlRi43JCQiI0xGKyQiM0wrKytrdHNIZkYuNyQkIiNNRiskIjM3KysrV3VpJjQlRi43JCQiI05GKyQiMzIrKysibyEzRGZGLjckJCIjT0YrJCIzOysrKzMlPiUqNCVGLjckJCIjUEYrJCIzYCoqKioqejs2MSNmRi43JCQiI1FGKyQiMz8rKys0MTIuVEYuNyQkIiNSRiskIjMqKioqKioqei41aiJmRi43JCQiI1NGKyQiMzorKysmeihlMVRGLjckJCIjVEYrJCIzVysrK24hcEAiZkYuNyQkIiNVRiskIjM2KysrJUh4KjRURi43JCQiI1ZGKyQiM0crKytmMD0zZkYuNyQkIiNXRiskIjMvKysrVl1DOFRGLjckJCIjWEYrJCIzQisrKy90TC9mRi43JCQiI1lGKyQiM3cqKioqKlxoJ1I7VEYuNyQkIiNaRiskIjNcKysrR0RqK2ZGLjckJCIjW0YrJCIzRisrK2lzVj5URi43JCQiI1xGKyQiMyEpKioqKioqbylmcSplRi43JCQiI11GKyQiM3YqKioqKio+PlBBVEYuNyQkIiNeRiskIjNxKioqKip6TThPKmVGLjckJCIjX0YrJCIzRCsrK1ZfP0RURi43JCQiI2BGKyQiM08rKytKdEchKmVGLjckJCIjYUYrJCIzKSoqKioqKlJpVHo3JUYuNyQkIiNiRiskIjMqKioqKioqKjRsMigpZUYuNyQkIiNjRiskIjMmKioqKioqej8mZUlURi43JCQiI2RGKyQiM1oqKioqKnAoZShSKWVGLjckJCIjZUYrJCIzeSoqKioqXCIqUko4JUYuNyQkIiNmRiskIjMlKSoqKioqenEhKTQpZUYuNyQkIiNnRiskIjMwKysrXCU0YzglRi43JCQiI2hGKyQiMzMrKytPbDN5ZUYuNyQkIiNpRiskIjMlKSoqKioqNEsoKno4JUYuNyQkIiNqRiskIjM2KysrTCIqR3ZlRi43JCQiI2tGKyQiMygpKioqKipSJm9JU1RGLjckJCIjbEYrJCIzJyoqKioqKjRdJWVzZUYuNyQkIiNtRiskIjM3KysrLDdhVVRGLjckJCIjbkYrJCIzWSoqKioqUiYpbypwZUYuNyQkIiNvRiskIjMlKSoqKioqPk4uWjklRi43JCQiI3BGKyQiMykqKioqKioqKWZRdSdlRi43JCQiI3FGKyQiM0MrKytZaHpZVEYuNyQkIiNyRiskIjNNKysrRi4qXCdlRi43JCQiI3NGKyQiMyopKioqKipIR0EpW1RGLjckJCIjdEYrJCIzdSoqKioqeiEzaWllRi43JCQiI3VGKyQiMysrKytGVnldVEYuNyQkIiN2RiskIjMoKioqKioqUidwS2dlRi43JCQiI3dGKyQiMzMrKytAWm9fVEYuNyQkIiN4RiskIjNWKysrd2U1ZWVGLjckJCIjeUYrJCIzISkqKioqKiopeURYOiVGLjckJCIjekYrJCIzXisrK2ZaJmYmZUYuNyQkIiMhKUYrJCIzRysrK1goNGo6JUYuNyQkIiMiKUYrJCIzQisrK3E0KFEmZUYuNyQkIiMjKUYrJCIzKikqKioqKkhxUSFlVEYuNyQkIiMkKUYrJCIzayoqKioqcCk+Jj0mZUYuNyQkIiMlKUYrJCIzeSoqKioqZm45KGZURi43JCQiIyYpRiskIjNaKysrLWEqKVxlRi43JCQiIycpRiskIjMiKSoqKioqemVSODslRi43JCQiIygpRiskIjMhKioqKioqcCMqKSp6JWVGLjckJCIjKSlGKyQiMyIpKioqKipcRjpIOyVGLjckJCIjKilGKyQiM3gqKioqKipwLjtZZUYuNyQkIiMhKkYrJCIzJCkqKioqKiopW1ZXOyVGLjckJCIjIipGKyQiM3QqKioqKnBreFYlZUYuNyQkIiMjKkYrJCIzJyoqKioqKnohZiNmOyVGLjckJCIjJCpGKyQiM10rKythKFtFJWVGLjckJCIjJSpGKyQiMyIqKioqKioqUVRPblRGLjckJCIjJipGKyQiM0ErKyshenI0JWVGLjckJCIjJypGKyQiMzErKys7KGYob1RGLjckJCIjKCpGKyQiMyIqKioqKioqR1xNUmVGLjckJCIjKSpGKyQiMzgrKyswVDZxVEYuNyQkIiMqKkYrJCIzUCsrK0Jrd1BlRi43JCQiJCsiRiskIjMhKioqKioqKj0oRzk8JUYuNyQkIiQsIkYrJCIzMysrKy5ZQk9lRi43JCQiJC0iRiskIjNGKysrL1xxc1RGLjckJCIkLiJGKyQiM3kqKioqKnoneXVNZUYuNyQkIiQvIkYrJCIzOisrK2JSJVI8JUYuNyQkIiQwIkYrJCIzOysrKyRvL0wkZUYuNyQkIiQxIkYrJCIzOysrKzlyOXZURi43JCQiJDIiRiskIjNjKioqKioqeU4hPiRlRi43JCQiJDMiRiskIjN5KioqKioqb2JKd1RGLjckJCIkNCJGKyQiMygpKioqKipIOVYwJGVGLjckJCIkNSJGKyQiMycpKioqKipSWV11PCVGLjckJCIkNiJGKyQiM00rKytDP0FIZUYuNyQkIiQ3IkYrJCIzdSoqKioqeipHYnlURi43JCQiJDgiRiskIjNJKysrRSpReiNlRi43JCQiJDkiRiskIjMlKioqKioqUiNSaXpURi43JCQiJDoiRiskIjNYKysrJ2YjcEVlRi43JCQiJDsiRiskIjM5KysrZVhtIT0lRi43JCQiJDwiRiskIjNoKioqKioqUT1bRGVGLjckJCIkPSJGKyQiMyMpKioqKipIeHY7PSVGLjckJCIkPiJGKyQiM0MrKysvYklDZUYuNyQkIiQ/IkYrJCIzKikqKioqKmZdZUU9JUYuNyQkIiRAIkYrJCIzbioqKioqXFtpSiNlRi43JCQiJEEiRiskIjMnKioqKioqemw4Tz0lRi43JCQiJEIiRiskIjNfKioqKipccl4/I2VGLjckJCIkQyJGKyQiMy0rKyssQGElPSVGLjckJCIkRCJGKyQiM2MrKytxQCg0I2VGLjckJCIkRSJGKyQiM0IrKytwWVcmPSVGLjckJCIkRiJGKyQiM08rKytvRyMqPmVGLjckJCIkRyJGKyQiM3oqKioqKlI7QWo9JUYuNyQkIiRIIkYrJCIzPSsrK2VHISo9ZUYuNyQkIiRJIkYrJCIzJyoqKioqKmZPdnI9JUYuNyQkIiRKIkYrJCIzcyoqKioqUkE2eiJlRi43JCQiJEsiRiskIjMhKSoqKioqPi0wISk9JUYuNyQkIiRMIkYrJCIzJCoqKioqKj40WnAiZUYuNyQkIiRNIkYrJCIzPysrK1o9IikpPSVGLjckJCIkTiJGKyQiM0YrKys1JzRnImVGLjckJCIkTyJGKyQiMzorKytTbGYqPSVGLjckJCIkUCJGKyQiM28qKioqKipmejQ6ZUYuNyQkIiRRIkYrJCIzcyoqKioqSHhmLj4lRi43JCQiJFIiRiskIjNrKioqKipcTjZVImVGLjckJCIkUyJGKyQiMyEqKioqKip6PS02PiVGLjckJCIkVCJGKyQiM2IrKytOIVxMImVGLjckJCIkVSJGKyQiMycpKioqKipIVEM9PiVGLjckJCIkViJGKyQiM24qKioqKipvLV43ZUYuNyQkIiRXIkYrJCIzPSsrK1hxXyM+JUYuNyQkIiRYIkYrJCIzXisrK1tWcDZlRi43JCQiJFkiRiskIjMjKioqKioqKnAxQCQ+JUYuNyQkIiRaIkYrJCIzIyoqKioqKioqZSs0ImVGLjckJCIkWyJGKyQiMyMqKioqKip6JWUoUT4lRi43JCQiJFwiRiskIjNeKioqKioqSCRHLCJlRi43JCQiJF0iRiskIjNHKysrQ0pfJT4lRi43JCQiJF4iRiskIjNGKysrT3BQNGVGLjckJCIkXyJGKyQiMykpKioqKio0LWBePiVGLjckJCIkYCJGKyQiMysrKysielgnM2VGLjckJCIkYSJGKyQiM0ErKytaZ3cmPiVGLjckJCIkYiJGKyQiM18qKioqKjRJTXohZUYuNyQkIiRjIkYrJCIzQysrKyNwaWo+JUYuNyQkIiRkIkYrJCIzRisrKyQqPUMyZUYuNyQkIiRlIkYrJCIzKysrK0ZNJXA+JUYuNyQkIiRmIkYrJCIzWCoqKioqKjQhb2whZUYuNyQkIiRnIkYrJCIzKSkqKioqKj5yM3Y+JUYuNyQkIiRoIkYrJCIzSCsrKzlAImYhZUYuNyQkIiRpIkYrJCIzIyoqKioqKjQqKmUhKT4lRi43JCQiJGoiRiskIjM6KysrIXB0XyFlRi43JCQiJGsiRiskIjM1KysrJHAlZik+JUYuNyQkIiRsIkYrJCIzWisrK0tBbC9lRi43JCQiJG0iRiskIjM/KysrSmk2Kj4lRi43JCQiJG4iRiskIjMnKioqKioqXERaUyFlRi43JCQiJG8iRiskIjMlKioqKioqNCxDJyo+JUYuNyQkIiRwIkYrJCIzaioqKioqekdlTSFlRi43JCQiJHEiRiskIjM7KysrPiU9LD8lRi43JCQiJHIiRiskIjNtKioqKipwKFspRyFlRi43JCQiJHMiRiskIjN5KioqKio0JCkqZitVRi43JCQiJHQiRiskIjNYKysrJWVFQiFlRi43JCQiJHUiRiskIjN3KioqKipSaG81PyVGLjckJCIkdiJGKyQiM1YrKyt1SHksZUYuNyQkIiR3IkYrJCIzJioqKioqKj43RDo/JUYuNyQkIiR4IkYrJCIzYSoqKioqUmtgNyFlRi43JCQiJHkiRiskIjM8KysrJnBwPj8lRi43JCQiJHoiRiskIjNkKioqKipmPVEyIWVGLjckJCIkIT1GKyQiM0IrKyttRVMtVUYuNyQkIiQiPUYrJCIzJCkqKioqKmZATy0hZUYuNyQkIiQjPUYrJCIzPysrK2VWI0c/JUYuNyQkIiQkPUYrJCIzSSsrK2V0dSp6JkYuNyQkIiQlPUYrJCIzdCoqKioqPjNOSz8lRi43JCQiJCY9RiskIjNrKioqKipwQ3IjKnomRi43JCQiJCc9RiskIjN4KioqKip6OE5PPyVGLjckJCIkKD1GKyQiMykpKioqKio+YDIpKXomRi43JCQiJCk9RiskIjMkKSoqKioqPiNbLS9VRi43JCQiJCo9RiskIjM7KysrbmVOKXomRi43JCQiJCE+RiskIjM0KysrPVdTL1VGLjckJCIkIj5GKyQiMysrKytDZiJ6eiZGLjckJCIkIz5GKyQiMyoqKioqKio0P3VaPyVGLjckJCIkJD5GKyQiM1MrKyshUShbKHomRi43JCQiJCU+RiskIjMiKSoqKioqZlZNXj8lRi43JCQiJCY+RiskIjNPKysrRipwcXomRi43JCQiJCc+RiskIjMwKysrJVEmWzBVRi43JCQiJCg+RiskIjNvKioqKipcRGpteiZGLjckJCIkKT5GKyQiMzsrKysrdCNlPyVGLjckJCIkKj5GKyQiMzcrKyt3cUUneiZGLjckJCIkKyNGKyQiMyUpKioqKio+VWdoPyVGLi0lJ0xFR0VOREc2I1EoQ3VydmV+MTYiLSUmQ09MT1JHNiYlJFJHQkckRmhuISIiJEYrRltbb0ZcW28tJSZUSVRMRUc2JFE3MS1sb2N1c34yLWFsbGVsZX5tb2RlbEZlam4tSSVGT05UR0Zlam42JUkmVElNRVNHRmVqbkklQk9MREdGZWpuRmJvLSUrQVhFU0xBQkVMU0c2JFEzZ2VuZXJhdGlvbn5udW1iZXJ+RmVqblEieEZlam4tJSZTVFlMRUc2IyUmUE9JTlRHRmFbby0lJVZJRVdHNiQ7JCEjU0ZbW28kIiQvI0YrOyQiMisrU0owRiNmUiEjPCQiMS0rJykqRz8lemchIzs=

g(x);

NiMqJCwmIiIiRiUqJkkoZXBzaWxvbkc2IkYlSSJ4R0YoIiIjRiUhIiI=

factor(g(g(x))-x);

NiMsJCooLCgqJkkoZXBzaWxvbkc2IiIiIkkieEdGKCIiI0YpKiZGKkYpRidGKSEiIkYpRilGKSwoKiZGJ0YpRioiIiRGKUYqRilGLUYpRiksKkYpRilGJkYrKiZGJ0YrRioiIiVGKUYnRilGLUYt

fsolve(subs(e=4,e*x^3+x-1),x);

NiMkIisrKysrXSEjNQ==

Frequency-dependent selection in one-locus two-allele diploid model (Gavrilets and Hastings, 1995, Proc.R.Soc.Lond.B, 261: 233-238). i) Derivation of the dynamics equation; ii) Analyses of conditions for existance and stability of equilibria; iii) Numerical solutions - intermittency and transient chaos restart:with(linalg): loading librariesF:=array(1..3,1..3,[[1,beta,alpha],[g,eta,g],[alpha,beta,1]]); matrix of coefficientswAA:=F[1,1]*p^2+F[1,2]*2*p*q+F[1,3]*q^2; fitness of AAwAa:=F[2,1]*p^2+F[2,2]*2*p*q+F[2,3]*q^2; fitness of Aawaa:=F[3,1]*p^2+F[3,2]*2*p*q+F[3,3]*q^2; fitness of aawmean:=wAA*p^2+wAa*2*p*q+waa*q^2; mean fitness of the populationwA:=wAA*p+wAa*q: wa:=wAa*p+waa*q: induced allele fitnessq:=1-p: f:=p*(wA-wmean)/wmean; allele frequency in the next generationfactor(%);numer(%)/(p*(2*p-1)*(p-1));factor(%+1-g);f:= p+var(p)*(p-q)*(1-g-Omega*var(p))/(1-2*(2-beta-g)*var(p)+2*Omega*(var(p))^2);allele frequency in the next generationS:=diff(f,p): eigenvaluesvar(p):=p*q;subs(p=0,S); eigenvalues at p=0subs(p=1,S); eigenvalues at p=1subs(p=1/2,S); eigenvalues at p=1/2factor(%);var(p):='var(p)':S;subs({diff(var(p),p)=q-p,var(p)=(1-g)/Omega},%);f:=x->x+x*(1-x)*(2*x-1)*(1-g-Omega*x*(1-x))/(1-2*(2-beta-g)*x*(1-x)+2*Omega*x^2*(1-x)^2);Omega:=C-4*g+(1-beta)^2: new parameter CT:=200: X:=array[0..T]: parameters and initial conditionsX[0]:= .551: C:=.1; g:=1.1; beta:=-2;for i from 1 to T do Y:=f(X[i-1]): X[i]:=Y od: plot( [seq([i,X[i]],i=0..T)], 0..T,0..1,style=point, labels=[ `generation number `, `p` ], title=`Frequency-dependent selection`, titlefont=[TIMES,BOLD,12],font=[TIMES,BOLD,12]);