function n_to_fix = num_gens_to_fixation(pop_size)
% simulate random mating between a population containing two alleles "1"
% and "2" of some gene
n_to_fix = 1; % number of generations it takes to fixation
n_alleles = 2; % number of loci at which the gene appears
p_1 = 0.5; % probability of allele "1" in the first generation
pop_curr = zeros(pop_size,n_alleles);
pop_next = zeros(pop_size,n_alleles);
% initialize the first generation of the population
for j=1:pop_size
for l=1:n_alleles
r = rand();
% randomly choose allele "1" or "2"
if r < p_1
pop_curr(j,l) = 1;
else
pop_curr(j,l) = 2;
end
end
end
% compute the number of allele "1"s in the current generation
n_1_curr = sum(sum(pop_curr==1));
% loop over subsequent generations of the simulation
while ~(n_1_curr==2*pop_size || n_1_curr==0)
% loop over the individuals in the population
for j=1:pop_size
% randomly choose two individuals to reproduce
J_1 = randi(pop_size,1); % first randomly chosen to reproduce
J_2 = randi(pop_size,1); % second randomly chosen to reproduce
% now choose one allele with equal probability from the first
% parent
r = rand();
if r < 0.5
pop_next(j,1) = pop_curr(J_1,1);
else
pop_next(j,1) = pop_curr(J_1,2);
end
% now choose one allele with equal probability from the second
% parent
r = rand();
if r < 0.5
pop_next(j,2) = pop_curr(J_2,1);
else
pop_next(j,2) = pop_curr(J_2,2);
end
end
pop_curr = pop_next;
n_1_curr = sum(sum(pop_curr==1));
n_to_fix = n_to_fix + 1;
end