% Copyright (c) 2013, 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.
% 

%% Simple Binding of a Ligand to a Receptor

% This example shows how to analyze a simple binding experiment.

% We start by loading the project 

close all
clear

sbioloadproject('../TOOLBOXES/BINDING_KINETICS/Simple_Binding');

% or alternatively:

% m1 = sbmlimport('Simple_Binding.xml');

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

%%
% The project contains a model, |m1|, of the simple binding of a ligand to
% a receptor. We can also load the model interactively by starting the
% SimBiology desktop with |simbiology| and opening the project file
% Simple_Binding.sbproj.

%%
% <<../Simple_Binding.jpg>>

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

%%
% Now we want to see what happens when we change some of the parameters.
% First, 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');
LR = sbioselect(m1,'Name','LR');
kon = sbioselect(m1,'Name','kon');
koff = sbioselect(m1,'Name','koff');

% Check the variable values
% C
% L
% R
% LR
% kon
% koff

%% Time Evolution of the system.
%
% Before we start the simulation we will set all the necessary parameters
% and variable values. For convenience of representation we will use the
% following values of kon and koff so that Kd = koff/kon = 10 µM and
% the concentration of the receptor is [R] = 1 µM:

kon.Value = 10;
koff.Value = 100;
R.InitialAmount = 1;
L.InitialAmount = 200;

Kd = koff.Value/kon.Value;
display(['Kd = ' num2str(Kd) ' µM']);

%%
% Now we will run a small simulation with the existing parameters. Let's
% first see what the stored parameters are:

configset_m1 = getconfigset(m1);
get(configset_m1)

% or

configset_m1

% We will use 5 ms as the stop time and will set a high value for the
% sensitivity of the solver.

set(configset_m1, 'StopTime', 0.005);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);
set(configset_m1, 'SolverType', 'ode23t');
% We can make sure the configuration has actually changed:
% configset_m1

% We run the simulation with the new StopTime.

Binding_Kinetics = sbiosimulate(m1);

%%
% We can get some information about the simulation.
% get(Binding_Kinetics)

% 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 array.

% openvar Binding_Kinetics.Time
% openvar Binding_Kinetics.Data

% 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 be more selective and plot only the 
% concentrations of the free and bound receptor.

% close all

%%
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:3))
legend('Free R','Bound R','Location','Best')
set(gca,'YLim',[-0.005 1.05]);
set(gca,'XLim',[-0.0002 0.0052]);
ylabel('Species Concentration (µM)')
xlabel('Time (s)')
title('Progress Curve with [L] = 200 µM');

%%
% We will repeat the simulation using a ligand concentration below the Kd.
% For this simulation we will take a longer stop time (50 msec).
close all
set(configset_m1, 'StopTime', 0.05);

L.InitialAmount = 1;
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:3))
legend('Free R','Bound R','Location','Best')
set(gca,'YLim',[-0.005 1.05]);
set(gca,'XLim',[-0.001 0.05]);
ylabel('Species Concentration (µM)')
xlabel('Time (s)')
title('Progress Curve with [L] = 1 µM');

%% Simulation of a binding curve.

close all
clear Bound Unbound

% 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 in order
% to completely saturate the receptor. Thus, we will cover the range of
% concentrations 1-200 µM for which we already know the amount of time
% we need to follow the reaction in order to achieve equilibrium.

High_Conc = 200;
conc_vec = logspace(log10(1),log10(High_Conc),35);
% openvar conc_vec

range = conc_vec <= 1 + High_Conc;
s_range = range; % This is the range we will use for the Scatchard plot
Bound = zeros(length(conc_vec),1);
Unbound = zeros(length(conc_vec),1);
time = Binding_Kinetics.Time;

% 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 the binding equilibrium.

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

% It is also very important that we store in different variables the
% initial amounts of each species; this will be very important later on
% when we loop through different concentrations of the ligand.

C_Init = 1;
R_Init = 1;
LR_Init = 0;

