(* ::Package:: *)
% script to simulate genetic drift resulting from random mating between
% members of a population containing two alleles, "1" and "2"
pop_size = 1000; % number of individuals in population
n_gen = 250; % number of generations to simulate
n_copies = 2; % each individual contains two copies of the gene
% set the probabilities of the alleles in the first generation
p_ 1 = 0.5; % probability of allele "1" in first generation
p_ 2 = 1 - p_ 1; % probability of allele "2" in first generation
pop_mat = zeros(n_gen,pop_size,n_copies);
% initialize the first generation of the simulation
for j=1:pop_size
for l=1:n_copies
r = rand();
if r < p_ 1
pop_mat(1,j,l) = 1;
else
pop_mat(1,j,l) = 2;
end
end
end
% now loop over the rest of the generations of the simulation
for i=2:n_gen
for j=1:pop_size
% randomly choose two individuals from the population to reproduce
J_ 1 = randi(pop_size,1); % first random individual
J_ 2 = randi(pop_size,1); % second random individual
% now choose which allele gets passed on from the first parent
% each allele is chosen with 50% probability
r = rand();
if r < 0.5
pop_mat(i,j,1) = pop_mat(i-1,J_ 1,1);
else
pop_mat(i,j,1) = pop_mat(i-1,J_ 1,2);
end
% now choose which allele gets passed on from the second parent
% again, each allele is chosen with 50% probability
r = rand();
if r < 0.5
pop_mat(i,j,2) = pop_mat(i-1,J_ 2,1);
else
pop_mat(i,j,2) = pop_mat(i-1,J_ 2,2);
end
end
end
% now the simulation is over! all that remains is to plot the results
% we'll use a function called imagesc to make a visual map of how the
% population changes over time.
subplot(2,1,1)
imagesc(sum (pop_mat,3)');
xlabel('Generation number')
ylabel('Individual in population')
% next we will plot how the different allele frequencies change over time
subplot(2,1,2)
% plot probability of allele 1
plot(sum (sum(pop_mat==1,3),2)/(pop_size*2),'b')
hold on
% plot probability of allele 2
plot(sum (sum(pop_mat==2,3),2)/(pop_size*2),'r')
xlabel('Generation number')
ylabel('Allele frequency')
ylim([0,1])
legend('Allele 1', 'Allele 2')
hold off
set(gcf,'units','normalized','outerposition',[0 0 1 1])