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

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

%% Path
close all
clear, clc

%% Declaration of global variables
global m1 conc_vec kon_R kon_R_1 koff_R koff_R_1 
global kon_R_B kon_R_1B koff_R_B koff_R_1B C R R_1 LR LR_1 LR_TOT 
global L C_Init R_Init R_1_Init

%%
% We start by loading the project 

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

% or alternatively:

% 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')
getequations(m1)

%% Cooperativity
% Notice how 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
%
%%
% 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_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

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

%% Positive cooperativity

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

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

%%
% 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.4 mseconds
Stop = 0.004;
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.

sbioplot(Binding_Kinetics)

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

%%
% 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);
% 
% sbioplot(Binding_Kinetics)

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

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

High_Conc = 600;
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;

% 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)
    error = (rand(3,1)-0.5);
    error(1) = 1 + error(1)*0.10;
    error(2) = 1 + error(2)*0.10;
    error(3) = 1 + error(3)*0.10;
    C.Capacity = C_Init*error(1);
    R.InitialAmount = R_Init*error(2)/C.Capacity;
    R_1.InitialAmount = R_1_Init*error(2)/C.Capacity;
    L.InitialAmount = conc_vec(i)*error(3)/C.Capacity;
    LR.InitialAmount = 0;    
    LR_1.InitialAmount = 0;    
    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

%%
Two_cooperative_binding_sites_1 = figure;
    set(Two_cooperative_binding_sites_1,'Units','normalized',...
        'Position',[0.1 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),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','c',...
             'MarkerSize',7.5);
plot(Unbound(range),Bound_R(range),':og',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','y',...
             'MarkerSize',7.5);
plot(Unbound(range),Bound_R_1(range),':ob',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','b',...
             'MarkerSize',2);

set(gca,'YLim',[-0.1 1.15]);
set(gca,'XLim',[-10 10+Unbound(end)]);

xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve ');

% 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. Since the sum of the two receptors is 1, the fractional
% saturation is equal to the 'bound' fraction.

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

f = fittype('(0.5*x/(a + x)) + (0.5*x/(b + x))');
[Hyperb,GOF] = fit(X,Y,f,'StartPoint',[3 30]);
 
plot(X,Hyperb(X),'-r');
% plot(X,Y,'o','MarkerEdgeColor','b','MarkerFaceColor','c','MarkerSize',7.5);
ylim([-0.1,1.15]); xlim([-10,Unbound(end)+10]);

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

legend('LR\_TOT','LR','LR\_1','Hyperb. fit','Location','Best');

%% Scatchard Plot 
    subplot2 = subplot(3,1,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(18:end) = false;
s_range_2 = true(1,35);
s_range_2(1:18) = false;

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

% Low 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');
[Scatchard_1,GOF_1] = fit(X1,Y1,f,'Robust','Bisquare','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);

% High 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);
hold on
f = fittype('a*x + b');
[Scatchard_2,GOF_2] = fit(X2,Y2,f,'Robust','Bisquare','StartPoint',[-1 1]);
% [Scatchard_2,GOF_2] = fit(X2,Y2,f,'StartPoint',[-1 1]);
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]);
xlim([0,1.05]);
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');

% Low affinity site (inferred)
Kd_1_corr = 1/(Scatchard_1_params(2)/Receptor_conc_2);

%% Hill Coefficient: Logistic fit
% 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.
%
% Hill coefficient: Hill plot
% This plot is derived from the Hill equation, in which
% the Hill coefficient 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.
%
%
% The traditional way of determining the Hill coefficient is with 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
%

    subplot3 = subplot(3,1,3,'Parent',figure(gcf));
    box(subplot3,'on');
    grid(subplot3,'on');
    hold(subplot3,'all');
    
FR = Bound/(R_Init + R_1_Init);
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(15:19),Y(15:19),f,'StartPoint',[2 -8],'Robust','on');
Hill_3
plot(Hill_3,'-g');

% ylim([-6 6]);
% xlim([0 10]);
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));

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

