restart:<Text-field firstindent="0.0" layout="Heading 1" leftmargin="0.0" rightmargin="0.0" style="Heading 1"><Font executable="false">Frequency-dependent selection: haploid model</Font></Text-field>plotting solutionseq:=2*s*p(t)*(1-p(t))*(P-p(t)): # dynamic equations:=1: P:=.5: t_max:=10: # parametersinit_cond:={seq([0,i/10],i=1..9)}; # initial conditionswith(DEtools): # load librariesDEplot(diff(p(t),t)=eq,p(t),0..t_max, init_cond, p=0..1,title=`Linear frequency-dependent selection in 1-locus 2-allele haploid population`); # plot<Text-field firstindent="0.0" layout="Heading 1" leftmargin="0.0" rightmargin="0.0" style="Heading 1"><Font executable="false">Coevolution of two haploid species: no within-species interactions</Font></Text-field>With no between-species interactions, the dynamics are described by dp1/dt = p1 q1 b (p2-q2), dp2/dt = p2 q2 c (p1-q1).The solutions of the former satisfy to (p1 q1)^c / (p2 q2)^b = const. The following plots these curves.restart:with(plots): z:= (x*(1-x))^c/(y*(1-y))^b;Let b and c have the same sign (competitive or cooperative interactions)b:=1: c:=1: contourplot(z, x=0.1..0.9,y=0.1..0.9);Let b and c have opposite signs (host-parasite, victim-exploiter type interactions)b:=-1: c:=1: contourplot(z, x=0.1..0.9,y=0.1..0.9);<Text-field firstindent="0.0" layout="Heading 1" leftmargin="0.0" rightmargin="0.0" style="Heading 1"><Font executable="false">Coevolution of two haploid species: general case</Font></Text-field>restart:Coevolutionary dynamics of two haploid species (linear symmetric frequency-dependent selection)eq:=[diff(x(t),t)=x(t)*(1-x(t))*(a*(2*x(t)-1)+b*(2*y(t)-1)), diff(y(t),t)=y(t)*(1-y(t))*(c*(2*x(t)-1)+d*(2*y(t)-1))]: # dynamic equationsa:=-1.0 : b:=1.2: c:=-1.3: d:=1.01: t_max:=500: # parametersinit_cond:={seq([0,1/4,j/4],j=1..3),seq([0,3/4,j/4],j=1..3)}; # initial conditionsinit_cond:={[0,1/4,3/4],[0,.5,.6]}:with(DEtools): # load librariesphaseportrait( eq, [x(t),y(t)], 0..t_max, init_cond,stepsize=.1,x=0..1, y=0..1, arrows=THIN, title=`Coevolutionary dynamics of two haploid species`); # plot<Text-field firstindent="0.0" layout="Heading 1" leftmargin="0.0" rightmargin="0.0" style="Heading 1"><Font executable="false"> Complete stability analysis and dynamics for the mimicry model</Font></Text-field>Dynamic equations f1:=x*(1-x)*(a*(2*x-1)+b*(2*y-1)); f2:=y*(1-y)*(c*(2*x-1)+d*(2*y-1)); Equilibriasolve({f1=0,f2=0},{x,y});factor(resultant(f1,f2,x));with(linalg): Define stability matrixS:=array([[diff(f1,x),diff(f1,y)],[diff(f2,x),diff(f2,y)]]); S1:=subs([x=0,y=0],evalm(S)): lambda1:=eigenvals(evalm(S1));S2:=subs([x=1,y=0],evalm(S)): lambda2:=eigenvals(evalm(S2));S3:=subs([x=0,y=1],evalm(S)): lambda3:=eigenvals( evalm(S3));S4:=subs([x=1,y=1],evalm(S)): lambda4:=eigenvals( evalm(S4));S5:=subs([x=0,y=(d+c)/(2*d)],evalm(S)): lambda5:=eigenvals( evalm(S5));S6:=subs([x=1,y=(d-c)/(2*d)],evalm(S)): lambda6:=eigenvals( evalm(S6));S7:=subs([x=(a+b)/(2*a),y=0],evalm(S)): lambda3:=eigenvals( evalm(S7));S8:=subs([x=(a-b)/(2*a),y=1],evalm(S)): lambda3:=eigenvals( evalm(S8));S9:=subs([x=1/2,y=1/2],evalm(S)): lambda9:=eigenvals( evalm(S9)); determ:=det(S9); tr:=factor(trace(S9));restart:Coevolutionary dynamics of two haploid species (linear symmetric frequency-dependent selection)eq:=[diff(x(t),t)=x(t)*(1-x(t))*(a*(2*x(t)-1)+b*(2*y(t)-1))+mu*(1-2*x(t)), diff(y(t),t)=y(t)*(1-y(t))*(c*(2*x(t)-1)+d*(2*y(t)-1))+mu*(1-2*y(t))]: # dynamic equationsa:=1.1 : b:=1.2: c:=-1.1: d:=1: mu:=.01: t_max:=500: # parametersinit_cond:={seq([0,1/4,j/4],j=1..3),seq([0,3/4,j/4],j=1..3)}; # initial conditionsinit_cond:={[0,1/4,3/4],[0,.5,.6]}:with(DEtools): # load librariesphaseportrait( eq, [x(t),y(t)], 0..t_max, init_cond,stepsize=.1,x=0..1, y=0..1, arrows=THIN, title=`Coevolutionary dynamics of two haploid species`); # plot<Text-field layout="Heading 1" style="Heading 1">Udovic model: diploid popualtion</Text-field>restart:F:=p+p*q*(w[A]-w[a])/wmean;w[A]:=w[AA]*p+w[Aa]*q;w[a]:=w[Aa]*p+w[aa]*q;wmean:=w[A]*p+w[a]*q;w[AA]:=f(p-q); w[Aa]:=h(p-q); w[aa]:=f(q-p);q:=1-p;factor(F-p);subs(p=1/2,%);diff(F,p);subs(p=1/2,%);factor(subs({f(0)=1,h(0)=1-s},%));factor(%);<Text-field firstindent="0.0" layout="Heading 1" leftmargin="0.0" rightmargin="0.0" style="Heading 1"><Font executable="false">Frequency-dependent selection: diploid population</Font></Text-field>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 coeeficientswAA:=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 alllele 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); # eigenvalue at p=1subs(p=1/2,S); # eigenvalue 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 conditions:X[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 `, `frequency` ], title=`Frequency-dependent selection`);