%Simulation by Marina Rafajlovic %This code is used for simulating the two-locus establishment model %presented in detail in Rafajlovic et al. (2016), Evolution. %It is assumed that a mutation of a fixed mutation-effect size %(denoted by epsilon below) lands at the undifferentiated locus (that has %allele-effect size denoted by Ys=0 below). %(Note that according to the notations used in the paper mentioned, this %corresponds to setting Yw=0.) %The second locus is assumed to be differentiated with alleles of effect %sizes Yi1 and -Yi1 (in notations used below) that are assumed to be in %migration-selection balance at the start of the simulation. %(Note that Yi1 is denoted by Ys in the paper mentioned.) RandStream.setDefaultStream (RandStream('mt19937ar','seed',sum(100*clock))); %parameters set by the user N=1000;%number of individuals in each population sigma=4;%selection parameter r=0.0005;%recombination rate between the two loci repeat=5*10^5;%number of independent realisations m=1/10;%migration rate epsilon=0.05;%mutation-effect size Ys=0;%allele-effect size at the undifferentiated locus Yi1=(0.1:0.2:0.9);%allele-effect size at the differentiated locus. %Varying Yi1 one varies the level of differentiation at this locus. alphai=Yi1-Ys; Yi=Yi1-alphai; theta1=2;%optimal phenotype in the first population theta2=-2;%optimal phenotype in the second populations %quantities of interest trep=zeros(repeat,length(Yi)); %trep saves times needed for a successful establishment of a mutation; if a %muation fails to establish successfully, the time is saved as 0 (this can %be changed to NaN alternatively) prep=zeros(1,length(Yi)); %prep saves the probability of successful establishment for each value of Yi1 pboth=zeros(1,length(Yi)); %pboth saves the probability that both loci have a polymorphism maintained %by the time maxTime (see below). boxesM=[0 m 1+m]; timeSteady=zeros(1,length(Yi)); %timeSteady saves the time needed to reach the intial migration-selection balance %under a deterministic approximation (prior to a mutation) maxTime=10^5; %maxTime sets the maximum time to trace the fate of the mutant alelle, and %to obtain the intial frequencies under a migration-selection balance %the following three matrices save the steady-state frequencies of the 8 %possible genotypes with 2 alelles at each locus freqH0=zeros(length(Yi),8); freqH0_old=10*ones(length(Yi),8); freqH1=zeros(length(Yi)*repeat,8); for I=1:length(Yi) I alpha=alphai(I); Y=Yi(I); p0=0; p11=0; count=0; steady=0; pop1=[1/3 1/3 0 0 1/3 0 0 0 0 0]; pop2=pop1; gametes1=zeros(1,4); %determining the intial condition i=0; while steady==0 i=i+1; %migration pop_new1=(1-m)*pop1+m*pop2; pop_new2=(1-m)*pop2+m*pop1; %subpopulation 1 fitness1=[4*Y+2*alpha, 2*Y, 2*Y+2*alpha, 0, -2*alpha, 0, -2*Y-2*alpha, 2*alpha, -2*Y, -4*Y-2*alpha]; fitness1=exp(-(theta1-fitness1).^2/(2*sigma^2)); fitness1=fitness1/(sum(pop_new1.*fitness1)); gametes1(1)=pop_new1(1)*fitness1(1)+pop_new1(2)*fitness1(2)/2+pop_new1(3)*fitness1(3)/2+pop_new1(4)*fitness1(4)*(1-r)/2+pop_new1(6)*fitness1(6)*r/2; gametes1(2)=pop_new1(2)*fitness1(2)/2+pop_new1(4)*fitness1(4)*r/2+pop_new1(5)*fitness1(5)+pop_new1(6)*fitness1(6)*(1-r)/2+pop_new1(7)*fitness1(7)/2; gametes1(3)=pop_new1(3)*fitness1(3)/2+pop_new1(4)*fitness1(4)*r/2+pop_new1(6)*fitness1(6)*(1-r)/2+pop_new1(8)*fitness1(8)+pop_new1(9)*fitness1(9)/2; gametes1(4)=pop_new1(4)*fitness1(4)*(1-r)/2+pop_new1(6)*fitness1(6)*r/2+pop_new1(7)*fitness1(7)/2+pop_new1(9)*fitness1(9)/2+pop_new1(10)*fitness1(10); freqH0(I,1)=gametes1(1); freqH0(I,2)=gametes1(2); freqH0(I,3)=gametes1(3); freqH0(I,4)=gametes1(4); %make a new population out of these: pop1(1)=gametes1(1)^2; pop1(2)=2*gametes1(1)*gametes1(2); pop1(3)=2*gametes1(1)*gametes1(3); pop1(4)=2*gametes1(1)*gametes1(4); pop1(5)=gametes1(2)^2; pop1(6)=2*gametes1(2)*gametes1(3); pop1(7)=2*gametes1(2)*gametes1(4); pop1(8)=gametes1(3)*gametes1(3); pop1(9)=2*gametes1(3)*gametes1(4); pop1(10)=gametes1(4)^2; fitness1=[4*Y+2*alpha, 2*Y, 2*Y+2*alpha, 0, -2*alpha, 0, -2*Y-2*alpha, 2*alpha, -2*Y, -4*Y-2*alpha]; fitness1=exp(-(theta2-fitness1).^2/(2*sigma^2)); fitness1=fitness1/(sum(pop_new2.*fitness1)); gametes1(1)=pop_new2(1)*fitness1(1)+pop_new2(2)*fitness1(2)/2+pop_new2(3)*fitness1(3)/2+pop_new2(4)*fitness1(4)*(1-r)/2+pop_new2(6)*fitness1(6)*r/2; gametes1(2)=pop_new2(2)*fitness1(2)/2+pop_new2(4)*fitness1(4)*r/2+pop_new2(5)*fitness1(5)+pop_new2(6)*fitness1(6)*(1-r)/2+pop_new2(7)*fitness1(7)/2; gametes1(3)=pop_new2(3)*fitness1(3)/2+pop_new2(4)*fitness1(4)*r/2+pop_new2(6)*fitness1(6)*(1-r)/2+pop_new2(8)*fitness1(8)+pop_new2(9)*fitness1(9)/2; gametes1(4)=pop_new2(4)*fitness1(4)*(1-r)/2+pop_new2(6)*fitness1(6)*r/2+pop_new2(7)*fitness1(7)/2+pop_new2(9)*fitness1(9)/2+pop_new2(10)*fitness1(10); pop2(1)=gametes1(1)^2; pop2(2)=2*gametes1(1)*gametes1(2); pop2(3)=2*gametes1(1)*gametes1(3); pop2(4)=2*gametes1(1)*gametes1(4); pop2(5)=gametes1(2)^2; pop2(6)=2*gametes1(2)*gametes1(3); pop2(7)=2*gametes1(2)*gametes1(4); pop2(8)=gametes1(3)*gametes1(3); pop2(9)=2*gametes1(3)*gametes1(4); pop2(10)=gametes1(4)^2; freqH0(I,5)=gametes1(1); freqH0(I,6)=gametes1(2); freqH0(I,7)=gametes1(3); freqH0(I,8)=gametes1(4); if i>1 delta=sum((freqH0(I,:)-freqH0_old(I,:)).^2); %condition for finding the steady state if delta<=8*10^(-8) count=count+1; if count==10^3 steady=1; timeSteady(1,I)=i; end else count=0; if i>maxTime steady=1; end end end freqH0_old(I,:)=freqH0(I,:); end %these are the starting (steady-state) frequencies of gametes %choose the starting population randomly and introduce a mutation pop1_st=pop1; pop2_st=pop2; for J=1:repeat a=rand(1,N); boxes=cumsum([0 pop1_st]); [~,b]=histc(a,boxes); pop1=[]; c=find(b==1); pop1(c,1:4)=kron([Y Y+alpha Y (Y+alpha)],ones(length(c),1)); c=find(b==2); pop1(c,1:4)=kron([Y Y+alpha Y -(Y+alpha)],ones(length(c),1)); c=find(b==3); pop1(c,1:4)=kron([Y Y+alpha -Y (Y+alpha)],ones(length(c),1)); c=find(b==4); pop1(c,1:4)=kron([Y Y+alpha -Y -(Y+alpha)],ones(length(c),1)); c=find(b==5); pop1(c,1:4)=kron([Y -(Y+alpha) Y -(Y+alpha)],ones(length(c),1)); c=find(b==6); pop1(c,1:4)=kron([Y -(Y+alpha) -Y (Y+alpha)],ones(length(c),1)); c=find(b==7); pop1(c,1:4)=kron([Y -(Y+alpha) -Y -(Y+alpha)],ones(length(c),1)); c=find(b==8); pop1(c,1:4)=kron([-Y Y+alpha -Y (Y+alpha)],ones(length(c),1)); c=find(b==9); pop1(c,1:4)=kron([-Y Y+alpha -Y -(Y+alpha)],ones(length(c),1)); c=find(b==10); pop1(c,1:4)=kron([-Y -(Y+alpha) -Y -(Y+alpha)],ones(length(c),1)); a=rand(1,N); boxes=cumsum([0 pop2_st]); [~,b]=histc(a,boxes); pop2=[]; c=find(b==1); pop2(c,1:4)=kron([Y Y+alpha Y (Y+alpha)],ones(length(c),1)); c=find(b==2); pop2(c,1:4)=kron([Y Y+alpha Y -(Y+alpha)],ones(length(c),1)); c=find(b==3); pop2(c,1:4)=kron([Y Y+alpha -Y (Y+alpha)],ones(length(c),1)); c=find(b==4); pop2(c,1:4)=kron([Y Y+alpha -Y -(Y+alpha)],ones(length(c),1)); c=find(b==5); pop2(c,1:4)=kron([Y -(Y+alpha) Y -(Y+alpha)],ones(length(c),1)); c=find(b==6); pop2(c,1:4)=kron([Y -(Y+alpha) -Y (Y+alpha)],ones(length(c),1)); c=find(b==7); pop2(c,1:4)=kron([Y -(Y+alpha) -Y -(Y+alpha)],ones(length(c),1)); c=find(b==8); pop2(c,1:4)=kron([-Y Y+alpha -Y (Y+alpha)],ones(length(c),1)); c=find(b==9); pop2(c,1:4)=kron([-Y Y+alpha -Y -(Y+alpha)],ones(length(c),1)); c=find(b==10); pop2(c,1:4)=kron([-Y -(Y+alpha) -Y -(Y+alpha)],ones(length(c),1)); t0=0; %now introduce one mutation of size epsilon on the first locus %here, this locus is set to be less diverged than the other one; %find one parent in the first (positive) subpopulation, that %carries allele of size Ys, and substitute this allele by the one %with size Ys+epsilon %trace the frequencies of alleles Ys and Ys+epsilon at the first %locus, until either Ys or Ys+epsilon goes extinct, or until %Ys+epsilon becomes most frequent in the first population, while Ys %is most frequent in the second population (successful %establishement event) a=[pop1(:,1); pop1(:,3)]; b=find(abs(a-Y)<=10^(-4)); c=floor(1+rand(1,1).*(length(b)-1)); b=b(c); if b<=N pop1(b,1)=Y+epsilon; else pop1(b-N,3)=Y+epsilon; end %this is the initial condition; %take the maximum time to follow the dynamics, or until loss or %establishment occurs; counter=0; while counter==0 t0=t0+1; %migration a=rand(2,N); [~,c]=histc(a(1,:),boxesM); c1=find(c==1); pop_new2=pop1(c1,:); pop1(c1,:)=[]; pop_new1=pop1; [~,c]=histc(a(2,:),boxesM); c1=find(c==1); pop_new1=[pop_new1;pop2(c1,:)]; pop2(c1,:)=[]; pop_new2=[pop_new2;pop2]; %subpopulation 1 fitness1=sum(pop_new1'); fitness1=exp(-(theta1-fitness1).^2/(2*sigma^2)); fitness1=fitness1/(sum(fitness1)); a=rand(1,2*N); [~,c]=histc(a,cumsum([0 fitness1])); gametes1=pop_new1(c,:); a=rand(1,2*N); [~,c]=histc(a,[0 r 1+r]); c1=find(c==1); c2=c1; gametes_new1=zeros(2*N,2); if length(c1)>0 a=rand(1,length(c1)); a1=find(a<=0.5); gametes_new1(c1(a1),:)=gametes1(c1(a1),[1 4]); c1(a1)=[]; gametes_new1(c1,:)=gametes1(c1,[3 2]); end c=1:2*N; c(c2)=[]; if length(c)>0 a=rand(1,length(c)); a1=find(a<=0.5); a2=c(a1); gametes_new1(a2,:)=gametes1(a2,[1 2]); c(a1)=[]; gametes_new1(c,:)=gametes1(c,[3 4]); end %subpopulation 2 fitness2=sum(pop_new2'); fitness2=exp(-(theta2-fitness2).^2/(2*sigma^2)); fitness2=fitness2/(sum(fitness2)); a=rand(1,2*N); [~,c]=histc(a,cumsum([0 fitness2])); gametes2=pop_new2(c,:); a=rand(1,2*N); [~,c]=histc(a,[0 r 1+r]); c1=find(c==1); gametes_new2=zeros(2*N,2); c2=c1; if length(c1)>0 a=rand(1,length(c1)); a1=find(a<=0.5); gametes_new2(c1(a1),:)=gametes2(c1(a1),[1 4]); c1(a1)=[]; gametes_new2(c1,:)=gametes2(c1,[3 2]); end c=1:2*N; c(c2)=[]; if length(c)>0 a=rand(1,length(c)); a1=find(a<=0.5); a2=c(a1); gametes_new2(a2,:)=gametes2(a2,[1 2]); c(a1)=[]; gametes_new2(c,:)=gametes2(c,[3 4]); end pop1=[gametes_new1(1:N,:) gametes_new1(N+1:end,:)]; pop2=[gametes_new2(1:N,:) gametes_new2(N+1:end,:)]; a=[pop1(:,1);pop1(:,3);pop2(:,1);pop2(:,3)]; c=length(find(abs(a-Y-epsilon)<10^(-4))); if c==0 counter=1; else a=[pop1(:,1);pop1(:,3)]; b=unique(a); [b1,~]=histc(a,b); b2=find(b1==max(b1)); if length(b2)==1 a1=b(b2);%most frequent allele in the first population else a1=NaN; end a=[pop2(:,1);pop2(:,3)]; b=unique(a); [b1,~]=histc(a,b); b2=find(b1==max(b1)); if length(b2)==1 a2=b(b2);%most frequent alelle in the second population else a2=NaN; end %the extent of divergence at the initially undifferentiated locus %is equal to a1-a2; if this is equal to epsilon, then the %mutation has established successfully; due to numerical %precision, take a convenient condition to check whether a1-a2 %is sufficinelty close to epsilon; since epsilon is about %10^(-2), the condition given next is sufficient if abs(a1-a2-epsilon)<=10^(-4) counter=1; trep(J,I)=t0; p0=p0+1; hapl1=sum(gametes_new1'); a=abs(hapl1-(2*Y+alpha+epsilon)); freqH1((I-1)*repeat+J,1)=length(find(a<=10^(-4)))/(2*N); a=abs(hapl1-(epsilon-alpha)); freqH1((I-1)*repeat+J,2)=length(find(a<=10^(-4)))/(2*N); a=abs(hapl1-(alpha)); freqH1((I-1)*repeat+J,3)=length(find(a<=10^(-4)))/(2*N); a=abs(hapl1-(-2*Y-alpha)); freqH1((I-1)*repeat+J,4)=length(find(a<=10^(-4)))/(2*N); hapl1=sum(gametes_new2'); a=abs(hapl1-(2*Y+alpha+epsilon)); freqH1((I-1)*repeat+J,5)=length(find(a<=10^(-4)))/(2*N); a=abs(hapl1-(-alpha+epsilon)); freqH1((I-1)*repeat+J,6)=length(find(a<=10^(-4)))/(2*N); a=abs(hapl1-(alpha)); freqH1((I-1)*repeat+J,7)=length(find(a<=10^(-4)))/(2*N); a=abs(hapl1-(-2*Y-alpha)); freqH1((I-1)*repeat+J,8)=length(find(a<=10^(-4)))/(2*N); else if t0>=maxTime counter=1; p11=p11+1; end end end end end prep(1,I)=p0/repeat; pboth(1,I)=p11/repeat; end %The quantities of interest (saved for further analyses) are, e.g. %parameters, time to reach a successful establishment (trep), the %probability of successful establishment (prep), etc.