% sexsim.m, Joan Roughgarden, February 22, 2002, touched August 12, 2004
% sexual population displacing asexual population in fluctuating environment
% coded for Matlab, Vers 5, Student Edition

s=2;
u=.0001;
v=1-u;
tmax=30;
ght=1;
rand('seed',0);

xa11=.25;
xa12=.5;
xa22=.25;
xs11=.25;
xs12=.5;
xs22=.25;
na=.5;
ns=.5;

nahis=[na];
nshis=[ns];

f0=figure;
f1=subplot(4,1,1);
f2=subplot(4,1,2);
f3=subplot(4,1,3);
f4=subplot(4,1,4);

for t=1:tmax

 subplot(f1)
 hold on
 axis([0 tmax 0 1])
 axis off
 title(['Sexual vs. Asexual Reproduction, Generation ' num2str(t)])
 text(.5,1,'Abundance')
 plot(0:t-1,nshis(1:t),'y')
 plot(0:t-1,nahis(1:t),'c')
 plot(t-1,nshis(t),'xy')
 plot(t-1,nahis(t),'oc')

 if rand>0.5 w11=s; else w11=1/s; end
 if rand>0.5 w12=s; else w12=1/s; end
 if rand>0.5 w22=s; else w22=1/s; end

 subplot(f2)
 hold off
 cla
 hold on
 axis([0 8 0 1.1*ght])
 axis off
 text(3.15,ght,'Environment')
 if (w11==s) plot(1,ght,'*g'); else plot(1,ght,'*r'); end
 if (w12==s) plot(2,ght,'*g'); else plot(2,ght,'*r'); end
 if (w22==s) plot(3,ght,'*g'); else plot(3,ght,'*r'); end
 if (w11==s) plot(5,ght,'*g'); else plot(5,ght,'*r'); end
 if (w12==s) plot(6,ght,'*g'); else plot(6,ght,'*r'); end
 if (w22==s) plot(7,ght,'*g'); else plot(7,ght,'*r'); end
 text(3.6,ght/2,'Initial')
 text(0,ght/2,'Sexual')
 text(7.1,ght/2,'Asexual')
 bar([1,2,3],[xs11,xs12,xs22],'y')
 bar([5,6,7],[xa11,xa12,xa22],'c')

 wa=xa11*w11+xa12*w12+xa22*w22;
 ws=xs11*w11+xs12*w12+xs22*w22;

 nap=wa*na;
 nsp=ws*ns;
 na=nap/(nap+nsp);
 ns=nsp/(nap+nsp);

 nahis=[nahis,na];
 nshis=[nshis,ns];

 xa11p=w11*xa11/wa;
 xa12p=w12*xa12/wa;
 xa22p=w22*xa22/wa;

 xs11p=w11*xs11/ws;
 xs12p=w12*xs12/ws;
 xs22p=w22*xs22/ws;

 xa11pp=v*v*xa11p+u*v*xa12p+u*u*xa22p;
 xa12pp=2*u*v*xa11p+(u*u+v*v)*xa12p+2*u*v*xa22p;
 xa22pp=u*u*xa11p+u*v*xa12p+v*v*xa22p;

 xs11pp=v*v*xs11p+u*v*xs12p+u*u*xs22p;
 xs12pp=2*u*v*xs11p+(u*u+v*v)*xs12p+2*u*v*xs22p;
 xs22pp=u*u*xs11p+u*v*xs12p+v*v*xs22p;

 subplot(f3)
 hold off
 cla
 hold on
 axis([0 8 0 3*ght/4])
 axis off
 text(3.05,ght/2,'After Selection')
 bar([1,2,3],[xs11pp,xs12pp,xs22pp],'y')
 bar([5,6,7],[xa11pp,xa12pp,xa22pp],'c')

 xa11=xa11pp;
 xa12=xa12pp;
 xa22=xa22pp;

 p=xs11pp+xs12pp/2;
 q=xs22pp+xs12pp/2;

 xs11=p*p;
 xs12=2*p*q;
 xs22=q*q;

 subplot(f4)
 hold off
 cla
 hold on
 axis([0 8 0 3*ght/4])
 axis off
 text(3.25,ght/2,'After Mating')
 bar([1,2,3],[xs11,xs12,xs22],'y')
 bar([5,6,7],[xa11,xa12,xa22],'c')
 text(3.3,-.2*ght,'Genotypes')

 pause

end;