% Copyright (c) 2015, 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 Independent Binding Sites of the Same Receptor 

global m1 conc_vec kon kon_R_1 koff koff_R_1 scale C R R_1 LR LR_1 LR_TOT 
global scale_R scale_R1 L C_Init

%%
% We start by loading the project 

% sbioloadproject('Two_Independent_Binding_Sites');

% or alternatively:

m1 = sbmlimport('../TOOLBOXES/BINDING_KINETICS/Two_Independent_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_Independent_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
% the by starting SimBiology desktop with |simbiology| and opening the
% project file Two_Independent_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)

%%
% 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 = 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');

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

% Before we start the simulation we set all the necessary parameters
% and variable values.

C.Capacity = 1;
L.InitialAmount = 3;
R.InitialAmount = 5;
R_1.InitialAmount = 5;
LR.InitialAmount = 0;
LR_1.InitialAmount = 0;
LR_TOT.InitialAmount = 0;
scale = 1;
scale_R = 1;
scale_R1 = 1;

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

% Here we store the initial parameters.
C_Init = 1;
R_Init = 5;
R_1_Init = 5;
LR_Init = 0;
LR_1_Init = 0;
LR_TOT_Init = 0;

kon.Value = 10;
koff.Value = 170; % Kd(R) = 17 µM
kon_R_1.Value = 10;
koff_R_1.Value = 3500; % Kd(R_1) = 350 µM

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

configset_m1 = getconfigset(m1);

% We will set the initial stop time at 50 mseconds
Stop = 0.05;
set(configset_m1, 'StopTime', Stop);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

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

% 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:5 )*scale)
legend('Free R','Bound R','Bound R\_1','Total Bound R','Location','Best')
set(gca,'YLim',[-0.05 6]);
set(gca,'XLim',[-0.001 1.05*Stop]);
ylabel('Species Concentration (µM)')
xlabel('Time (s)')
title('Progress Curve with [L] = 3 microM');

%%
% We repeat considering the other possible extreme in the concentration of
% the ligand.

% close all

C.Capacity = 1;
L.InitialAmount = 1000;
R.InitialAmount = 0.5;
R_1.InitialAmount = 0.5;
LR.InitialAmount = 0;
LR_1.InitialAmount = 0;
LR_TOT.InitialAmount = 0;

kon.Value = 10;
koff.Value = 170; % Kd(R) = 17 µM
kon_R_1.Value = 10;
koff_R_1.Value = 3500; % Kd(R_1) = 350 µM

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

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

% We recall that Kd(R) = 17 µM and Kd(R_1) = 350 µ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):

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

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;

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

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

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.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)*scale;
    Bound_R(i) = Binding_Kinetics.Data(end,3)*scale_R;
    Bound_R_1(i) = Binding_Kinetics.Data(end,4)*scale_R1;
    Unbound(i) = Binding_Kinetics.Data(end,1);
end

Binding_Kinetics.DataNames

plot(conc_vec,Bound,'or',conc_vec,Bound,'-b'); hold on
plot(conc_vec,Bound_R,'--m',conc_vec,Bound_R_1,'--c'); hold off


%% Choice of parameters to refine.
% We need to set some initial value for the refinement of parameters: we
% will refine koff, koff_R1, scale_R, scale_R1.

p1 = 200;
p2 = 2000;
p3 = 1;
p4 = 1;
pin=[p1;p2;p3;p4];

%% Non-linear global fit of all the progress curves.
% Here we set the vector of experimental values against which we refine our
% simulated product vector:
Product_vector = Bound;
Conc_vector = conc_vec;
 
% Here we use the nlinfit function from the Statistics Toolbox. This is the
% best choice if the toolbox is available.
% options=statset('TolX',1e-10,'TolFun',1e-10,'RobustWgtFun','bisquare',...
%     'Display','iter','MaxIter',500);
options=statset('TolX',1e-10,'TolFun',1e-10,'Display','iter',...
    'MaxIter',200);
[u,R1,J1,Cov1,MSE1] = ...
    nlinfit(Conc_vector,Product_vector,'two_independent_binding_sites_fit',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
[fcurve1,delta1] = nlpredci('two_independent_binding_sites_fit',Conc_vector,u,R1,'Covar',Cov1);
u_ci = nlparci(u,R1,'covar',Cov1);

%%
% Alternative solution with lsqcurvefit: this function usually has a wider 
% radius of convergence, but is slower.
options = ...
optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-8,'TolX',1e-8);
[u,sos,res,~,~,~,J] = ...
lsqcurvefit(@two_independent_binding_sites_fit,pin,Conc_vector,Product_vector,...
[100;1000;0.5;0.5],[300;3000;1.5;1.5],options); 

%%
% Alternative solution with the generic Matlab function fminsearch (uses
% Simplex search rather than derivatives: it minimizes the scalar output of
% a function. For this reason the model function includes the calculation
% of the sum of squared residuals, and therefore there is no Jacobian).
% This is the slowest method, but in some cases it may be able to find the
% global minimum better than other methods.

modelfun = @(pin) sum((two_independent_binding_sites_fit(pin,Conc_vector)-Product_vector).^2);
options = optimset('Display','iter', 'TolFun',1e-6, 'TolX',1e-6);
[u,fval,exitflag,output] = fminsearch(modelfun,pin,options);

%% Regenerate binding curves
koff.Value = u(1);
koff_R_1.Value = u(2);
scale_R = u(3);
scale_R1 = u(4);
Bound_R = zeros(length(conc_vec),1);
Bound_R_1 = zeros(length(conc_vec),1);

for i = 1:length(conc_vec)
    set(L,'InitialAmount',conc_vec(i));
    Binding_Kinetics = sbiosimulate(m1);
    Bound_R(i) = Binding_Kinetics.Data(end,3);
    Bound_R_1(i) = Binding_Kinetics.Data(end,4);
    
end
Bound = Bound_R*scale_R + Bound_R_1*scale_R1;    

% figure;plot(conc_vec,fcurve1,'-b',conc_vec,Product_vector,'or');hold on
figure;plot(conc_vec,Bound,'-b',conc_vec,Product_vector,'or');hold on
plot(conc_vec,Bound_R'*scale_R,'--c',conc_vec,Bound_R_1'*scale_R1,'--m'); hold off


%%
close all
% save(savefile);