%% Binding of a Ligand to Two Independent Binding Sites of the Same Receptor 

% This example shows how to analyze a binding experiment in which a ligand
% binds to two independent sites on a receptor.

% Copyright 2015 Domenico Gatti.

%%
% We will start by loading a project containing the model of the binding of
% a ligand to 2 different sites (with different Kd's) in the same receptor.
% The best way to simulate two independent binding sites in a single
% receptor is to equate the system to two receptors each with one binding
% site. Then, the total concentration of the receptor is given by the sum
% of the two receptors. This is shown in the following SimBiology project:

addpath(genpath('../TOOLBOXES/BINDING_KINETICS'));
sbioloadproject('../TOOLBOXES/BINDING_KINETICS/Two_Independent_Binding_Sites');
 
% We can also load the model interactively by starting the SimBiology
% desktop and opening the '.sbproj' file, or by importing the data from the
% xml file:

m1 = sbmlimport('../TOOLBOXES/BINDING_KINETICS/Two_Independent_Binding_Sites.xml');

% The project contains a model, m1, of the binding of a ligand to two
% receptors R and R1 with different Kd's.

sbioselect(m1,'Type','compartment') 
sbioselect(m1,'Type','species')
sbioselect(m1,'Type','parameter')
sbioselect(m1,'Type','reaction')
getequations(m1)
 
% We extract the parameters from the model and save them as variables in
% the workspace.

C = sbioselect(m1,'Name','Cell');
L = sbioselect(m1,'Name','L');
R = sbioselect(m1,'Name','R');
R_1 = sbioselect(m1,'Name','R_1');
LR = sbioselect(m1,'Name','LR');
LR_1 = sbioselect(m1,'Name','LR_1');
LR_TOT = sbioselect(m1,'Name','LR_TOT');
kon = sbioselect(m1,'Name','kon');
koff = sbioselect(m1,'Name','koff');
kon_R_1 = sbioselect(m1,'Name','kon_R_1');
koff_R_1 = sbioselect(m1,'Name','koff_R_1');
 
% Before we start the simulation we set all the necessary parameters and
% variable values.
 
C.Capacity = 1;
R.InitialAmountUnits = 'micromole/liter';
R_1.InitialAmountUnits = 'micromole/liter';
L.InitialAmountUnits = 'micromole/liter';
LR.InitialAmountUnits = 'micromole/liter';
LR_1.InitialAmountUnits = 'micromole/liter';
 
% First we set a very high concentration of the ligand to observe the
% complete saturation of the receptor.

L.InitialAmount = 6000;
R.InitialAmount = 0.5;
R_1.InitialAmount = 0.5;
LR.InitialAmount = 0;
LR_1.InitialAmount = 0;
LR_TOT.InitialAmount = 0;
 
% The following is the default rule: both LR and LR_1 absorb the same in
% the spectrophotometer.
 
set(m1.rule,'Rule','LR_TOT = LR + LR_1');
 
% What happens if we change the rule? For example, LR_1 could absorb less
% than LR. set(m1.rule,'Rule','LR_TOT = 0.5*LR + 1*LR_1');
 
% Here we store the initial parameters.

C_Init = 1;
R_Init = 0.5;
R_1_Init = 0.5;
LR_Init = 0;
LR_1_Init = 0;
LR_TOT_Init = 0;
 
kon.Value = 10;
koff.Value = 100; % Kd(R) = 10 µM
kon_R_1.Value = 10;
koff_R_1.Value = 3000; % Kd(R_1) = 300 µM
 
% We also get some information on the configuration parameters of the
% simulation, and we set the initial stop time. For example, we will use 2
% ms as the stop time and will set the solver to ODE23t (trapezoidal
% solver, moderatively stiff) or ODE15s (RK solver, stiff).
 
configset_m1 = getconfigset(m1);
configset_m1
 
Stop = 0.0002;
set(configset_m1, 'StopTime', Stop);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);
set(configset_m1, 'SolverType', 'ode15s');

% Now we can run a small simulation with the existing parameters. Only 2
% milliseconds (msec) are necessary for the system to reach the equilibrium
% (steady-state). We recall here that since we are simulating the reaction
% from the ideal time in which there is an instantaneous mixing of ligand
% and receptor (time 0) until the time in which the equilibrium is reached
% (steady-state), the result is a reaction curve that mimics an experiment
% of pre-steady-state kinetics.