for i = 1:length(conc_vec)
    set(R,'InitialAmount',1);    
    set(L,'InitialAmount',conc_vec(i));
    set(LR,'InitialAmount',0);  % not necessary, but OK to be sure is 0
    Binding_Kinetics = sbiosimulate(m1);
    Bound(i) = Binding_Kinetics.Data(end,3);       % 3rd data column
    Unbound(i) = Binding_Kinetics.Data(end,1);     % 1st data column 
end

%% Binding Curve with error-free data
% Here we plot all the results obtained with theoretical data.
% 
FIG_3 = figure;
    set(FIG_3,'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');

% Plotting versus the free [L]
plot(Unbound(range),Bound(range),'--or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5);
% plot(conc_vec(range),Unbound(range),'--ob',...
%              'MarkerEdgeColor','k',...
%              'MarkerFaceColor','c',...
%              'MarkerSize',5);
set(gca,'YLim',[-0.05 1.1]);
set(gca,'XLim',[-5 205]);
         
legend('Bound Ligand','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve with error-free data');

%% Scatchard Plot
% Now we can choose various ways of analyzing the data. We start with a
% Scatchard plot: this is a plot of the ratio of concentrations of bound
% ligand to unbound ligand versus the bound ligand concentration:
%
% Scatchard Plot:
%
% Y = [Bound]/[Unbound]
%
% X = [Bound] 
% 
% Slope = -1/Kd.
%
% X axis intercept = [R]
%
% Sometimes the binding data do not form a straight line when plotted in a
% Scatchard plot. One case in which this can happen is when ligand bound to
% substrate is not allowed to achieve equilibrium before the binding is
% measured (see below). Another case is when there is more than one binding
% site on the same receptor (see example "Two_Independent_Binding_Sites").
%
    subplot2 = subplot(3,2,2,'Parent',figure(gcf));
    box(subplot2,'on');
    grid(subplot2,'on');
    hold(subplot2,'all');

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

% First we fit:
f = fittype('a*x + b');
[Scatchard,GOF] = fit(X,Y,f,'StartPoint',[-10 50]);

% Then we plot the fit and the data:
plot(Scatchard,'-r');
plot(X,Y,'s',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);
legend('Fit','Scatchard');
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
% ylim([0,Scatchard(0)]);

% Retrieve the fit parameters.
Scatchard_params = coeffvalues(Scatchard);
Kd = -1/Scatchard_params(1);
Receptor_conc = -Scatchard_params(2)/Scatchard_params(1);
string1 = 'Kd = ';
string2 = num2str(Kd,'%6.2f\n');
string3 = ' µM';
% Create textbox
annotation(FIG_3,'textbox',...
    [0.8 0.85 0.09 0.025],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);
title('Scatchard Plot with error-free data');

%% Data with noise
% Of course this is what we would obtain with perfect data. But the
% situation could be much worse if there is noise in the data.

% close all

% for i = 1:length(conc_vec)
%     error = (rand(3,1)-0.5);
%     error(1) = 1 + error(1)*0.10;
%     error(2) = 1 + error(2)*0.05;
%     error(3) = 1 + error(3)*0.05;
%     set(C,'Capacity',C_Init*error(1));          % Stored initial value!
%     set(R,'InitialAmount',R_Init*error(2));     % Stored initial value!
%     set(L,'InitialAmount',conc_vec(i)*error(3));
%     set(LR,'InitialAmount',LR_Init);            % Stored initial value!
%     Binding_Kinetics = sbiosimulate(m1);
%     Bound(i) = Binding_Kinetics.Data(end,3);
%     Unbound(i) = Binding_Kinetics.Data(end,1);
% end

for i = 1:length(conc_vec)
    error = (rand(2,1)-0.5);
    error(1) = 1 + error(1)*0.15;
    error(2) = 1 + error(2)*0.15;
    set(R,'InitialAmount',R_Init*error(1));     % Stored initial value!
    set(L,'InitialAmount',conc_vec(i)*error(2));
    set(LR,'InitialAmount',LR_Init);            % Stored initial value!
    Binding_Kinetics = sbiosimulate(m1);
    Bound(i) = Binding_Kinetics.Data(end,3);
    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);

set(gca,'YLim',[-0.05 1.1]);
set(gca,'XLim',[-5 205]);
legend('Bound Ligand','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve with random-error data');