%%
% % An alternative way is with a direct 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.
% %
% % If there is only one component, we have 4 unknown: 
% % a = Kd, b = nH, c = max(y), d = min(y); 
% % y = (c - d)/(1 + (a/x)b) (1 components)
% % If there are two components we have a total of 8 unknown: 
% %  y = (c - d)/(1 + (a/x)b) + (g - h)/(1 + (e/x)f) (2 components)
% 
% clear X Y 
% HILL_4 = figure;
%     set(HILL_4,'Units','normalized','Position',[0.5 0.4 0.5 0.6],...
%         'Name','Binding Curve Fit');clf
%     subplot1 = subplot(2,1,1,'Parent',figure(gcf));
%     box(subplot1,'on');
%     grid(subplot1,'on');
%     hold(subplot1,'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_4,GOF_4] = 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_4,GOF_4] = fit(X,Y,f,'StartPoint',[ Kd_1 mean([nH_1 nH_2]) 0.5 0 Kd_2 0.5 0]);
% Hill_4
% GOF_4
% 
% 
% plot(X,Y,'o',...
%              'MarkerEdgeColor','b',...
%              'MarkerFaceColor','c',...
%              'MarkerSize',5);
% plot(Hill_4,'-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.
% 
%     subplot2 = subplot(2,1,2,'Parent',figure(gcf),'XScale','log',...
%         'XMinorTick','on','XMinorGrid','on');
%     box(subplot2,'on');
%     grid(subplot2,'on');
%     hold(subplot2,'all');
% 
% % 2 components (8parameters) fit:
% % Parameters    
% hpar = coeffvalues(Hill_4);
% 
% % Theoretical values for nH = 2
% for i = 1:length(Unbound)
% H4(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound(i))^2) + ...
%     (hpar(7) - hpar(8))/(1 + (hpar(5)./Unbound(i))^2);
% end
% semilogx(Unbound,H4,'-r','Parent',subplot2);
% 
% % Theoretical values for nH = 1
% for i = 1:length(Unbound)
% H4(i) = (hpar(3) - hpar(4))/(1 + (hpar(1)./Unbound(i))^1) + ...
%     (hpar(7) - hpar(8))/(1 + (hpar(5)./Unbound(i))^1);
% end
% semilogx(Unbound,H4,'-g','Parent',subplot2);
% 
% % Experimental values
% semilogx(Unbound,Hill_4(Unbound),'-b','Parent',subplot2);
% semilogx(X,Y,'o',...
%              'MarkerEdgeColor','b',...
%              'MarkerFaceColor','c',...
%              'MarkerSize',5,'Parent',subplot2);
% 
% ylim([-0.05 1.1]);
% 
% % Parameters    
% hpar = coeffvalues(Hill_4);
% 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');
% legend('nH = 2','nH = 1','location','Best');
% 
% xlabel('[Free Ligand]');
% ylabel('[Bound Ligand]');
% title('Logistic fit as Dose-Response curve')

%% Direct fit of progress curves: choice of parameters to refine.
close all

% We need to set some initial value for the refinement of parameters: we
% will refine koff_R, koff_RB.

% We take the Kd values from the Hill plot, and we multiply by the expected
% kon to obtain the values of koff.
p1 = Kd_1*10;
p2 = Kd_2*10;
pin=[p1;p2];

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

% Check vectors are the same dimensions:
% if size(Product_vector)~= size(Conc_vector)
%     Product_vector = Product_vector';
% end
%%
% Solution 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_cooperative_binding_sites_fit,pin,Conc_vector,Product_vector,...
[],[],options);

nobs = length(Conc_vector);
nvars = length(u);

MSE = sos/(nobs-nvars)
RMSE = sqrt(MSE)
J = full(J);
Cov = inv(J'*J).*MSE;
[Corr,sigma] = corrcov(Cov);
uci = [u-1.96*sigma u+1.96*sigma] 
u

% Here we use the nlinfit function from the Statistics Toolbox. This is the
% best choice if the toolbox is available. Can be used after lsqcurvefit to
% first narrow the radius of convergence.

% options=statset('TolX',1e-10,'TolFun',1e-10,'RobustWgtFun','bisquare',...
%     'Display','iter','MaxIter',500);
% if exist('u','var')
%     pin = u;
% end
% options=statset('TolX',1e-9,'TolFun',1e-9,'Display','iter',...
%     'MaxIter',200);
% [u,res1,J1,Cov1,MSE1] = ...
%     nlinfit(Conc_vector,Product_vector,'two_cooperative_binding_sites_fit',pin,options);
% [Corr1,sigma1] = corrcov(Cov1);
% sos1 = res1'*res1
% [fcurve1,delta1] = nlpredci('two_cooperative_binding_sites_fit',Conc_vector,u,R1,'Covar',Cov1);
% u_ci = nlparci(u,res1,'covar',Cov1);

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

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

for i = 1:length(conc_vec)
    C.Capacity = C_Init;
    R.InitialAmount = R_Init;
    R_1.InitialAmount = R_1_Init;
    L.InitialAmount = conc_vec(i);
    LR.InitialAmount = 0;    
    LR_1.InitialAmount = 0;    
    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

% 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','--c',conc_vec,Bound_R_1','--m'); 

%%
% We can compare the residual from this fit with that we would obtain taking
% the Hill plot solution:
koff_R.Value = Kd_1*10;
koff_R_1.Value = Kd_1*10;
koff_R_B.Value = Kd_2*10;
koff_R_1B.Value = Kd_2*10;

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

for i = 1:length(conc_vec)
    C.Capacity = C_Init;
    R.InitialAmount = R_Init;
    R_1.InitialAmount = R_1_Init;
    L.InitialAmount = conc_vec(i);
    LR.InitialAmount = 0;    
    LR_1.InitialAmount = 0;    
    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

R_hill = (Bound-Product_vector);
sos_hill = R_hill'*R_hill
sos

plot(conc_vec,Bound,'--c');hold off
legend('Direct Fit','Observed','Hill plot fit')
xlabel('[L] microM  ')
ylabel('Fractional saturation  ')

%%
% close all
save(savefile);

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