Binding_Kinetics = sbiosimulate(m1);

% We can get some information about the simulation. If we want to see the
% actual numbers, the time points are in the array Binding_Kinetics.Time
% and the the time-course of the simulation is in the
% Binding_Kinetics.Data. The names of the species in the reaction are in
% Binding_Kinetics.DataNames.

get(Binding_Kinetics)
openvar Binding_Kinetics.Time
openvar Binding_Kinetics.Data
openvar Binding_Kinetics.DataNames
 
% Species are logged in the following order: L R LR LR_1 LR_TOT R_1. We can
% obtain a fast plot of all the species in the simulation using the
% command:

% sbioplot(Binding_Kinetics)

% but it makes more sense to plot only some of the species:

%%
Two_binding_sites_1 = figure;
set(Two_binding_sites_1,'Units','normalized','Position',[0.4 0.0 0.6 1.0],...
    'Name','Binding Curve ');
subplot1 = subplot(2,1,1,'Parent',figure(gcf));
box(subplot1,'on');
grid(subplot1,'on');
hold(subplot1,'all');
plot(Binding_Kinetics.Time,Binding_Kinetics.Data(:,[2 6 3:5] ))
legend('Free R','Free R\_1','Bound R','Bound R\_1','Total Bound R','Location','Best')
set(gca,'YLim',[-0.05 1.05]);set(gca,'XLim',[-0.00001 1.05*Stop]);
ylabel('Species Concentration (µM)');xlabel('Time (s)')
title('Progress Curve with [L] = 6 mM');

% At this concentration R_1 is only slightly less saturated than R.
% However, it is even more important to know what is the saturation level
% and how long it takes for the reaction to reach equilibrium at very low
% concentration of the ligand.
 
C.Capacity = 1;
L.InitialAmount = 1;
R.InitialAmount = 0.5;
R_1.InitialAmount = 0.5;
LR.InitialAmount = 0;
LR_1.InitialAmount = 0;
LR_TOT.InitialAmount = 0;
 
Stop = 0.05;
set(configset_m1, 'StopTime', Stop);
 
Binding_Kinetics = sbiosimulate(m1);

subplot2 = subplot(2,1,2,'Parent',figure(gcf));
box(subplot2,'on');
grid(subplot2,'on');
hold(subplot2,'all');
plot(Binding_Kinetics.Time,Binding_Kinetics.Data(:,[1:2 6 3:5]))
legend('Free L','Free R','Free R\_1','Bound R','Bound R\_1','Total Bound R','Location','Best')
set(gca,'YLim',[-0.05 1.05]);set(gca,'XLim',[-1E-4 1.05*Stop]);
ylabel('Species Concentration (µM)');xlabel('Time (s)')
title('Progress Curve with [L] = 1 µM');
 
% It is quite apparent that at low concentration of ligand only the
% high-affinity site (receptor R) is saturated.

%% Generation of data with random errors
close all
clear Bound Unbound Bound_R Bound_R_1 

% Now we are ready to simulate an entire equilibrium binding experiment. We
% recall that Kd(R) = 10 µM and Kd(R_1) = 300 µM. In order to examine the
% effect of different concentrations of the ligand it is convenient to
% write a 'loop'. First we define a vector with all the concentrations we
% want to explore spaced logarithmically. It is a good idea to cover a
% range from 1/10 to 20-30 times the Kd value(s) in order to completely
% saturate the receptor(s). Thus, we will cover the range of concentrations
% 1-6000 µM. We also extend the simulation time to 0.5 seconds in order to
% be absolutely certain that under all concentrations of the ligand we have
% reached equilibrium.

Stop = 0.5;
set(configset_m1, 'StopTime', Stop);

% Finally, to make the simulation more realistic we can add a normally
% distributed random error to the volume of the binding assay and to the
% concentrations of receptor and ligand:

High_Conc = 6000;
conc_vec = logspace(log10(1),log10(High_Conc),35);
range = conc_vec <= 1 + High_Conc;

