% 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.
%
%% Build negative cooperativity binding curve
close all
clear, clc

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

savefile = 'Problem_Binding_Curve.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;
R_Total_Init = R_Init + R_1_Init;

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

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

%% Skip this section if starting from external file 
% 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)
    error = (rand(3,1)-0.5);
    error(1) = 1 + error(1)*0.1;
    error(2) = 1 + error(2)*0.06;
    error(3) = 1 + error(3)*0.06;
    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

time = Binding_Kinetics.Time;

%% External data
close all
clear Bound Unbound Bound_R Bound_R_1

data = dlmread('../DATABASE/Problem_Binding_Curve.txt');
conc_vec = data(:,1)';
Bound = data(:,2);

%% Logistic fit.

% Reset Unbound to original concentration vector

Y = Bound;
X = conc_vec';
max_y = Bound(end);
min_y = Bound(1);
diff = abs(Bound - repmat(max_y/2,35,1))
[~,Kd_ind] = min(diff)
Kd = X(Kd_ind);

% We use a 4 parameter logistic fit:
         
f = fittype('(c - d)/(1 + (a/x)^b)');
[Log_1,GOF_1] = fit(X,Y,f,'StartPoint',[ Kd 1 max_y min_y ],...
        'Lower',[0.6*Kd 0.2 0.6*max_y 0.6*min_y],'Upper',[1.4*Kd 2 1.4*max_y 1.4*min_y ]);
Log_1
GOF_1

LOGISTIC_FIT = figure
plot(X,Log_1(X),'-r');hold on
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
ylim([-0.05 1.1]);
legend('Binding Curve','Logistic Fit','location','Best');

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

% Create textbox
annotation(LOGISTIC_FIT,'textbox',...
    [0.5 0.35 0.25 0.06],...
    'String',{[string1 string2 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

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

%% Rescaling of Bound and Unbound
% Bound becomes the fractional saturation
Bound = Bound/hpar(3);
Unbound = conc_vec' - Bound*R_Total_Init;
range = true(length(Bound),1);

%%
% Here we plot all the results.
%
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);

set(gca,'YLim',[-0.05 1.1]);
set(gca,'XLim',[-100 100+conc_vec(end)]);
         
legend('LR\_TOT','Location','Best');
xlabel('[Free Ligand]');
ylabel('[Bound Ligand]');
title('Binding Curve');

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

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

legend('Scatchard','Fit\_1: points','Fit\_1: Kd\_1','Fit\_2: points','Fit\_2: Kd\_2')
string1 =  'Kd\_1\_orig = ';
string2 = num2str(Kd_1_orig,'%5.2e\n');
string1b = 'Kd\_1       = ';
string2b = 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(NEG_COOPERATIVE_BINDING,'textbox',...
    [0.78 0.76 0.16 0.061],...
    'String',{[string1 string2 string0;...
    string1b string2b string0; string3 string4 string0]},...
    'FontWeight','bold',...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 0.800000011920929],...
    'Color',[1 0 0]);

title('Scatchard Plot');

Kd_1_scatchard = Kd_1;
Kd_2_scatchard = Kd_2;
%%
% Hyperbola fitting.

clear X Y
%

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

Kd_1_hyperb = Kd_2;
Kd_2_hyperb = Kd_3;

%%
% 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(3:13),Y(3:13),f,'StartPoint',[1 -4]);
Hill_1
plot(Hill_1,'-b');
[Hill_2,GOF_2] = fit(X(24:32),Y(24:32),f,'StartPoint',[1 -4],'Robust','Bisquare');
Hill_2
plot(Hill_2,'-r');
[Hill_3,GOF_3] = fit(X(14:22),Y(14:22),f,'StartPoint',[1 -4]);
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');

Kd_1_hill = Kd_1;
Kd_2_hill = Kd_2;

%% Logistic fit.
% A Dose-Response Logistic fit will give the Hill coefficient. This fit is
% based on the following equation:

% Response = (Max-Min) / (1+(Kd/Dose)nH)

% in which nH = Hill coefficient. If there is only one component, we have 4
% unknown:

% a = Kd, b = nH, c = max(y), d = min(y);

% if there are two components we have a total of 8 unknown:

% y = (c - d)/(1 + (a/x)b) (1 components) y = (c - d)/(1 + (a/x)b) + (g - h)/(1 + (e/x)f) (2 components)

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.2 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: this is a plot of the
% response (the bound ligand concentration), or the percentage of the
% maximal response (the "fractional saturation"), against the log of the
% dose (the free ligand concentration). This plot is 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 sigmoidal curve. After the logistic fit is carried out (previous
% cell), the data (and its fit, blue curve) can then be compared with the
% theoretical Logistic Dose-Response curves expected for nH = 2, nH = 1, nH
% = 0.5.

    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 a logistic fit. However,
% it is more difficult to get correct values for the low and high affinity
% Kd's.

%% 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 hyperbola plot, and we multiply by the expected
% kon to obtain the values of koff.
p1 = Kd_1_hyperb*10;
p2 = Kd_2_hyperb*10;
pin=[p1;p2;1];

%% 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 set the configuration of the solver
Stop = 10;
set(configset_m1, 'StopTime', Stop);
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.e-9);

% 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_2,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_2',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);
Kd_1_direct = u(1)/10
Kd_2_direct = u(2)/10

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

Bound = Bound*u(3);
Bound_R = Bound_R*u(3);
Bound_R_1 = Bound_R_1*u(3);


%%
% We can compare the residual from this fit with that we would obtain taking
% the best starting solution:

figure;plot(Unbound,Bound,'-b',Unbound,Product_vector,'or');hold on
% plot(Unbound,Bound_R','--k',Unbound,Bound_R_1','--m'); 

plot(Unbound,Hyperb2(Unbound),'--g');
plot(Unbound,Hill_1(Unbound),'--m');

hold off
legend('Direct Fit','Observed','Hyperbola fit','Logistic fit','Location','Best')
xlabel('[L] microM  ')
ylabel('Fractional saturation  ')

res_1 = (Product_vector - Bound);
sos_res_1 = res_1'*res_1
res_2 = (Product_vector - Hyperb2(Unbound));
sos_res_2 = res_2'*res_2
res_3 = (Product_vector - Hill_1(Unbound));
sos_res_3 = res_2'*res_3
cum_sos_res_1 = cumsum(res_1.*res_1);
cum_sos_res_2 = cumsum(res_2.*res_2);
cum_sos_res_3 = cumsum(res_3.*res_3);
figure;plot(Unbound,cum_sos_res_1,'-b',Unbound,cum_sos_res_2,'-g',Unbound,cum_sos_res_3,'-m')
%%
close all
% save(savefile);
% dlmwrite('../DATABASE/Problem_Binding_Curve.txt',[Conc_vector' Product_vector],'delimiter','\t', ...
%          'precision', 8)  
%% Conclusion
%
% 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. Logistic fit and dose-response curves both show good
% performance in detecting the presence of negative cooperativity. Direct
% kinetic fit provide the best Kd values.
%


