%In this script, we'll compute the probability of ligand being bound in a
%cooperative system. In running this script, try to change the energy of
%interaction (dEi) to see how we can produce the Hill function from a
%statistical mechanical view point. We'll begin by declaring some
%variables.
dEb = -5; %Energy of ligand binding to receptor (in kT)
dEi = -2; %Energy of interaction for two ligands bound to receptor (in kT).
L = logspace(-4, -1, 1000); %The range of ligand concentration (L / Omega).
%To simply our computation, we'll define the partition function (Z) first.
Z = 1 + 2 .* L .* exp(-dEb) + L.^2 * exp(-(2*dEb + dEi)); %Partition func.
%Now, let's enumerate the probabilities.
Pempty = 1 ./ Z; %No ligand bound.
Pone = L * exp(-dEb) ./ Z; %One ligand being bound
Pboth = L.^2 * exp(-(2 * dEb + dEi)) ./ Z; %Two ligands bound.
%Now we can plot the probabilities as a function of ligand concentration.
figure(1)
semilogx(L, Pempty, 'r-');
hold on %So we can plot several lines on one set of axes.
semilogx(L, Pone, 'b-');
semilogx(L, Pboth, 'k-');
xlabel('Ligand concentration (L/Omega)');
ylabel('Probability of state');
legend('Pempty', 'Pone', 'Pboth'); %MATLAB will automatically place the legend.
hold off %Stop plotting on the same set of axes.
%We should also compute the average number of bound ligands at each
%concentration. What should it converge to?
Nbound = 2 * Pone + 2 * Pboth;
%Plot the average number of bound ligands.
figure(2)
plot(L, Nbound, 'k-')
xlabel('L / Omega');
ylabel('Average number bound');