clear all
%This script will numerically integrate the probability of a receptor being
%bound as a function of time. While the solution can be computed
%analytically, this will be a good example of how to do this type of
%numerical approximation. Let's begin by defining some variables.
%Define the variables.
K_on = 0.001; %On rate for ligand to receptor (in s^-1)
K_off = [0.001 0.002 0.004 0.008 0.016]; %The range of K_off we'll use (in s^-1).
valueSaturated = zeros(1, length(K_off)); %Vector to add the value at which we saturated.
%Now we can perform the numerical integration. Since we'll be cycling
%through a range of K_off, we'll need two layers to the for loop. In
%addition, we'll create two figures showing the various trajectories.
figure(1)
for j=1:length(K_off) %For each entry in K_off,
dt = 0.01 / (K_on + min(K_off)); %Time step that' we'll integrate by.
time = 0:dt:500*dt; %Total length of the integration (in s^-1).
Pbound = zeros(1,length(time)); %Vector to which we'll at Pbound.
for i=1:length(Pbound) - 1
%We'll first compute the change in Pbound.
dPbound = (1- Pbound(i))*K_on*dt - Pbound(i)*K_off(j) *dt;
%Now we'll simply add the change in Pbound to the previous Pbound value.
Pbound(i+1) = Pbound(i) + dPbound;
end %Finish the ith round of integration.
%Now we add the saturated value to our empty vector.
valueSaturated(j) = max(Pbound);
%Add a plot of the j'th trajectory versus time.
plot(time, Pbound, '-', 'LineWidth', 3); %As a line of width 3.
hold on; %Don't clear the plot each time we make the plot.
end %End the entire integration.
hold off %Prevent anything else from being plotted.
xlabel('Time (s)');
ylabel('P_bound');
legend('K_off = 0.001', 'K_off = 0.002', 'K_off = 0.004', 'K_off = 0.008',...
'K_off = 0.016', 'Location', 'northwest');
%Let's also take a look at the saturated values, and then plot them as a
%function of K_off.
valueSaturated
figure(2)
plot(K_off, valueSaturated, 'k-o') %As black circles connected by lines.
xlabel('K_off (s^-1)');
ylabel('Saturation Value');