clear all %Don't worry about this command.
% This script simulates a Moran process to describe the genetic drift in a
% population of fixed size. In this model, the total number of individuals
% is kept constant. In this simulation, the population consists of two
% populations, m and n. At each time step along a trajectory, an individual
% is chosen from each population and one is randomly selected to die while
% the other lives and reproduces. This keeps the population fixed in size,
% although the proportion of m or n in the population might change.
%To show the power of genetic drift, we'll run the simulation a number of
%times. First, let's start by defnining some parameters.
timeSteps = 1200; %Number of times we'll randomly pick individuals.
nSimulations = 100; %Number of times we'll run the simulation.
timeOfExtinction = zeros(nSimulations, 1); %Where we'll put the time of ext.
timeOfFixation = zeros(nSimulations, 1); %Where we'll put the time of fix.
%Since we are also making a plot, we'll instantiate a figure.
figure(1)
%Now let's do the simulation!
for x=1:nSimulations
%We can randomly choose which individual lives or dies by generating a
%random number between 0 and 1. If this number is greater than 0.500
% an individual from population n will die while m will live and
% reproduce. If the number is less than 0.500, the opposite will occur.
%Let's generate a list of these random numbers. This is much faster
%than generating a random number at each step of the simulation.
randomSelection = rand(1, timeSteps); %1D matrix.
time = linspace(1, timeSteps, timeSteps);
%We could also include mutation in this model as well. This is
%temporarily commented out.
%mutationProb = rand(1,1000);
%As mentioned above, we can can choose a number in the interval [0,1]
%below which a member from population n will reproduce.
%Change these to see the effect of selection.
probN = 0.5; %Probability N will reproduce.
probM = 1.0 - probN; %Probability N will reproduce.
%Now we need to set the number of each member in the population.
sizeN = 15; %Genetic drift is very important at small pop. sizes.
sizeM = 15;
%Now we can make a list so we can keep track of how many individuals of
%each genotype are present at each time step.
N_t = zeros(timeSteps, 1);
M_t = zeros(timeSteps, 1);
%Now we can interate throughout our random number sequence and decide
%at each point which individual dies and which reproduces.
for i=1:length(randomSelection)
%At everytime point the total population must be kept constant. If
%the number of n or m individuals reaches zero, it will remain at
%zero unless a mutation occurs turning a member from the other
%population into an m or n.
if sizeN > 0 & sizeN < (sizeN + sizeM)
%If these conditions are met, the simulation will proceed with
%selecting reproducers.
if randomSelection(i) < probN
sizeN = sizeN + 1; %Increase n population by 1
sizeM = sizeM - 1; %Decrease m population by 1
N_t(i) = sizeN; %Add current # of n individuals to N_t.
M_t(i) = sizeM; %Add current # of m individuals to M_t
else
sizeN = sizeN - 1;
sizeM = sizeM + 1;
N_t(i) = sizeN;
M_T(i) = sizeM;
end
%If the number of n individuals is 0 or 100, they will stay at this
%number for the rest of the simulation.
else
N_t(i) = sizeN;
M_t(i) = sizeM;
%We should keep track of when a population is fixed or has gone
%extinct.
if sizeN == (sizeM + sizeN) & (N_t(i-1) - N_t(i-2))~=0
timeOfFixation(x) = i-1;
elseif sizeN == 0 & (N_t(i-1) - N_t(i-2)) ~= 0
timeOfExtinction(i) = i-1;
end
end
end
%Now we can plot the trace of the number of n individuals as a function
%of time.
plot(time, N_t, '-', 'linewidth', 1)
hold on
end
%Add labels to our plot so it's clear what we are looking at.
ylabel('Number of "n" individuals')
xlabel('Time (arbitrary units)')
hold off
%Now let's plot a histogram of when the N genotypes reached fixation. Zeros
%in this array represent events where the N individuals died or when
%neither population died. We should first drop those zeros leaving an array
%with only times where fixation was reached.
timeOfFixation(timeOfFixation==0) = [];
%Now we'll make the histogram.
figure(2)
hist(timeOfFixation)
%Let's label the histogram.
xlabel('Time of Fixation')
ylabel('Number of Fixation Events')