%% Dose-Response Curve from data with noise
% A Dose-Response curve is a plot of the response (the bound ligand
% concentration), or the percentage of the maximal response, against the
% log of the dose (the free ligand concentration).

    subplot4 = subplot(3,2,4,'Parent',figure(gcf),'XScale','log',...
        'XMinorTick','on','XMinorGrid','on');
    box(subplot4,'on');
    grid(subplot4,'on');
    hold(subplot4,'all');
    
    semilogx(Unbound(range),Bound(range),'--or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',subplot4);

    set(gca,'YLim',[-0.05 1]);
    set(gca,'XLim',[1 210]);
    vline([9.8 9.8],{'g-','LineWidth',2});

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

% We will use Dose-Response curves extensively in our examples. They are
% particularly useful because at a glance we can see if we have saturated
% the receptor (the part of the curve corresponding to the highest
% concentrations of ligand should be approximately flat), and we can get
% the Kd from the middle of the sigmoid curve.

%% Scatchard Plot from data with noise

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

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

% First we fit the data by linear regression (least square fit) to the
% equation of a straight line.
f = fittype('a*x + b');
[Scatchard,GOF] = fit(X,Y,f,'StartPoint',[-10 50]);

% Then we plot the fit and the data:
plot(Scatchard,'-r');
plot(X,Y,'s',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);

% ylim([0,Scatchard(0)]);
legend('Fit','Scatchard');
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
title('Scatchard plot with random-error data');

% Retrieve the fit parameters.
Scatchard
GOF
% or simply:
Scatchard_params = coeffvalues(Scatchard);
Kd = -1/Scatchard_params(1);
Receptor_conc = -Scatchard_params(2)/Scatchard_params(1);

string1 = 'Kd = ';
string2 = num2str(Kd,'%6.2f\n');
string3 = ' µM';

% string2 = num2str(Kd,5);
% Create textbox
annotation(FIG_3,'textbox',...
    [0.36 0.25 0.09 0.025],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

%% Hyperbola fitting
%
% However, rather than using a Scatchard plot we could fit directly the
% plot of the fractional saturation of the receptor (FR) vs the amount of 
% free ligand (L). For this purpose we will use the hyperbolic function
% we presented at the beginning of the example:
%
% [LR] = [R0]*[L]/([L]+Kd)
%
% which allows to determine the concentration of the complex at equilibrium
% provided the Kd, the total concentration of receptor [R0], and the
% concentration of the free ligand [L] are known. We rearrange the function
% in terms of fractional saturation:
%
% [LR]/[R0] = [L]/([L]+Kd)
%

FR = Bound/R_Init;          % Stored initial value!
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');

% As usual, first we fit:         
f = fittype('x/(a + x)');
[Hyperb,GOF] = fit(X,Y,f,'StartPoint',10);

% and then we plot everything:
plot(X,Hyperb(X),'-r');
ylim([0,1.1]);
xlim([-5 205]);
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);

legend('Fit','Rbound/Rtotal','Location','Best');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

% Retrieve the fit parameters.
Hyperb
GOF
% or simply:
Hyperb_params = coeffvalues(Hyperb);
Kd = Hyperb_params(1);

string1 = 'Kd = ';
string2 = num2str(Kd,'%6.2f\n');
string3 = ' µM';
% Create textbox
annotation(FIG_3,'textbox',...
    [0.8 0.25 0.087 0.025],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);
title('Hyperbola Fit with random-error data');

%% Effect of different kon and koff at constant Kd.
% We may wonder what happens if we change the rate constants in our model
% keeping the ratio between them (the Kd) constant.

close all

kon.Value = 0.0001;
koff.Value = 0.001;
R.InitialAmount = 1;
L.InitialAmount = 200;

set(configset_m1, 'StopTime', 500); 
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

Binding_Kinetics = sbiosimulate(m1);

