% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% Binding of a Ligand to Two Cooperative Binding Sites in the Same Receptor 

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

% Copyright 2015 Domenico Gatti.


%% Aims
%
% * Learn how to represent a binding experiment with two cooperative
%   binding sites.
%
% * Observe the effect of cooperativity on the binding curve in the 
%   low-concentration range.
% 
% * Learn how to recognize the distortion produced in the Scatchard plot by 
%   the presence of two cooperative binding sites.
%

%% Background
% The best way to simulate two cooperative binding sites in a single
% receptor is to equate the system to two receptors each with one binding
% site. Then, the Kd of one site is represented as a function of the Kd in
% the other site. Also in this case the total conncentration of the receptor
% is given by the sum of the two receptors.
%%
% We start by loading the project 

% sbioloadproject('Two_Cooperative_Binding_Sites');

% or alternatively:

close all
clear, clc

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

% and adding some programs to the path:
% addpath(genpath('vhline'));

% We also set a file name to save our simulation results:
savefile = 'Two_Cooperative_Binding_Sites.mat';

%%
% The project contains a model, |m1|, of the binding of a ligand to two
% receptors with different Kd's. We can also load the model interactively
% by starting SimBiology desktop with |simbiology| and opening the
% project file Two_Cooperative_Binding_Sites.sbproj.

%%
% Get information about the model.
sbioselect(m1,'Type','compartment')
sbioselect(m1,'Type','species')
sbioselect(m1,'Type','parameter')
sbioselect(m1,'Type','reaction')

% Let's also take a look at the equations that describe the model:
getequations(m1)

%% Cooperativity
% Notice how positive cooperativity between the two site is implemented in
% just two lines of code:
%  
% ReactionFlux1 = (kon_R*[L]*[R] - koff_R*[LR])*([R_1]/([R_1]+[LR_1])) - koff_R_B*[LR]*([LR_1]/([R_1]+[LR_1])))*Cell
%
% ReactionFlux2 = (kon_R_1*[L]*[R_1] - koff_R_1*[LR_1])*([R]/([R]+[LR])) - koff_R_1B*[LR_1]*([LR]/([R]+[LR])))*Cell
%
% Without cooperativity it would be:
%
% ReactionFlux1 = (kon_R*[L]*[R] - koff_R*[LR])*Cell
%
% ReactionFlux2 = (kon_R_1*[L]*[R_1] - koff_R_1*[LR_1])*Cell
%
%%
% As usual, we will 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_R = sbioselect(m1,'Name','kon_R');
koff_R = sbioselect(m1,'Name','koff_R');
kon_R_B = sbioselect(m1,'Name','kon_R_B');
koff_R_B = sbioselect(m1,'Name','koff_R_B');
kon_R_1 = sbioselect(m1,'Name','kon_R_1');
koff_R_1 = sbioselect(m1,'Name','koff_R_1');
kon_R_1B = sbioselect(m1,'Name','kon_R_1B');
koff_R_1B = sbioselect(m1,'Name','koff_R_1B');

%% Time evolution of the system
% We will simulate the time evolution of this system. 

% Before we start the simulation we will 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';

L.InitialAmount = 600;
R.InitialAmount = 0.5;
R_1.InitialAmount = 0.5;
LR.InitialAmount = 0;
LR_1.InitialAmount = 0;
LR_TOT.InitialAmount = 0;

% 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;

% Rate law for R:
% set(m1.reactions(1), 'ReactionRate',...
%   'kon_R*L*R - koff_R*LR');
set(m1.reactions(1), 'ReactionRate',...
  'kon_R*L*R - koff_R*LR*(R_1/(R_1+LR_1)) - koff_R_B*LR*(LR_1/(R_1+LR_1))');

% Rate law for R_1:
% set(m1.reactions(2), 'ReactionRate',...
%   'kon_R_1*L*R_1 - koff_R_1*LR_1');
set(m1.reactions(2), 'ReactionRate',...
  'kon_R_1*L*R_1 - koff_R_1*LR_1*(R/(R+LR)) - koff_R_1B*LR_1*(LR/(R+LR))');