Bound = zeros(length(conc_vec),1);
Bound_R = zeros(length(conc_vec),1);
Bound_R_1 = zeros(length(conc_vec),1);
Unbound = zeros(length(conc_vec),1);
time = Binding_Kinetics.Time;
 
for i = 1:length(conc_vec)
    error = (rand(3,1)-0.5);
    error(1) = 1 + error(1)*0.15;
    error(2) = 1 + error(2)*0.15;
    error(3) = 1 + error(3)*0.15;
    set(C,'Capacity',C_Init*error(1));
    set(R,'InitialAmount',R_Init*error(2)/C.Capacity);
    set(R_1,'InitialAmount',R_1_Init*error(2)/C.Capacity);
    set(L,'InitialAmount',conc_vec(i)*error(3)/C.Capacity);
    set(LR,'InitialAmount',0);    
    set(LR_1,'InitialAmount',0);    
    set(LR_TOT,'InitialAmount',0);    
    Binding_Kinetics = sbiosimulate(m1);
    Bound(i) = Binding_Kinetics.Data(end,5);
    Bound_R(i) = Binding_Kinetics.Data(end,3);
    Bound_R_1(i) = Binding_Kinetics.Data(end,4);
    Unbound(i) = Binding_Kinetics.Data(end,1);
end

%% Binding curve
% Assuming these synthetic data are a realistic representation of a real
% equilibrium binding experiment, we can now derive the binding parameters
% from them using traditional tool, and see if they correspond well to the
% parameters used to generate the data. First we plot the binding curve:
% each point is the value of the bound ligand at steady-state in a
% simulation in which the ligand concentration is changed to progressively
% higher values.
clear X Y X1 X2 Y1 Y2

Two_binding_sites_2 = figure;
    set(Two_binding_sites_2,'Units','normalized','Position',[0.4 0.0 0.6 1.0],...
        'Name','Binding Curve Fit');clf
