%This program performs fertilisation of gametes function [out_pop,Nv,Pself,Ntot] = DiploidFertilisation(inp_pop,NV) %Declare global variables global M global loci global self %Preallocate arrays mate = cell(1,M+1); out_pop = cell(1,M+1); Pself = zeros(1,M+1); %Run function for k = 1:M+1 %Shuffle gametes mate{1,k} = inp_pop{1,k}(:,randperm(NV(k))); %Remove odd gamete H = floor(NV(k)/2); %Ellimination of selfers E = zeros(1,2*H); %Fraction of selfers P = zeros(1,2*H); %Identify selfers for l = 1:H E(2*l) = rand*(mate{1,k}(end,2*l) == mate{1,k}(end,2*l-1))>self; E(2*l-1) = E(2*l); P(2*l) = (mate{1,k}(end,2*l) == mate{1,k}(end,2*l-1))*(1-E(2*l)); P(2*l-1) = P(2*l); end %Accept mate-pairs that fulfil the selfing/no selfing criterion mate{1,k} = mate{1,k}(1:end-1,E==0); %Set new population size after elimination of selfing pairs H = sum(1-E)/2; %Preallocate array for the next diploid generation population out_pop{1,k} = zeros(2,loci,H); %Create the next generation of local diploid population for l = 1:H out_pop{1,k}(1,:,l) = mate{1,k}(:,2*l-1); out_pop{1,k}(2,:,l) = mate{1,k}(:,2*l); end %Calculate the fraction of selfing Pself(k) = max(0,sum(P)/(2*H)); %Set next local population size NV(k) = H; end %Next generation local population sizes Nv = NV; %Next generation total population size Ntot = sum(Nv); end