FIG_4 = figure;
    set(FIG_4,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','Simple Binding');clf
    
axes1 = axes('Parent',FIG_4,...
    'Position',[0.08 0.1 0.9 0.82]);
hold(axes1,'all');
box(axes1,'on');
plot(Binding_Kinetics.Time,Binding_Kinetics.Data(:,2:3))
legend('Free R','Bound R','Location','Best')
set(gca,'YLim',[-0.005 1.05]);
set(gca,'XLim',[-0.001 500.5]);
ylabel('Species Concentration (µM)')
xlabel('Time (s)')
title('Progress Curve with [L] = 200 µM');

% Notice that the curves look identical to the ones we obtained with faster
% rate constants, but the time required to reach equilibrium is now ~ 500
% seconds !!

%%
% Likewise we can repeat the simulation with a small amount of ligand.

kon.Value = 0.0001;
koff.Value = 0.001;
R.InitialAmount = 1;
L.InitialAmount = 1;

set(configset_m1, 'StopTime', 5000);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

Binding_Kinetics = sbiosimulate(m1);

FIG_5 = figure;
    set(FIG_5,'Units','normalized','Position',[0.6 0.6 0.4 0.4],...
        'Name','Simple Binding');clf
    
axes1 = axes('Parent',FIG_5,...
    'Position',[0.08 0.1 0.9 0.82]);
hold(axes1,'all');
box(axes1,'on');
plot(Binding_Kinetics.Time,Binding_Kinetics.Data(:,2:3))
legend('Free R','Bound R','Location','Best')
set(gca,'YLim',[-0.005 1.05]);
set(gca,'XLim',[-0.001 5000.5]);
ylabel('Species Concentration (µM)')
xlabel('Time (s)')
title('Progress Curve with [L] = 1 µM');

% The time required to reach equilibrium is 5000 seconds !!

%% Failure to reach equilibrium
%
% FIG_4 and FIG_5 look exactly identical to FIG_1 and FIG_2, the only
% difference being that the time scale is 100,000 times bigger. The ratio
% between the two rate constants, and therefore the Kd, has not changed:
% from the experimental point of view this means only that we would have to
% wait up to 2 hours for each point to make sure that a chemical
% equilibrium has been reached. 
%
% We may wonder how big our error in the Kd determination would be if we
% waited only 2 minutes per point. In order to answer this question we
% repeat everything as before changing only the stop time.

set(configset_m1, 'StopTime', 120);

for i = 1:length(conc_vec)
    set(R,'InitialAmount',1);    
    set(L,'InitialAmount',conc_vec(i));
    set(LR,'InitialAmount',0);  % not necessary, but OK to be sure is 0
    Binding_Kinetics = sbiosimulate(m1);
    Bound(i) = Binding_Kinetics.Data(end,3);       % 3rd data column
    Unbound(i) = Binding_Kinetics.Data(end,1);     % 1st data column 
end

%%
% We will now repeat some of the analyses we carried out in the first part
% of the example. To simplify things we use only error-free data.
%

s_range = range;
s_range_1 = true(1,35);
s_range_1(10:end) = false;
s_range_2 = true(1,35);
s_range_2(1:31) = false;

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

% Plotting versus the free [L]
plot(Unbound(range),Bound(range),'--or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5);

set(gca,'YLim',[-0.05 1.0]);
set(gca,'XLim',[-5 205]);
         
legend('Bound Ligand','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve');

% Now we can choose various ways of analyzing the data. We start with a
% Scatchard plot of Bound/Unbound vs Bound. With this plot slope = -1/Kd and 
% the intercept on the X axis is the concentration of the receptor.

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

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

% First we fit the data corresponding to the entire range.

f = fittype('a*x + b');
[Scatchard_1,GOF_1] = fit(X1,Y1,f,'StartPoint',[-10 50]);
plot(Scatchard_1,'-c');

%ylim([Scatchard_1(1),Scatchard_1(0)]);
xlim([0 1]);
% Retrieve the fit parameters.
Scatchard_1
GOF_1
Scatchard_1_params = coeffvalues(Scatchard_1);
Kd1 = -1/Scatchard_1_params(1);
Receptor_conc_1 = -Scatchard_1_params(2)/Scatchard_1_params(1);
string1 = 'Kd\_1 = ';
string2 = num2str(Kd1,'%5.2e\n');
string3 = ' µM';

% Then we fit the data corresponding to the left hand side part of the plot.

X2 = Bound(s_range_1);
Y2 = Bound(s_range_1)./Unbound(s_range_1);

f = fittype('a*x + b');
[Scatchard_2,GOF_2] = fit(X2,Y2,f,'StartPoint',[-0.01 0.05]);
plot(Scatchard_2,'-b');
% ylim([0,Scatchard_1(0)]);
ylim([0.004,Scatchard_2(0)]);
% 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);
string4 = 'Kd\_2 = ';
string5 = num2str(Kd_2,'%5.2e\n');

% And the data corresponding to the right hand side part of the plot.

X3 = Bound(s_range_2);
Y3 = Bound(s_range_2)./Unbound(s_range_2);

f = fittype('a*x + b');
[Scatchard_3,GOF_3] = fit(X3,Y3,f,'StartPoint',[-0.003 0.003]);
plot(Scatchard_3,'-r');
xlim([0,1]);

% Retrieve the fit parameters.
Scatchard_3
GOF_3
Scatchard_3_params = coeffvalues(Scatchard_3);
Kd_3 = -1/Scatchard_3_params(1);
Receptor_conc_3 = -Scatchard_3_params(2)/Scatchard_3_params(1);
string6 = 'Kd\_3 = ';
string7 = num2str(Kd_3,'%5.2e\n');
 
% Finally we plot also the Scatchard data points:
plot(X1,Y1,'-ys',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);

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

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

legend('Fit\_1','Fit\_2','Fit\_3','Scatchard');
title('Scatchard Plot');

% We can also try an hyperbola fitting.

FR = Bound/R_Init;          % Stored initial value!
X = Unbound(s_range);
Y = FR(s_range);

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

% As usual, first we fit:         
f = fittype('b*x/(a + x)');
[Hyperb0,GOF0] = fit(X,Y,f,'StartPoint',[10 1]);
Hyperb0

% Notice how b is larger than 1, which makes no chemical sense: this is
% also a strong indication that the data points were not collected at
% equilibrium.

f = fittype('x/(a + x)');
[Hyperb,GOF] = fit(X,Y,f,'StartPoint',[10]);

% and then we plot everything:
plot(Unbound,Hyperb0(Unbound),'-y');
plot(Unbound,Hyperb(Unbound),'-r');
ylim([0,1.1]);
xlim([-5 205]);
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);