subplot1 = subplot(3,1,1,'Parent',figure(gcf));
box(subplot1,'on');grid(subplot1,'on');hold(subplot1,'all');
plot(Unbound(range),Bound(range),'--or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5);
plot(Unbound(range),Bound_R(range),':xg',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5);
plot(Unbound(range),Bound_R_1(range),':+b',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5);
set(gca,'YLim',[-0.05 1.15]);set(gca,'XLim',[-100 Unbound(end)+100]);
legend('Bound Total','Bound R','Bound R\_1','Location','Best');
xlabel('[Free Ligand]');ylabel('[Bound Ligand]');title('Binding Curve ');
 
%% Scatchard fit
% We can determine the binding parameters with a Scatchard plot of the
% ratio of [bound]/[unbound (free)] vs the [bound] ligand concentrations.
% In this case we have:

% Y = [Bound]/[Free] X = [Bound] Slope = -1/Kd. X axis intercept = [R0]
% (total receptor conc.)

% First we fit the data corresponding to the high affinity receptor by
% linear regression (least square fit) to the equation of a straight line.
% We recall here that since slope = -1/Kd, if the Kd is small (high
% affinity) the slope is going to be steeper: this is the left-hand side of
% the Scatchard plot. We are going to fit a straight line through the first
% 9 points of the Scatchard. Then, separately we fit the data corresponding
% to the low affinity receptor (right hand side of the plot).
 
% Ranges of concentrations that will be used in the fit of the Scatchard
% plot.
s_range = range;
s_range_1 = true(1,35);s_range_1(19:end) = false;
s_range_2 = true(1,35);s_range_2(1:22) = false;

subplot2 = subplot(3,1,2,'Parent',figure(gcf));
box(subplot2,'on');grid(subplot2,'on');hold(subplot2,'all');
X = Bound(s_range);
Y = Bound(s_range)./Unbound(s_range);
plot(X,Y,'s','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5);

% High affinity range
X1 = Bound(s_range_1);
Y1 = Bound(s_range_1)./Unbound(s_range_1);
 
f = fittype('a*x + b');
[Scatchard_1,GOF_1] = fit(X1,Y1,f,'StartPoint',[-0.01 0.05]);
plot(Scatchard_1,'-b');
ylim([0,Scatchard_1(0)]);
xlim([0,1.1]);

Scatchard_1
GOF_1
Scatchard_1_params = coeffvalues(Scatchard_1);
Kd_1 = -1/Scatchard_1_params(1);
Receptor_conc_1 = -Scatchard_1_params(2)/Scatchard_1_params(1);
 
% Low affinity range
X2 = Bound(s_range_2);
Y2 = Bound(s_range_2)./Unbound(s_range_2);
 
f = fittype('a*x + b');
[Scatchard_2,GOF_2] = fit(X2,Y2,f,'StartPoint',[-0.03 0.03]);
plot(Scatchard_2,'-r');
 
Scatchard_2
GOF_2
Scatchard_2_params = coeffvalues(Scatchard_2);
Kd_2 = -1/Scatchard_2_params(1);
Receptor_conc_2 = -Scatchard_2_params(2)/Scatchard_2_params(1);
 
% This result is likely to be incorrect for the low affinity site because
% of the contribution from the high affinity site.
 
legend('Scatchard','Fit\_1','Fit\_2');
string0 = ' µM';
string1 = 'Kd\_1 = ';
string2 = num2str(Kd_1,'%6.2e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%6.2e\n');
 
annotation(Two_binding_sites_2,'textbox',...
    [0.58 0.48 0.12 0.09],...    
    'String',{[string1 string2 string0; string3 string4 string0;]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.8],...
    'Color',[1 0 0]);
 
xlabel('[Bound Ligand]');ylabel('[Bound Ligand]/[Free Ligand]');title('Scatchard plot ');

%% Hyperbola fit
% Alternatively,  we can fit directly the plot of the fractional saturation
% of the receptor (FR) vs the amount of free ligand (L) with the hyperbolic
% function we presented at the beginning of the example, rearranged in
% terms of fractional saturation:

% [LR] = [R0]*[L]/([L]+Kd)
% [LR]/[R0] = [L]/([L]+Kd)

% Two hyperbolas are fitted each contributing 50% to the curve. For a
% single hyperbola it would be: f = fittype('x/(a + x)'); For two
% hyperbolas we use a two-component fit with 2 unknowns. 'a' and 'b' are
% the two Kd's. We will use the values of Kd_1 and Kd_2 found by the
% Scatchard plot as the starting points for the fits. We can also include
% some lower and upper bounds for the fitting parameters.

FR = Bound/(R_Init + R_1_Init);X = Unbound(s_range);Y = FR(s_range);
 
subplot3 = subplot(3,1,3,'Parent',figure(gcf));
box(subplot3,'on');grid(subplot3,'on');hold(subplot3,'all');
 
f = fittype('(0.5*x/(a + x)) + (0.5*x/(b + x))');
% [Hyperb,GOF] = fit(X,Y,f,'StartPoint',[Kd_1 Kd_2],...
%     'Lower',[0.3*Kd_1 0.3*Kd_2],'Upper',[1.7*Kd_1 1.7*Kd_2]);
[Hyperb,GOF] = fit(X,Y,f,'StartPoint',[Kd_1 Kd_2]);
 
plot(X,Hyperb(X),'-r');
plot(X,Y,'o','MarkerEdgeColor','b','MarkerFaceColor','c','MarkerSize',5);
ylim([0,1.15]); xlim([-100,Unbound(end)+100]);
legend('Rbound/Rtotal','Fit','Location','Best');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

Hyperb
GOF
Hyperb_params = coeffvalues(Hyperb);
Kd_1 = Hyperb_params(1);
Kd_2 = Hyperb_params(2);

string0 = ' µM';
string1 = 'Kd\_1 = ';
string2 = num2str(Kd_1,'%6.3e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%6.3e\n');

% Create textbox
annotation(Two_binding_sites_2,'textbox',...
    [0.54 0.18 0.12 0.08],...
    'String',{[string1 string2 string0 ; string3 string4 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.8],...
    'Color',[1 0 0]);
title('Hyperbola fit ');
 
%% Conclusion
% The presence of multiple binding sites can be easily inferred from the
% presence of a biphasic Scatchard plot. However, low affinity sites are
% not determined well by Scatchard. In this case, the best determination of
% the Kd's of the individual binding sites is obtained by fitting a 2 (or
% more components) hyperbola to the fractional saturation plot.
 
