clear all %Don't worry about this.
%This script simulates the Luria-Delbruck fluctuation test. In this
%simulation, we will measure the quantitative and qualitative differences
%between "induced mutation" and "inherited mutation" of a population. These
%two models show different behavior.
%We'll first model induced mutation. In this model, a certain proportion of
%cells will mutate upon exposure to some toxic substance (i.e. antibiotic)
%and will become resistant. This mutation will only be experienced by that
%cell and the trait will not be passed to others in the population.
%Lets first define some parameters for the simulation. For the case of
%induced mutation, we'll have to define the rate at which the cells will
%mutate, the total number of cells in the population, and the number of
%simulations to run.
mutRate = 1e-5; %Probability of a cell mutating.
cells0 = 100; %Initial number of cells.
nDivisions = 10; %Number of cell divisions.
nCells = cells0 * 2^nDivisions; %Final number of cells.
nSimulations = 1000; %Total number of simulations to be run.
%Now we'll run the simulation.
for i=1:nSimulations
%To identify which cells will be mutated, we'll generate a 1D matrix of
%random numbers on the interval (0,1) the same length as our number of
%cells. We'll say that only those random numbers with a value less than
%our mutation rate have actually mutated. For each simulation, we'll
%sum the number of cells that have actually mutated and add it to our
%adaptiveOutcomes array.
randomChance = rand(nCells, 1); %Generate random numbers
adaptiveMutants = sum(randomChance < mutRate); %Find the mutants.
%Now we just need to add it to an array we'll read later.
adaptiveOutcomes(i) = adaptiveMutants;
end
%Let's take a look at the histogram of the values.
figure(1)
%We'll use the normal 'hist' command, but we'll change the binning. We want as
%many bins as the maximum number of mutants and have them spaced by one at a time.
hist(adaptiveOutcomes, 0:1:max(adaptiveOutcomes))
xlabel('Number of mutants');
ylabel('Count');
title('Case of adaptive mutation');
%This is a pretty narrow distribution! At most, we only have 5-6 mutants
%throughout the whole simulation (1000 trials!). Let's look at some of the
%statistics of the distribution.
meanAdaptive = mean(adaptiveOutcomes)
varianceAdaptive = var(adaptiveOutcomes)
fanoAdaptive = varianceAdaptive / meanAdaptive;
%If we run this simulation many times, we'll see that we occasionally get a fano
%factor > 1. The mean and the variance are very close to the same value,
%however, so we can say that the Fano factor is about 1.
%Now let's simulate the case in which mutations are inherited. We have to be a
%little clever with our programming to make this work. We have to take into
%account that the daughters of mutated cells will also be mutated. Because of
%this fact, we will need to iterate throughout all of the divisons, rather than
%just through the final number of cells.
for i=1:nSimulations
%We'll set up our initial matrix of cells. None of them are mutated, so
%we'll assign them as zero.
cells = zeros(cells0, 1);
%Now we'll iterate through all of the cell divisions.
for j=1:nDivisions
%Now we'll perform the cell division.
cellsDup = cells; %Duplicate the cells.
%Now generate random numbers and decide who gets mutated.
mutators = rand(length(cellsDup), 1) < mutRate;
%Now we'll mark anybody in the cellsDup array who mutated with a
%1 and anyone in the cells matrix who was mutated with a 1.
cellsDup(mutators | cells) = 1; % '|' is 'or'.
%Now merge the mutators and the original cells together.
cells = cat(1, cellsDup, cells);
end
%Now add the number of mutated cells to our randomOutcomes matrix.
randomOutcomes(i) = sum(cells);
end
%Let's make another histogram.
figure(2)
hist(randomOutcomes, 0:1:max(randomOutcomes)); %Using custom binning.
xlabel('Number of mutants');
ylabel('Count')
title('Case of random mutation')
%That looks pretty different than the adaptive case! Let's take a look at some
%statistics.
meanRandom = mean(randomOutcomes)
varianceRandom = var(randomOutcomes)
fanoRandom = varianceRandom / meanRandom
%That's a much larger Fano factor! This is certainly not in the Poisson regime
%(where mean = variance). Let's show this by plotting a Poisson distribution
%over the histograms.
figure(1)
hold on
adaptivePoisson = poisspdf(0:1:max(adaptiveOutcomes), nCells * mutRate) * 1000;
plot(0:1:max(adaptiveOutcomes), adaptivePoisson, 'r-'); %Poisson as red line.
legend('Distribution', 'Poisson');
figure(2)
hold on
randomPoisson = poisspdf(0:1:max(randomOutcomes), nCells * mutRate) * 1000;
plot(0:1:max(randomOutcomes), randomPoisson, 'r-'); %Poisson as red line.
legend('Distribution', 'Poisson');
%By zooming in on the random mutation case, we can see that the Poisson is a
%pretty terrible fit of the distribution!
%In class, we developed a theoretical prediction of the relationship between the
%number of mutated cells and the mutation rate. This relationship is given as
%
% R(t) = a*N_t * ln(a*Nt)
%
%where R is the number of resistant cells, a is the mutation rate, and N_t is
%the number of cells at time (t). Let's modify our above procedure to vary the
%mutation rate and observe the effect on the number of resistant cells. Rather
%than only making one iteration, we'll just iterate over a range of mutation rates.
nSimulations = 50; %Simulations for each iteration. For speed, < 1000.
mutRate = logspace(-6, -3);
%We'll start with the case of adaptive mutation.
for i=1:length(mutRate)
for j=1:nSimulations
%For the sake of my fingers, we'll write our above adaptive mutation
%function as a single line.
adaptiveOutcomes2(j) = sum(rand(nCells,1) < mutRate(i));
end
%Now, add the mean of the number of resistant cells to a matrix.
adaptiveOutcomesVariedMut(i) = mean(adaptiveOutcomes2);
end
%Now we can do this same procedure for the random mutation case.
for i=1:length(mutRate)
for j=1:nSimulations
%Make the starting cells.
cells = zeros(cells0, 1);
for k=1:nDivisions
cellsDup = cells; %Duplicate the cells.
%Find the mutators.
mutators = rand(length(cellsDup),1) < mutRate(i);
%Add mutators to the matrix.
cellsDup(mutators | cells) = 1;
%Combine the parents and daughters.
cells = cat(1, cellsDup, cells);
end
%Find the total number of mutants in each simulation.
randomOutcomes2(j) = sum(cells);
end
%Find the mean number of mutated cells at mutRate(i).
randomOutcomesVariedMut(i) = mean(randomOutcomes2);
end
%Now we should plot these simulations along with our theoretical prediction.
%Let's write the theoretical prediction first.
theoMuts = mutRate .* nCells .* log(mutRate .* nCells);
%Now we can plot everything together.
figure(3)
plot(mutRate, adaptiveOutcomesVariedMut, 'ro');
hold on
plot(mutRate, randomOutcomesVariedMut, 'bo');
plot(mutRate, theoMuts, 'k-'); %Theoretical curve.
xlabel('Mutation rate');
ylabel('Number of resistant cells');
legend('Adaptive case', 'Random case', 'Theoretical prediction');
%Unsurprisingly, our theoretical curve agrees quite nicely with the case of
%random mutations, but not at all in the case of adaptive mutation. In fact, our
%theoretical prediction is a slight underestimate of the simulation. This is
%understandable because of reasons described in class.
%How does the mutation rate scale with the number of final cells? We could solve
%the above equation by hand, but the mathetmatics gets complicated. If you
%happen to have the 'solve' matlab package installed, we can solve it
%trivially.
%We'll start by defining the symbols we'll use.
syms R, a, Nt
%Define the equation we'd like to solve.
eqn = R == a*Nt * log(a*Nt); %'==' in this case is NOT a logical operator.
%Now we just tell Matlab to solve it.
sol = solve(eqn, a); %Solve the equation with respect to a.
sol;
%The equation spits out the following result
%
% a = R / (Nt * wrightOmega(-log(1/R))).
%
%The 'wrightOmega' is a specific mathematical function which you can look up on
%your own. Let's see how this varies as we change the number of cells. What is
%the mutation rate needed to get 50 mutants at a range of N_final.
R = 50; %Number of mutants we want.
N_final = logspace(2, 6, 1000);
mutRate = R ./ (N_final .* wrightOmega(-log(1./R)));
%Now we'll plot it.
figure(4)
loglog(mutRate, N_final, 'b-'); %As a blue line.
xlabel('Mutation rate');
ylabel('N final needed for 50 mutants.')