legend('Fit with scale > 1','Fit with scale = 1','Rbound/Rtotal','Location','Best');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');

% Retrieve the fit parameters.
Hyperb
GOF
% or simply:
Hyperb_params = coeffvalues(Hyperb);
Kd = Hyperb_params(1);

string1 = 'Kd = ';
string2 = num2str(Kd,'%6.2f\n');
string3 = ' µM';
% Create textbox
annotation(FIG_6,'textbox',...
    [0.16 0.35 0.087 0.025],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);
title('Hyperbola Fit');

% 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');
    
    semilogx(Unbound(range),Bound(range),'--or',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',5,'Parent',subplot4);

    set(gca,'YLim',[-0.05 1]);
    set(gca,'XLim',[1 210]);
    vline([60 60],{'g-','LineWidth',2});

legend('Bound Ligand','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Dose-Response Curve');

%% Conclusion
%
% A 'curved down' or 'convex' Scatchard plot is very often an indication
% that the experimental points were not collected at equilibrium. In this
% case neither the 'global' fit of the Scatchard plot (Fit_1) nor a more
% selective fit of the left hand side (Fit_2) or the right hand side
% (Fit_3) of the Scatchard plot yield a correct Kd value.
%
% A hyperbola fitting and a dose-response curve also fail to provide the
% correct Kd value. If the hyperbola fit is carried out using a scale
% factor, it catches well all the points, but the scaling factor refines to
% a value larger than 1. This is also indication that the points were not
% collected at equilibrium.
%
% Another possible explanation for a convex Scatchard is that there is more
% than one binding site, and there is cooperativity between the sites (see
% example 'Two Cooperative Binding Sites'). In the example 'Two Independent
% Binding Sites' we will also see a case in which the Scatchard curve is
% 'curved up' or concave.
%
%%
close all
% save(savefile);