% Modulated site R: (Kd(R) = Kd(R_1) 
% Without cooperativity 
kon_R.Value = 10;
koff_R.Value = 300; % Kd(R) = 30 µM
% With cooperativity
kon_R_B.Value = 10; 
koff_R_B.Value = 30; % Kd(R)cooper = 3 µM

Kd = koff_R.Value/kon_R.Value;
display(['Kd = ' num2str(Kd) ' µM']);

% Modulated site R1:
% Without cooperativity 
kon_R_1.Value = 10;
koff_R_1.Value = 300; % Kd(R_1) = 30 µM
% With cooperativity
kon_R_1B.Value = 10; 
koff_R_1B.Value = 30; % Kd(R_1)cooper = 3 µM

Kd_cooper = koff_R_1B.Value/kon_R_1B.Value;
display(['Kd_cooper = ' num2str(Kd_cooper) ' µM']);

%%
% We also need to get some information on the configuration parameters of
% the simulation.

configset_m1 = getconfigset(m1);
configset_m1
% We set the initial stop time at 0.2 mseconds
Stop = 0.002;
set(configset_m1, 'StopTime', Stop);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

% Here we simulate.
Binding_Kinetics = sbiosimulate(m1);

% Let's get the order in which the species are logged:
Binding_Kinetics.DataNames

% Species are logged in the following order: L R LR LR_1 LR_TOT R_1. We
% will be plotting the free receptor R and the total bound receptor LR_TOT.

%%
FIG_1 = figure;
    set(FIG_1,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','Simple Binding');clf
    
axes1 = axes('Parent',FIG_1,...
    'Position',[0.08 0.1 0.9 0.82]);
hold(axes1,'all');
box(axes1,'on');
plot(Binding_Kinetics.Time,Binding_Kinetics.Data(:,[2 6 3 4 5]))
legend('R','R\_1','LR','LR\_1','LR\_TOT','Location',[0.8 ...
    0.6 0.1 0.2])
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');

% Notice how LR and LR_1 rise together. This would not be the case if there
% was not cross-cooperativity between the two sites.

%%
%
% set(gca,'XLim',[-0.00005 0.001],'YLim',[-0.05 0.7]);
%
% legend('R','R_1','LR','LR\_1','LR\_TOT','Location','Best');
%

%%
% We repeat considering the other possible extreme (1 µM) in the
% concentration of the ligand.

close all

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);

%%
FIG_2 = figure;
    set(FIG_2,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','Simple Binding');clf
    
axes1 = axes('Parent',FIG_2,...
    'Position',[0.08 0.1 0.9 0.82]);
hold(axes1,'all');
box(axes1,'on');
plot(Binding_Kinetics.Time,Binding_Kinetics.Data(:,[2 6 3 4 5]))
legend('R','R\_1','LR','LR\_1','LR\_TOT','Location','Best')
set(gca,'YLim',[-0.05 0.55]);
set(gca,'XLim',[-0.0001 1.05*Stop]);
ylabel('Species Concentration (µM)')
xlabel('Time (s)')
title('Progress Curve with [L] = 10 µM');

%% Simulation of the binding curve
close all
clear Bound Unbound Bound_R Bound_R_1

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;

% 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):

High_Conc = 600;
conc_vec = logspace(log10(1),log10(High_Conc),35);
% openvar conc_vec
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);
% Bound_R_2 = zeros(length(conc_vec),1);
Unbound = zeros(length(conc_vec),1);
time = Binding_Kinetics.Time;

% We also extend the simulation time to 1 second in order to be certain
% that under all concentrations of the ligand we have reached the binding
% equilibrium.

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

for i = 1:length(conc_vec)
    set(R,'InitialAmount',0.5);    
    set(R_1,'InitialAmount',0.5);    
    set(LR,'InitialAmount',0);    
    set(LR_1,'InitialAmount',0);    
    set(LR_TOT,'InitialAmount',0);    
    set(L,'InitialAmount',conc_vec(i));
    Binding_Kinetics = sbiosimulate(m1);
    % In this case we will take the 5th colums, which is the total bound R.
    Unbound(i) = Binding_Kinetics.Data(end,1);
    Bound_R(i) = Binding_Kinetics.Data(end,3);
    Bound_R_1(i) = Binding_Kinetics.Data(end,4);
    Bound(i) = Binding_Kinetics.Data(end,5);
end

%%
% Here we plot all the results obtained with theoretical data.
%
clear X Y X1 X2 Y1 Y2

COOPERATIVE_BINDING = figure;
    set(COOPERATIVE_BINDING,'Units','normalized','Position',[0.4 0.0 0.6 1.0],...
        'Name','Binding Curve Fit');clf
    subplot1 = subplot(3,2,1,'Parent',figure(gcf));
    box(subplot1,'on');
    grid(subplot1,'on');
    hold(subplot1,'all');

% Dose-Response plots (percentage of maximal response against log of the
% ligand concentration):

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.1]);
set(gca,'XLim',[-5 5+conc_vec(end)]);
         
legend('LR\_TOT','LR','LR\_1','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve with error-free data');

%% Sigmoidicity of the binding curve in the low concentration range
% No sigmoidicity of the binding curve in the low concentration range.

set(gca,'XLim',[-5 50],'YLim',[-0.05 0.5]);
legend('LR\_TOT','LR','LR\_1','Location','Best');

%%
% and back.
set(gca,'XLim',[-5 600],'YLim',[-0.05 1.1]);
legend('LR\_TOT','LR','LR\_1','Location','Best');

%% Scatchard plot

    subplot2 = subplot(3,2,2,'Parent',figure(gcf));
    box(subplot2,'on');
    grid(subplot2,'on');
    hold(subplot2,'all');

s_range = range;    
s_range_1 = true(1,35);
s_range_1(15:end) = false;
s_range_2 = true(1,35);
s_range_2(1:29) = false;
        
X = Bound(s_range);
Y = Bound(s_range)./Unbound(s_range);
plot(X,Y,'-ys',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
title('Scatchard Plot');

% First we fit the data corresponding to the left hand side of the plot:
% this represents the low affinity stage of the binding.

X1 = Bound(s_range_1);
Y1 = Bound(s_range_1)./Unbound(s_range_1);

% Set the plot limits
ylim([0,0.04]);
xlim([0,1.1]);

f = fittype('a*x + b');
[Scatchard_1,GOF_1] = fit(X1,Y1,f,'StartPoint',[-0.004 0.035]);
plot(Scatchard_1,'-b');
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

% Retrieve the fit parameters.
Scatchard_1
GOF_1
% or simply:
Scatchard_1_params = coeffvalues(Scatchard_1);
Kd_1 = -1/Scatchard_1_params(1);
Receptor_conc_1 = -Scatchard_1_params(2)/Scatchard_1_params(1);
         
% Next we fit the data corresponding to the right hand side of the plot:
% this represents the high affinity stage of the binding.

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.3 0.3]);
plot(Scatchard_2,'-r');
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

% Retrieve the fit parameters.
Scatchard_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);
% 

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

xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

% Create textbox
annotation(COOPERATIVE_BINDING,'textbox',...
    [0.59 0.8 0.12 0.040],...
    'String',{[string1 string2 string0; string3 string4 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);
title('Scatchard Plot');
legend('Scatchard','Fit\_1: Kd\_1','Fit\_2: Kd\_2','Location','Best');

% We can obtain a good measure of the low affinity (unmodulated) Kd_1 from
% the end-point of the Scatchard plot:

Kd_1_best = Receptor_conc_2/Scatchard_1_params(2);

%% Data with noise
%
% close all
clear Bound Unbound Bound_R Bound_R_1 X Y X1 X2 Y1 Y2

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.05;
    error(2) = 1 + error(2)*0.05;
    error(3) = 1 + error(3)*0.05;
    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

%% 
    subplot3 = subplot(3,2,3,'Parent',figure(gcf));
    box(subplot3,'on');
    grid(subplot3,'on');
    hold(subplot3,'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.1]);
set(gca,'XLim',[-5 5+conc_vec(end)]);

legend('LR\_TOT','LR','LR\_1','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve with random error-data');

%%
% Scatchard Plot 
    subplot4 = subplot(3,2,4,'Parent',figure(gcf));
    box(subplot4,'on');
    grid(subplot4,'on');
    hold(subplot4,'all');

X = Bound(s_range);
Y = Bound(s_range)./Unbound(s_range);
plot(X,Y,'-ys',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);

%--------------------------------------------------------------------------
% We repeat everything as we did for the theoretical data.
%--------------------------------------------------------------------------

s_range = range;    
s_range_1 = true(1,35);
s_range_1(16:end) = false;
s_range_2 = true(1,35);
s_range_2(1:20) = false;

% High affinity range
X1 = Bound(s_range_1);
Y1 = Bound(s_range_1)./Unbound(s_range_1);

plot(X1,Y1,'yo',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','none',...
             'MarkerSize',8);

f = fittype('a*x + b');
% options = fitoptions(f);
% options.StartPoint = [-0.004 0.035];
[Scatchard_1,GOF_1] = fit(X1,Y1,f,'Robust','Bisquare','StartPoint',[-0.004 0.035]);
plot(Scatchard_1,'-b');

% Retrieve the fit parameters.
Scatchard_1
GOF_1
% or simply:
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);

plot(X2,Y2,'yo',...
             'MarkerEdgeColor','r',...
             'MarkerFaceColor','none',...
             'MarkerSize',8);

f = fittype('a*x + b');
[Scatchard_2,GOF_2] = fit(X2,Y2,f,'Robust','Bisquare','StartPoint',[-0.3 0.5]);
plot(Scatchard_2,'-r');

% Retrieve the fit parameters.
Scatchard_2
GOF_2
% or simply:
Scatchard_2_params = coeffvalues(Scatchard_2);
Kd_2 = -1/Scatchard_2_params(1);
Receptor_conc_2 = -Scatchard_2_params(2)/Scatchard_2_params(1);

% Axes limits
ylim([0,0.04]);

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

annotation(COOPERATIVE_BINDING,'textbox',...
    [0.6 0.5 0.12 0.040],...
    'String',{[string1 string2 string0; string3 string4 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);


xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
title('Scatchard plot with random-error data');
legend('Scatchard','Fit\_1: points','Fit\_1: Kd\_1','Fit\_2: points','Fit\_2: Kd\_2','Location','Best');

%% Dose Response Curve
    subplot5 = subplot(3,2,5,'Parent',figure(gcf),'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
    box(subplot5,'on');
    grid(subplot5,'on');
    hold(subplot5,'all');
    semilogx(Unbound(range),Bound(range),'--or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',subplot5);
    semilogx(Unbound(range),Bound_R(range),':xg',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',subplot5);
    semilogx(Unbound(range),Bound_R_1(range),':+b',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',subplot5);

set(gca,'YLim',[-0.05 1.1]);
set(gca,'XLim',[1 10+conc_vec(end)]);

% Expected first site:
vline([3 3],{'g--','LineWidth',2});
% Expected second site:
vline([30 30],{'c--','LineWidth',2});
% Apparent Kd
vline([16.6 16.6],{'m-','LineWidth',2});

legend('LR\_TOT','LR','LR\_1','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Dose-Response Curve with random error-data');

%% Hyperbola fitting
clear X Y
% As usual, we fit directly the plot of the fractional saturation of the
% receptor (FR) vs the amount of free ligand (L).

FR = Bound/(R_Init + R_1_Init);
X = Unbound(s_range);
Y = FR(s_range);

subplot6 = subplot(3,2,6,'Parent',figure(gcf));
box(subplot6,'on');
grid(subplot6,'on');
hold(subplot6,'all');

% First we try to fit with just one component.
f = fittype('b*(x/(a + x))');
[Hyperb0,GOF0] = fit(X,Y,f,'StartPoint',[Kd_2 1]);

% Retrieve the fit parameters.
Hyperb0
GOF0
% or simply:
Hyperb0_params = coeffvalues(Hyperb0);
Kd_0 = Hyperb0_params(1);

% Then we use a two-component fit with 4 unknowns.         
f = fittype('b*(x/(a + x)) + d*(x/(c + x))');
% [Hyperb1,GOF1] = fit(X,Y,f,'StartPoint',[Kd_2 0.5 Kd_1 0.5]);
[Hyperb1,GOF] = fit(X,Y,f,'StartPoint',[Kd_1 0.5 Kd_2 0.5],...
    'Lower',[0.6*Kd_1 0.3 0.6*Kd_2 0.3],'Upper',[1.4*Kd_1 0.7 1.4*Kd_2 0.7]);

plot(X,Hyperb0(X),'-y');
plot(X,Hyperb1(X),'-r');
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);

ylim([-0.05,1.1]); 
xlim([-50,Unbound(end)+50]);
legend('1 Comp. Fit: Kd\_1','2 Comp. Fit: Kd\_2, Kd\_3','Rbound/Rtotal','Location','SouthEast');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

% Retrieve the fit parameters.
Hyperb1
GOF
% or simply:
Hyperb1_params = coeffvalues(Hyperb1);
Kd_2 = Hyperb1_params(1);
Kd_3 = Hyperb1_params(3);

string1 = 'Kd\_1 = ';
string2 = num2str(Kd_0,'%6.2e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%6.2e\n');
string5 = 'Kd\_3 = ';
string6 = num2str(Kd_3,'%6.2e\n');
string0 = ' µM';

% Create textbox
annotation(COOPERATIVE_BINDING,'textbox',...
    [0.77 0.21 0.12 0.05],...
    'String',{[string1 string2 string0; string3 string4 string0;...
    string5 string6 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

title('Hyperbola fit with random-error data');

%% Hill coefficient
% It is used to determine the degree of cooperativity of the ligand
% binding to the receptor. It is derived from the Hill equation, in which
% it is the 'n':
%
% log(theta/(1-theta)) = n*log[L] - log(Kd)
%
% n>1 - Positive cooperativity: binding of the ligand to one site increases
%       the affinity of the other site for the ligand.
%
% n<1 - Negative cooperativity: binding of the ligand to one site decreases
%       the affinity of the other site for the ligand.
%
% n=1 - Non-cooperativity: binding sites are independent.
%
%
% Two ways of determining the Hill coefficient:
%
% * 1. Traditional Hill plot. 
%
% It is a plot of the Hill equation in the form:
%
% log[y/(1-y)] against log[x] 
% 
% where y is the fractional saturation with ligand and [x] is the free
% (unbound) ligand concentration. The slope is the Hill coefficient, and
% the intercept with the ordinate (Y axis) is equal to -log(Kd).
%
% IMPORTANT: it is often reported that the intercept with the X axis
% represents the Kd because at Y = 0 theta = 0.5. However, this is only
% true if n = 1. In the more general case, solving for the x intercept
% (when the y axis variable is 0) in the Hill equations gives:
% 
% 0 = n*log[L] -log(Kd)
%
% n*log[L] = log(Kd)
%
% log[L] = (log(Kd))/n
%

HILL_1 = figure;
X = log(Unbound(s_range));
Y = log(FR(s_range)./(1 - FR(s_range)));
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
xlabel('log[Free Ligand]');
ylabel('log(Frac. Sat./(1-Frac. Sat.)');
title('Hill plot');
         
         
hold on
f = fittype('a*x + b');
[Hill_1,GOF_1] = fit(X(1:12),Y(1:12),f,'StartPoint',[1 -4],'Robust','on');
Hill_1
plot(Hill_1,'-b');
[Hill_2,GOF_2] = fit(X(20:35),Y(20:35),f,'StartPoint',[1 1],'Robust','on');
Hill_2
plot(Hill_2,'-r');
[Hill_3,GOF_3] = fit(X(14:19),Y(14:19),f,'StartPoint',[2 -8],'Robust','on');
Hill_3
plot(Hill_3,'-g');

ylim([-6 6]);
legend('Binding Curve','Hill Fit1: Kd\_1','Hill Fit2: Kd\_2',...
    'Hill Fit3: Kd\_3','Location','NorthWest');

% Retrieve the fit parameters.
Hill_1_params = coeffvalues(Hill_1);
Hill_2_params = coeffvalues(Hill_2);
Hill_3_params = coeffvalues(Hill_3);
nH_1 = Hill_1_params(1);
nH_2 = Hill_2_params(1);
nH_3 = Hill_3_params(1);
Kd_1 = exp(-Hill_1_params(2));
Kd_2 = exp(-Hill_2_params(2));
Kd_3 = exp(-Hill_3_params(2));

string1 = 'nH\_1 = ';
string2 = num2str(nH_1,'%6.3e\n');
string3 = 'nH\_2 = ';
string4 = num2str(nH_2,'%6.3e\n');
string5 = 'nH\_3 = ';
string6 = num2str(nH_3,'%6.3e\n');

% Create textbox
annotation(HILL_1,'textbox',...
    [0.65 0.4 0.22 0.13],...
    'String',{[string1 string2 ; string3 string4 ; string5 string6]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

string1 = 'Kd\_1 = ';
string2 = num2str(Kd_1,'%6.2e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%6.2e\n');
string5 = 'Kd\_3 = ';
string6 = num2str(Kd_3,'%6.2e\n');

% Create textbox
annotation(HILL_1,'textbox',...
    [0.60 0.2 0.24 0.13],...
    'String',{[string1 string2 string0; string3 string4 string0; string5 string6 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

xlabel('log[Free Ligand]');
ylabel('log(Frac. Sat./(1-Frac. Sat.)');
title('Hill plot');

%%
% * 2. An alternative way is with with a logistic fit of 
% the data. A Dose-Response Logistic fit will give the Hill coefficient.
% This fit is based on the following equation:
%
% Response = (Max-Min)/(1+(Kd/x)^n  
%
% in which n = nH = Hill coefficient.
%
% Notice that theoretically nH is not supposed to be any higher than the
% number of binding sites which are occupied simultaneously. In the
% equation below a = Kd and b = nH.
%
clear X Y X1 X2 Y1 Y2 H H2

HILL_2 = figure;
    set(HILL_2,'Units','normalized','Position',[0.5 0.4 0.5 0.6],...
        'Name','Binding Curve Fit');clf
    subplot1 = subplot(2,2,1,'Parent',figure(gcf));
    box(subplot1,'on');
    grid(subplot1,'on');
    hold(subplot1,'all');

Y = Bound(s_range);
X = Unbound(s_range);

% 1 component (4 parameters) logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b)');
[Hill_1,GOF_1] = fit(X,Y,f,'StartPoint',[ 30 2 1 0 ]);
Hill_1
GOF_1

plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5,'Parent',subplot1);
plot(Hill_1,'-r');

% xlim([200 100]);
ylim([-0.05 1.1]);
legend('Binding Curve','Logistic Fit','location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit with 1 component');

%%
% Here we replot the fit as a Dose-Response curve.

    subplot2 = subplot(2,2,2,'Parent',figure(gcf),'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
    box(subplot2,'on');
    grid(subplot2,'on');
    hold(subplot2,'all');

% Parameters    
hpar = coeffvalues(Hill_1);

% 1 component (4 parameters) fit:

% Theoretical values for nH = 2
for i = 1:length(Unbound)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound(i))^2);
end
semilogx(Unbound,H2,'-r','Parent',subplot2);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
ylim([-0.05 1.1]);
title('Logistic fit as Dose-Response curve')

% Theoretical values for nH = 1
for i = 1:length(Unbound)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound(i))^1);
end
semilogx(Unbound,H2,'-g','Parent',subplot2);

% Experimental values
semilogx(Unbound,Hill_1(Unbound),'-b','Parent',subplot2);
semilogx(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5,'Parent',subplot2);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
ylim([-0.05 1.1]);

% Parameters    
hpar = coeffvalues(Hill_1);
string1 = 'nH_e_x_p = ';
string2 = num2str(hpar(2),'%4.2f\n');

legend('nH = 2','nH = 1',[string1 string2],'location','Best');

string3 = 'Kd = ';
string4 = num2str(hpar(1),'%4.2f\n');

% Create textbox
annotation(HILL_2,'textbox',...
    [0.61 0.75 0.12 0.035],...
    'String',{[string3 string4 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit as Dose-Response curve')

%%
    subplot3 = subplot(2,2,3,'Parent',figure(gcf));
    box(subplot3,'on');
    grid(subplot3,'on');
    hold(subplot3,'all');

Y = Bound(s_range);
X = Unbound(s_range);

% 2 components (8 parameters) logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b) + (g - h)/(1 + (e/x)^f)');
[Hill_2,GOF_2] = fit(X,Y,f,'StartPoint',[ Kd_1 nH_1 0.5 0 Kd_2 nH_2 0.5 0]);
% f = fittype('(c - d)/(1 + (a/x)^b) + (f - g)/(1 + (e/x)^b)');
% [Hill_2,GOF_2] = fit(X,Y,f,'StartPoint',[ Kd_1 nH_3 0.5 0 Kd_2 0.5 0]);
Hill_2
GOF_2

plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5,'Parent',subplot3);
plot(Hill_2,'-r');

% xlim([200 100]);
ylim([-0.05 1.1]);
legend('Binding Curve','Logistic Fit','location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit with 2 components');

%%
% Here we replot the fit as a Dose-Response curve.

    subplot4 = subplot(2,2,4,'Parent',figure(gcf),'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
    box(subplot4,'on');
    grid(subplot4,'on');
    hold(subplot4,'all');

% 2 components (8parameters) fit:
% Parameters    
hpar = coeffvalues(Hill_2);

% Theoretical values for nH = 2
for i = 1:length(Unbound)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound(i))^2) + ...
    (hpar(7) - hpar(8))/(1 + (hpar(5)./Unbound(i))^2);
end
semilogx(Unbound,H2,'-r','Parent',subplot4);

% Theoretical values for nH = 1
for i = 1:length(Unbound)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound(i))^1) + ...
    (hpar(7) - hpar(8))/(1 + (hpar(5)./Unbound(i))^1);
end
semilogx(Unbound,H2,'-g','Parent',subplot4);

% Experimental values
semilogx(Unbound,Hill_2(Unbound),'-b','Parent',subplot4);
semilogx(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5,'Parent',subplot4);

ylim([-0.05 1.1]);

% Parameters    
hpar = coeffvalues(Hill_2);
string1 = 'nH_e_x_p_1 = ';
string2 = num2str(hpar(2),'%4.2f\n');
string3 = 'nH_e_x_p_2 = ';
string4 = num2str(hpar(6),'%4.2f\n');

legend('nH = 2','nH = 1',[string1 string2 ; string3 string4],'location','Best');

xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit as Dose-Response curve')

%% Can you guess what negative cooperativity looks like?
% It is very easy to check what the binding curve would look if there was
% negative instead of positive cooperativity: in this case we will increase
% the value of koff_R_B and koff_R_1B rather than decreasing it.
%

% Modulated site R: (Kd(R) = Kd(R_1) 
% Without cooperativity 
kon_R.Value = 10;
koff_R.Value = 300; % Kd(R) = 30 µM
% With cooperativity
kon_R_B.Value = 10; 
koff_R_B.Value = 3000; % Kd(R)cooper = 300 µM

Kd = koff_R.Value/kon_R.Value;
display(['Kd = ' num2str(Kd) ' µM']);

% Modulated site R1:
% Without cooperativity 
kon_R_1.Value = 10;
koff_R_1.Value = 300; % Kd(R_1) = 30 µM
% With cooperativity
kon_R_1B.Value = 10; 
koff_R_1B.Value = 3000; % Kd(R_1)cooper = 300 µM

Kd_cooper = koff_R_1B.Value/kon_R_1B.Value;
display(['Kd_cooper = ' num2str(Kd_cooper) ' µM']);

%% 
% Simulation of the binding curve
close all
clear Bound Unbound Bound_R Bound_R_1

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;

% Usual loop:

High_Conc = 6000;
conc_vec = logspace(log10(1),log10(High_Conc),35);
% openvar conc_vec
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);
Bound_R_2 = zeros(length(conc_vec),1);
Unbound = zeros(length(conc_vec),1);

% We also extend the simulation time to 10 seconds in order to be certain
% that under all concentrations of the ligand we have reached the binding
% equilibrium.

Stop = 10;
set(configset_m1, 'StopTime', Stop);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

for i = 1:length(conc_vec)
    set(R,'InitialAmount',0.5);    
    set(R_1,'InitialAmount',0.5);    
    set(LR,'InitialAmount',0);    
    set(LR_1,'InitialAmount',0);    
    set(LR_TOT,'InitialAmount',0);    
    set(L,'InitialAmount',conc_vec(i));
    Binding_Kinetics = sbiosimulate(m1);
    % In this case we will take the 5th colums, which is the total bound R.
    Unbound(i) = Binding_Kinetics.Data(end,1);
    Bound_R(i) = Binding_Kinetics.Data(end,3);
    Bound_R_1(i) = Binding_Kinetics.Data(end,4);
    Bound(i) = Binding_Kinetics.Data(end,5);
end

time = Binding_Kinetics.Time;

%%
% Here we plot all the results obtained with theoretical data.
%
clear X Y X1 X2 Y1 Y2

NEG_COOPERATIVE_BINDING = figure;
    set(NEG_COOPERATIVE_BINDING,'Units','normalized','Position',[0.4 0.0 0.6 1.0],...
        'Name','Binding Curve Fit');clf
    subplot1 = subplot(3,2,1,'Parent',figure(gcf));
    box(subplot1,'on');
    grid(subplot1,'on');
    hold(subplot1,'all');

% Dose-Response plots

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.1]);
set(gca,'XLim',[-5 5+conc_vec(end)]);
         
legend('LR\_TOT','LR','LR\_1','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve');

%% 
% Notice how there is no sigmoidicity in the low concentration range of the
% binding curve.

set(gca,'XLim',[-5 300],'YLim',[-0.05 0.4]);
legend('LR\_TOT','LR','LR\_1','Location','Best');

%%
% and back:

set(gca,'YLim',[-0.05 1.1]);
set(gca,'XLim',[-5 5+conc_vec(end)]);
legend('LR\_TOT','LR','LR\_1','Location','Best');

%% 
% Scatchard plot

    subplot2 = subplot(3,2,2,'Parent',figure(gcf));
    box(subplot2,'on');
    grid(subplot2,'on');
    hold(subplot2,'all');

s_range = range;    
s_range_1 = true(1,35);
s_range_1(7:end) = false;
s_range_2 = true(1,35);
s_range_2(1:26) = false;
        
X = Bound(s_range);
Y = Bound(s_range)./Unbound(s_range);
plot(X,Y,'-ys',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
title('Scatchard Plot');

% First we fit the data corresponding to the high affinity receptor:
X1 = Bound(s_range_1);
Y1 = Bound(s_range_1)./Unbound(s_range_1);
plot(X1,Y1,'yo',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','none',...
             'MarkerSize',8);
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

f = fittype('a*x + b');
[Scatchard_1,GOF_1] = fit(X1,Y1,f,'StartPoint',[-0.006 0.03]);
plot(Scatchard_1,'-b');
xlim([0,1.1]);
ylim([0,0.03]);

% Retrieve the fit parameters.
Scatchard_1
GOF_1
% or simply:
Scatchard_1_params = coeffvalues(Scatchard_1);
Kd_1 = -1/Scatchard_1_params(1);
Receptor_conc_1 = -Scatchard_1_params(2)/Scatchard_1_params(1);

% Next we fit the data corresponding to the low affinity receptor:
X2 = Bound(s_range_2);
Y2 = Bound(s_range_2)./Unbound(s_range_2);
plot(X2,Y2,'yo',...
             'MarkerEdgeColor','r',...
             'MarkerFaceColor','none',...
             'MarkerSize',8);
xlim([0,1.1]);
ylim([0,0.03]);
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

f = fittype('a*x + b');
[Scatchard_2,GOF_2] = fit(X2,Y2,f,'StartPoint',[-0.004 0.005]);
plot(Scatchard_2,'-r');

% Retrieve the fit parameters.
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);

% Usual correction: 
%
LR_formed = (Unbound(s_range_2)*R_Init)./(Unbound(s_range_2)+Kd_1);
X3 = Bound(s_range_2) - LR_formed;
Y3 = (Bound(s_range_2)- LR_formed)./Unbound(s_range_2);

plot(X3,Y3,':yo',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
         
f = fittype('a*x + b');
[Scatchard_3,GOF_3] = fit(X3,Y3,f,'StartPoint',[-0.001 0.001]);
plot(Scatchard_3,'-c');

% Retrieve the fit parameters.
Scatchard_3
GOF_3
% or simply:
Scatchard_3_params = coeffvalues(Scatchard_3);
Kd_3 = -1/Scatchard_3_params(1);
Receptor_conc_3 = -Scatchard_3_params(2)/Scatchard_3_params(1);


% Modified correction: we take the high affinity Kd from the end-point of
% the Scatchard Plot
%
Kd_1_mod = Receptor_conc_2/Scatchard_1_params(2);

LR_formed = (Unbound(s_range_2)*R_Init)./(Unbound(s_range_2)+Kd_1_mod);
X4 = Bound(s_range_2) - LR_formed;
Y4 = (Bound(s_range_2)- LR_formed)./Unbound(s_range_2);

plot(X4,Y4,':yo',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','m',...
             'MarkerSize',5);
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
         
f = fittype('a*x + b');
[Scatchard_4,GOF_4] = fit(X4,Y4,f,'StartPoint',[-0.001 0.001]);
plot(Scatchard_4,'-m');
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

% Retrieve the fit parameters.
Scatchard_4
GOF_4
% or simply:
Scatchard_4_params = coeffvalues(Scatchard_4);
Kd_4 = -1/Scatchard_4_params(1);
Receptor_conc_4 = -Scatchard_4_params(2)/Scatchard_4_params(1);


legend('Scatchard','Fit\_1: points','Fit\_1: Kd\_1','Fit\_2: points','Fit\_2: Kd\_2','Corrected Scatchard',...
     'Corrected Fit\_2: Kd\_3','Corrected Fit\_3: Kd\_4')
string1 = 'Kd\_1 = ';
string2 = num2str(Kd_1,'%5.2e\n');
string1b = 'Kd\_1\_mod = ';
string2b = num2str(Kd_1_mod,'%5.2e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%5.2e\n');
string5 = 'Kd\_3 = ';
string6 = num2str(Kd_3,'%5.2e\n');
string0 = ' µM';
string7 = 'Kd\_4 = ';
string8 = num2str(Kd_4,'%5.2e\n');
string0 = ' µM';

xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');

% Create textbox
annotation(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.78 0.76 0.12 0.057],...
    'String',{[string1 string2 string0;...
    string3 string4 string0;...
    string5 string6 string0; string7 string8 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);
annotation(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.60 0.85 0.14 0.036],...
    'String',{[string1b string2b string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

title('Scatchard Plot');

%%
% Hyperbola fitting.

clear X Y
%

FR = Bound/(R_Init + R_1_Init);
X = Unbound(s_range);
Y = FR(s_range);

subplot3 = subplot(3,2,3,'Parent',figure(gcf));
box(subplot3,'on');
grid(subplot3,'on');
hold(subplot3,'all');

plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
         
% First we try to fit with just one component.
f = fittype('b*(x/(a + x))');
[Hyperb1,GOF1] = fit(X,Y,f,'StartPoint',[300 1]);
plot(X,Hyperb1(X),'-y');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

% Retrieve the fit parameters.
Hyperb1
GOF1
% or simply:
Hyperb1_params = coeffvalues(Hyperb1);
Kd_1 = Hyperb1_params(1);

% Then we use a two-component fit with 4 unknowns.         
f = fittype('b*(x/(a + x)) + d*(x/(c + x))');
[Hyperb2,GOF2] = fit(X,Y,f,'StartPoint',[Kd_1 0.5 Kd_2 0.5],...
    'Lower',[0.6*Kd_1 0.3 0.6*Kd_2 0.3],'Upper',[1.4*Kd_1 0.7 1.4*Kd_2 0.7]);
plot(X,Hyperb2(X),'-r');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

ylim([-0.05,1.05]); xlim([-50,Unbound(end)+50]);
legend('Rbound/Rtotal','1 Comp. Fit: Kd\_1','2 Comp. Fit: Kd\_2, Kd\_3','Location','SouthEast');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

% Retrieve the fit parameters.
Hyperb2
GOF2
% or simply:
Hyperb2_params = coeffvalues(Hyperb2);
Kd_2 = Hyperb2_params(1);
Kd_3 = Hyperb2_params(3);

string1 = 'Kd\_1 = ';
string2 = num2str(Kd_1,'%6.2e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%6.2e\n');
string5 = 'Kd\_3 = ';
string6 = num2str(Kd_3,'%6.2e\n');
string0 = ' µM';

% Create textbox
annotation(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.30 0.50 0.12 0.05],...
    'String',{[string1 string2 string0; string3 string4 string0;...
    string5 string6 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

title('Hyperbola fit ');

%%
% HILL plot
subplot4 = subplot(3,2,4,'Parent',figure(gcf));
box(subplot4,'on');
grid(subplot4,'on');
hold(subplot4,'all');

X = log(Unbound(s_range));
Y = log(FR(s_range)./(1 - FR(s_range)));
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
xlabel('log[Free Ligand]');
ylabel('log(Frac. Sat./(1-Frac. Sat.)');
title('Hill plot');
         
hold on
f = fittype('a*x + b');
[Hill_1,GOF_1] = fit(X(1:4),Y(1:4),f,'StartPoint',[1 -4]);
Hill_1
plot(Hill_1,'-b');
[Hill_2,GOF_2] = fit(X(25:35),Y(25:35),f,'StartPoint',[1 -5]);
Hill_2
plot(Hill_2,'-r');
[Hill_3,GOF_3] = fit(X(10:22),Y(10:22),f,'StartPoint',[0.6 -3]);
Hill_3
plot(Hill_3,'-g');
xlabel('log[Free Ligand]');
ylabel('log(Frac. Sat./(1-Frac. Sat.)');

ylim([-4 4]);
legend('Binding Curve','Hill Fit1: Kd\_1','Hill Fit2: Kd\_2','Hill Fit3: Kd\_3','Location','NorthWest');

% Retrieve the fit parameters.
Hill_1_params = coeffvalues(Hill_1);
Hill_2_params = coeffvalues(Hill_2);
Hill_3_params = coeffvalues(Hill_3);
nH_1 = Hill_1_params(1);
nH_2 = Hill_2_params(1);
nH_3 = Hill_3_params(1);
Kd_1 = exp(-Hill_1_params(2));
Kd_2 = exp(-Hill_2_params(2));
Kd_3 = exp(-Hill_3_params(2));

string1 = 'nH\_1 = ';
string2 = num2str(nH_1,'%5.2e\n');
string3 = 'nH\_2 = ';
string4 = num2str(nH_2,'%5.2e\n');
string5 = 'nH\_3 = ';
string6 = num2str(nH_3,'%5.2e\n');

% Create textbox
annotation(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.77 0.45 0.10 0.045],...
    'String',{[string1 string2 ; string3 string4 ; string5 string6]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

string1 = 'Kd\_1 = ';
string2 = num2str(Kd_1,'%6.2e\n');
string3 = 'Kd\_2 = ';
string4 = num2str(Kd_2,'%6.2e\n');
string5 = 'Kd\_3 = ';
string6 = num2str(Kd_3,'%6.2e\n');

% Create textbox
annotation(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.59 0.50 0.12 0.045],...
    'String',{[string1 string2 string0; string3 string4 string0; string5 string6 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);


xlabel('log[Free Ligand]');
ylabel('log(Frac. Sat./(1-Frac. Sat.)');
title('Hill plot');

%%
% Logistic fit.
%
clear X Y X1 X2 Y1 Y2 H H2

    subplot5 = subplot(3,2,5,'Parent',figure(gcf));
    box(subplot5,'on');
    grid(subplot5,'on');
    hold(subplot5,'all');

Y = Bound(s_range);
X = Unbound(s_range);

% We use a 4 parameter logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b)');
[Hill_1,GOF_1] = fit(X,Y,f,'StartPoint',[ 300 0.5 0.85 0.05 ],...
        'Lower',[0.6*300 0.1 0.6 0.0],'Upper',[1.4*300 0.9 1.1 0.1]);
Hill_1
GOF_1

plot(X,Hill_1(X),'-r');
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5,'Parent',subplot5);
ylim([-0.05 1.1]);
legend('Binding Curve','Logistic Fit','location','Best');

% Parameters    
hpar = coeffvalues(Hill_1);
string1 = 'Kd\_1 = ';
string2 = num2str(hpar(1),'%5.2e\n');

% Create textbox
annotation(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.3 0.15 0.12 0.03],...
    'String',{[string1 string2 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit');

%%
% Here we replot the fit as a Dose-Response curve.

    subplot6 = subplot(3,2,6,'Parent',figure(gcf),'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
    box(subplot6,'on');
    grid(subplot6,'on');
    hold(subplot6,'all');

% Experimental values
semilogx(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5,'Parent',subplot6);

ylim([-0.05 1.2]);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit as Dose-Response curve')
    
% Parameters    
hpar = coeffvalues(Hill_1);
Unbound_1 = Unbound(s_range);

% Theoretical values for nH = 2
for i = 1:length(Unbound_1)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound_1(i))^2);
end
semilogx(Unbound_1,H2,'-r','Parent',subplot6);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');

% Theoretical values for nH = 1
for i = 1:length(Unbound_1)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound_1(i))^1);
end
semilogx(Unbound_1,H2,'-g','Parent',subplot6);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');

% Theoretical values for nH = 0.5
for i = 1:length(Unbound_1)
H2(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound_1(i))^0.5);
end
semilogx(Unbound_1,H2,'-m','Parent',subplot6);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');

% Experimental values.
semilogx(Unbound_1,Hill_1(Unbound_1),'-b','Parent',subplot6);
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');


% Parameters    
hpar = coeffvalues(Hill_1);
string1 = 'nH_e_x_p = ';
string2 = num2str(hpar(2),'%4.2f\n');

legend('Dose/Response','nH = 2','nH = 1','nH = 0.5',[string1 string2],'location','Best');

string1 = 'Kd\_1 = ';
string2 = num2str(hpar(1),'%6.3e\n');

% Create textbox
% annotation(HILL_2,'textbox',...
%     [0.75 0.75 0.13 0.035],...
%     'String',{[string1 string2]});

xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Logistic fit as Dose-Response curve')

%%
% Thus, we can easily distinguish between a case of two independent binding
% sites and a case of negative cooperativity using either a Hill plot or a
% logistic fit. However, it is more difficult to get correct values for the
% low and high affinity Kd's.

%% Conclusion
%
% # The presence of positive cooperativity between binding sites can be
% inferred from the presence of a 'convex' Scatchard plot. In this case, an
% accurate determination of the Kd's of the individual binding sites may be
% difficult with standard methods like Schatchard or hyperbola fitting. It
% is often useful to carry out a logistic fit of the binding data and to
% compare its determination of the Hill coefficient with that derived from
% the Hill plot.
% # Negative cooperativity gives origin to a "concave" Scatchard, somewhat
% similar to that previously obtained with two independent binding sites,
% but more markedly represented by a single curve rather than by two
% distinct phases. None of the different analysis methods give very good
% values for the Kd's of the sites. Hill plot and logistic fit both show
% good performance in detecting the presence of negative cooperativity.

%%
close all
% save(savefile);