% This script simulates Simpson's paraox in a mixed population of bacteria.
% Some of these bacteria (producers) generate a small molecule which induces
% antibiotic resistance. This molecule (autoinducer, AI) is released into these
% environment where another population of bacteria (nonproducers) can use it as
% well. These nonproducers do not synthesize the molecule and thereby suffer no
% cost from making it, yet reap all the benefit. Under specific conditions,
% Simpson's paradox is achieved in which the number of producers decreases
% locally, but increases globally. This script simulates these conditions.
clear all % Don't worry about this.
%Start by setting up the series of equations and parameters.
n0 = 4e7; %Initial number of cells in the population
K = 4e10; %Carrying capacity of the environment.
Km = 3; %[AI] where growth is half maximal (in uM)
kappa = 1.05; %Growth advantage to be a nonproducer.
alpha = 3e-13; %Rate of AI synthesis by producers (in umol/ (cell * min)).
sMax = 0.025; %Maximum growth rate (inverse minutes).
sMin = sMax / 10.0; %Minimum growth rate (inverse minutes).
pg = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 1.0]; %Initial proportion of
%producers in population.
time = 780; %Length of the simulation in minutes.
%Set up empty arrays to which we'll add our results.
pFinal = zeros(length(pg), 1); %Final # producers for each pg.
pInitial = zeros(length(pg), 1); %Initial # producers for each pg.
totalCells = zeros(length(pg), 1); %Total number of cells for each pg.
%Perform the simulation.
for i = 1:length(pg)
%Set up the initial conditions (i.e. time 0).
P_0 = pg(i) * n0; %Initial number of producers.
NP_0 = (1 - pg(i)) * n0; %Initial number of Nonproducers.
AI_0 = 30.0 .* (pg(i)/100.0); %Initial concentration of AI (diluted 100x).
s_0 = sMin + (sMax - sMin) * (AI_0 / (AI_0 + Km)); %Growth rate
%Compute the growth rate, [Producers], [Nonproducers], and [AI] at
% time = 1 min.
AI(1) = AI_0 + alpha * P_0;
P(1) = P_0 + s_0 * P_0 * (1 - ((P_0 + NP_0)/K));
NP(1) = NP_0 + s_0 * kappa * NP_0 * (1 - ((P_0 + NP_0)/K));
s(1) = sMin + (sMax - sMin) * (AI(1) / (AI(1) + Km));
%Perform the simulation for the rest of the time points.
for t = 2:time
AI(t) = AI(t-1) + alpha * P(t-1);
s(t) = sMin + (sMax - sMin) * (AI(t)/ (AI(t) + Km));
P(t) = P(t-1) + s(t-1) * P(t-1) * (1 - (P(t-1) + NP(t-1))/K);
NP(t) = NP(t-1) +...
s(t-1) * kappa * NP(t-1) * (1 - (P(t-1) + NP(t-1))/K);
end
%Extract the initial number of producers, final number of producers,
%and total number of cells (nonproducers + producers).
pInitial(i) = P_0;
pFinal(i) = P(time);
totalCells(i) = NP(time) + P(time);
end
%Compute the normalized growth factor, wg.
wg = pFinal ./ pInitial; %wg for each population.
wg_norm = wg ./ max(wg); %wg normalized to the case where pg=1.0
pFinalProportion = pFinal ./ totalCells;
%Make the plot.
plot(pg, wg_norm, 'r.', 'MarkerSize', 12)
hold on
plot(pFinalProportion, wg_norm, 'k.', 'MarkerSize', 12)
xlabel('Proportion of Producers')
ylabel('normalized wg')
axis([0.0 1.0 0.0 1.0])
legend('Inital Proportion', 'Final Proportion', 'Location', 'northwest')
hold off
%Compute the global change in producer population.
pOrig = sum(pg) ./ length(pg);
pPrime = sum(wg .* pFinalProportion) ./ sum(wg);
deltaP = pPrime - pOrig;
disp(['Global change in Producers: ' num2str(deltaP)])