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

%% Very Simple Binding of a Ligand to a Receptor

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

close all
clear, clc

sbioloadproject('../TOOLBOXES/BINDING_KINETICS/Simple_Binding');
 
% or alternatively:
 
% m1 = sbmlimport('TOOLBOXES/BINDING_KINETICS/Simple_Binding.xml');

% We get information about the model:

sbioselect(m1,'Type','compartment')
sbioselect(m1,'Type','species')
sbioselect(m1,'Type','parameter')
sbioselect(m1,'Type','reaction')
getequations(m1)

% It is convenient to 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');

% Before we start the simulation we can change all the necessary parameters
% and variables to values different from those originally stored in the
% project. For convenience of representation we will use values of kon and
% koff so that Kd = koff/kon = 10 µM:

kon.Value = 10;
koff.Value = 100;
R.InitialAmount = 1;
L.InitialAmount = 200;
 
Kd = koff.Value/kon.Value;
display(['Kd = ' num2str(Kd) ' µM']);

% We can also obtain the stored parameters for the simulation:

configset_m1 = getconfigset(m1);
get(configset_m1)
configset_m1
 
% and change them to something else. For example, we will use 5 ms as the
% stop time and will set the solver to ODE23t (trapezoidal solver,
% moderatively stiff) or ODE15s (RK solver, stiff).
 
set(configset_m1, 'StopTime', 0.005);
set(configset_m1, 'SolverType', 'ode15s');
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

% Now we can run a small simulation with the existing parameters. Only 5
% 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
% array.

get(Binding_Kinetics)
openvar Binding_Kinetics.Time
openvar Binding_Kinetics.Data
 
% sbioplot(Binding_Kinetics)
 
Simple_binding = figure;
set(Simple_binding,'Units','normalized',...
'Position',[0.6 0.6 0.4 0.4],'Name','Simple Binding');
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');

% Now we are ready to simulate an entire binding curve. In order to examine
% the effect of different concentrations of the ligand 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.

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

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 equilibrium.
 
Stop = 0.5;
set(configset_m1, 'StopTime', Stop);

% We store in different variables the initial amounts of each species; this
% becomes very important later on when we loop through different
% concentrations of the ligand. 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.
 
R_Init = 1;
LR_Init = 0;

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

%%
Binding_analysis = figure;
set(Binding_analysis,'Units','normalized','Position',[0.4 0.0 0.6 1.0],...
    'Name','Binding Curve ');

% Hyperbola fit
FR = Bound/R_Init;          % Stored initial value!
X = Unbound;
Y = FR;
 
subplot1 = subplot(2,1,1,'Parent',figure(gcf));
box(subplot1,'on');
grid(subplot1,'on');
hold(subplot1,'all');
 
f = fittype('x/(a + x)');
[Hyperb,GOF] = fit(X,Y,f,'StartPoint',10);
 
plot(X,Hyperb(X),'-r');
ylim([0,1.1]);
xlim([-5 215]);
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
 
legend('Fit','Rbound/Rtotal','Location','Best');
xlabel('[Free Ligand]');
ylabel('Fractional saturation');
title('Hyperbola fit ');

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

% Scatchard fit    
subplot2 = subplot(2,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);

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

plot(Scatchard,'-r');
plot(X,Y,'s',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','g',...
             'MarkerSize',5);
legend('Fit','Scatchard');
xlabel('[Bound Ligand]');
ylabel('[Bound Ligand]/[Free Ligand]');
title('Scatchard fit ');

% 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(Binding_analysis,'textbox',...
    [0.8 0.85 0.09 0.025],...
    'String',{[string1 string2 string3]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.8],...
    'Color',[1 0 0]);
