%Main function for the range expansion simulations function [lattice,N1,N2,Ntot,Vg,P,clinehist,selfpenalty,z,nohap,meanfitness] = DiploidMain(mu,recrate,K,N,theta,inplattice,TL) global M global loci global alleffect n = length(TL); if n > 5 %Use first row for shorter simulations and second row for longer (>7 %days) simulations. %LD_snaps = [1,max(2,round(n/80)),max(3,round(n/20)),max(4,round(n/5)),max(5,round(n/2)),n]; LD_snaps = [1,max(2,round(n/40)),max(3,round(n/10)),max(4,round(n/2.5)),n]; else LD_snaps = 1:n; end LD_count = 1; cline = zeros(loci,M+1); clinehist = zeros(loci,M+1,length(LD_snaps)); lattice = cell(length(LD_snaps),M+1); haplotypes = cell(1,M+1); meanfitness = zeros(n,M+1); nohap = zeros(n,M+1); z = zeros(n,M+1); Vg = zeros(n,M+1); N1 = zeros(n,M+1); N2 = zeros(n,M+1); P = zeros(n,M+1); selfpenalty = zeros(n,M+1); for l = 1:n for k = 1:TL(l) %Step 1: Selection and creation of gametes [gametes,Ngam,r] = DiploidSelection(inplattice,N,theta,K,recrate); %Step 2: Mutations gametes = DiploidMutation(gametes,mu,Ngam); %Step 3: Fertilisation [inplattice,Np,Pself,Ntot] = DiploidFertilisation(gametes,Ngam); %Step 4: Dispersal [inplattice,N] = DiploidMigration(inplattice,Np,Ntot); end %Local population size N1(l,:) = N; N2(l,:) = Np; for in = 1:M+1 %Penalty for disallowing selfing (0 if selfing allowed at no cost) selfpenalty(l,in) = (Np(in) - floor(Ngam(in)/2))/Np(in); %Clines in allele frequencies cline(:,in) = sum(sum(inplattice{1,in},1),3)/(2*alleffect*N(in)); %Each haplotype haplotypes{1,in} = cat(2,reshape(inplattice{1,in}(1,:,:),size(inplattice{1,in},2),size(inplattice{1,in},3)),reshape(inplattice{1,in}(2,:,:),size(inplattice{1,in},2),size(inplattice{1,in},3))); %Unique haplotypes nohap(l,in) = size(unique(haplotypes{1,in}','rows'),1); %The average growth rate in local population meanfitness(l,in) = mean(2*exp(squeeze(r{1,in}))); end %Local average phenotype z(l,:) = mean(cline); %Local genetic variance (LE component) Vg(l,:) = 2*alleffect^2*sum(cline.*(1-cline)); %Fraction of selfing (0 if selfing completely disallowed) P(l,:) = Pself; if l == LD_snaps(LD_count) %All genes clinehist(:,:,LD_count) = cline; for in = 1:M+1 lattice{LD_count,in} = logical(inplattice{1,in}); end LD_count = LD_count + 1; end end end