%This script will compute the statistical mechanical probabilities of a
%polymerase binding to a promoter which can be blocked by a repressor. To
%begin, we'll have to define our parameters.
P = 3E3; %Number of polymerases in the cell.
R = logspace(-1,4, 100); %Range of the number of repressors in the cell.
Nns = 5E6; %Nonspecific binding sites of polymerase (~ size of genome).
dEp = -8; %Energy of polymerase binding to the promoter (units of kT).
dEr = -16; %Energy of the repressor binding the operator (units of kT).
%To simply our computations, we can define the denominator (the partition
%function) separately from the probabilities.
Z = 1 + (P / Nns) * exp(-dEp) + (R ./ Nns)*exp(-dEr);
%Now we'll compute the probabilities of each promoter state.
Pempty = 1 ./ Z; %Nothing bound (no transcription)
Prepressor = (R ./ Nns) * exp(-dEr) ./ Z; %Repressor bound (no transcription)
Pbound = (P / Nns) * exp(-dEp) ./ Z; %Polymerase bound (transcription)
sumP = Pempty + Prepressor + Pbound; %Sum of all probabilites (should be =1)
%Now we should plot our probabilities on a semilog scale as we care about
%a range of different repressor numbers. To do this, we can just call
%the semilogx function.
semilogx(R, Pempty, 'r-') %As a red line.
hold on %This allows many plots to be made on the same axes.
semilogx(R, Prepressor, 'b-') %As a blue line.
semilogx(R, Pbound, 'k-') %As a black line.
semilogx(R, sumP, 'm-') %As a magenta line (not green because we also used red).
ylim([0 1.2]) %Change the limits so we can see the summed term.
hold off %Stop plotting on the same axes.
%Now we just have to apply a legend and axis labels. Northwest tells the
%the legend to be put in the upper lefthand corner of the plot.
legend('Pempty', 'Prepressor', 'Pbound', 'Location','northwest');
xlabel('Number of repressors per cell'); %Apply a x label.
ylabel('Probability of state'); %Apply a y label.