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

%% Two Substrates Bi-Bi Michaelis Menten Kinetics: Simulation and Analysis.

% Add to the path the enhanced scripts for drawing lines at specified
% coordinates.
% addpath(genpath('TOOLBOXES/ENZYME_KINETICS'));

global m1
global Time_points ntimes
global s1_cons s1_conc_vec s2_cons s2_conc_vec
global C C_init E E_init S1 S2 ES1 ES2 EP1 EP2 P1 P2
global kon1 koff1 kon2 koff2 kcatf kcatr kon3 koff3 kon4 koff4 

% Load the project
% sbioloadproject('Simple_Bi_Bi_Michaelis_Menten');
sbioloadproject('TOOLBOXES/ENZYME_KINETICS/Bi_Bi_Michaelis_Menten');

% Get information about the model.
sbioselect(m1,'Type','compartment')
sbioselect(m1,'Type','species')
sbioselect(m1,'Type','parameter')
sbioselect(m1,'Type','reaction')

% We also want to know the differential equations used in the model. This
% information may be very usefult for special applications (e.g., global
% fitting of multiple curves)
equations = getequations(m1);

% First, let's get the configuration of the simulation.
configset_m1 = getconfigset(m1);
get(configset_m1)

% Now we run a small simulation with the existing parameters.
bi_bi_kinetics = sbiosimulate(m1);

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

% We can plot the simulation
% sbioplot(bi_bi_kinetics)

% If we want to see the actual numbers, the time points are in the array
% bi_bi_kinetics.Time and the the time-course of the simulation is in the 
% bi_bi_kinetics.Data array
% openvar bi_bi_kinetics.Time
% openvar bi_bi_kinetics.Data

%%
close all

% Then we extend the simulation time, we set the solver and request a
% smaller tolerance
stop_time = 4000;

set(configset_m1, 'StopTime', stop_time);
set(configset_m1,'SolverType','ode15s');
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.0e-9);

% And we configure the solver so that states are logged only at
% specific times.
time_vec = logspace(log10(1),log10(stop_time),35);

Time_points = [0 time_vec]';

ntimes = size(Time_points,1);

set(configset_m1.SolverOptions, 'OutputTimes',Time_points);

get(configset_m1)

% We rerun the simulation with the new StopTime.
bi_bi_kinetics = sbiosimulate(m1);

% We can plot the  new simulation
% 
sbioplot(bi_bi_kinetics)

% bi_bi_kinetics.Data has the following fields:
% S1 E ES1 EP1 S2 ES2 EP2 P2 P1
bi_bi_kinetics.DataNames 

% It looks OK. 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');
E = sbioselect(m1,'Name','E');
S1 = sbioselect(m1,'Name','S1');
S2 = sbioselect(m1,'Name','S2');
ES1 = sbioselect(m1,'Name','ES1');
ES2 = sbioselect(m1,'Name','ES2');
EP1 = sbioselect(m1,'Name','EP1');
EP2 = sbioselect(m1,'Name','EP2');
P1 = sbioselect(m1,'Name','P1');
P2 = sbioselect(m1,'Name','P2');
kon1 = sbioselect(m1,'Name','kon1');
koff1 = sbioselect(m1,'Name','koff1');
kon2 = sbioselect(m1,'Name','kon2');
koff2 = sbioselect(m1,'Name','koff2');
kcatf = sbioselect(m1,'Name','kcatf');
kcatr = sbioselect(m1,'Name','kcatr');
kon3 = sbioselect(m1,'Name','kon3');
koff3 = sbioselect(m1,'Name','koff3');
kon4 = sbioselect(m1,'Name','kon4');
koff4 = sbioselect(m1,'Name','koff4');

% We store the initial value for the enzyme and the cell volume.
C_init = C.capacity;
E_init = E.InitialAmount;

%%
close all

% Now we can set different values. We notice that Kd1 = koff1/kon1 = 200/10
% = 20 microM, so we can try a concentration range from 1/10 the Kd to 20
% times the Kd, that is from 2 to 400 microM. Since Kd2 = koff2/kon2 =
% 50/10 = 5 microM, so we can try a concentration range from 0.5 to 100
% microM.

% As usual, in order to examine the effect of different concentrations of
% the two substrates we write a 'loop'. First we define a vector with all
% the concentrations we want to explore spaced logarithmically.

s1_cons = 10;
s2_cons = 10;
s1_conc_vec = logspace(log10(2),log10(400),s1_cons);
s2_conc_vec = logspace(log10(0.5),log10(100),s2_cons);
s2_cons = length(s2_conc_vec);

Initial_Rate = zeros(s1_cons,s2_cons);
Product_mat = nan(ntimes,s1_cons,s2_cons);
Time_mat = nan(ntimes,s1_cons);
String_mat_s1 = cell(s1_cons,1);
String_mat_s2 = cell(s2_cons,1);
Init_mat = nan(s1_cons,4,s2_cons);

for k = 1:s2_cons
    % Here we just get the legend entries.
    String_mat_s2{k,1} = [num2str(s2_conc_vec(k),'%6.2f') ' µM']; 
    for i = 1:s1_cons
        String_mat_s1{i,1} = [num2str(s1_conc_vec(i),'%6.2f') ' µM'];         
        error = (rand(4,1)-0.5);
        error(1) = 1 + error(1)*0.05;
        error(2) = 1 + error(2)*0.025;
        error(3) = 1 + error(3)*0.025;
        error(4) = 1 + error(4)*0.025;        
        % error = ones(4,1); % no error in any set up (perfect data).
        % error(4) = 1; % no error in the addition of the inhibitor.
        set(C,'Capacity',C_init*error(1));
        set(E,'InitialAmount',E_init*error(2)/C.Capacity);
        set(S1,'InitialAmount',s1_conc_vec(i)*error(3)/C.Capacity);
        set(S2,'InitialAmount',s2_conc_vec(k)*error(4)/C.Capacity);
        Init_mat(i,1,k) = C.Capacity;
        Init_mat(i,2,k) = E.InitialAmount;
        Init_mat(i,3,k) = S1.InitialAmount;
        Init_mat(i,4,k) = S2.InitialAmount;
        
        clear error
        
        bi_bi_kinetics = sbiosimulate(m1);
        % Find the P2 product (8th field in the data array) concentration 
        % at 1 second.
        Product = bi_bi_kinetics.Data(2,8);
    	Time = bi_bi_kinetics.Time(2);
        % Fill the 2D array with times and product values for each run.
        Product_mat(:,i,k) = bi_bi_kinetics.Data(:,8);
        Time_mat(:,i) = bi_bi_kinetics.Time(:,1);
        % Initial_Rate(i,k) = Product/Time;
    end
end

%% Calculation of all initial rates

% Initial_Rate = zeros(s1_cons,s2_cons);

for k = 1:s2_cons
    Initial_Rate(:,k) = (Product_mat(2,:,k)./Time_mat(2,:))';
end

%%
TWO_SUBSTRATES_MM_1 = figure;
    set(TWO_SUBSTRATES_MM_1,'Units','normalized','Position',[0.0 0.4 1.0 0.6],...
        'Name','Two Substrtates MM Kinetics: linear regression');clf
    
% Now we can choose various ways of analyzing the data. We start by
% plotting the initial velocity of the reaction for each initial
% concentration of the substrate at the different concentration of 
% inhibitor.
    

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

%     plot(s1_conc_vec,Initial_Rate(:,1:2:s2_cons));
%     legend({String_mat{1:2:s2_cons}},'Location','Best')
    plot(s1_conc_vec,Initial_Rate);
%    legend({String_mat{1:2:s2_cons}},'Location','Best')
    legend(String_mat_s2,'Location','Best')
    xlim([-10,1.1*max(s1_conc_vec)]);
    ylim([0,1.1*max(Initial_Rate(:))]);
    xlabel('[S1] (µM)');
    ylabel('Vo (µM/sec)');

% The following is the traditional Lineweaver-Burke plot.
%
    subplot2 = subplot(1,3,2,'Parent',figure(gcf));
    box(subplot2,'on');
    % grid(subplot4,'on');
    hold(subplot2,'all');

Vmax1 = zeros(s2_cons,1);
Km1 = zeros(s2_cons,1);
slope1 = zeros(s2_cons,1); 
R2 = zeros(s2_cons,1);
LB_x = 1./s1_conc_vec;
LB_y = 1./Initial_Rate;
a = zeros(s2_cons,1);
b = zeros(s2_cons,1);
        
colors = get(gca,'ColorOrder');
colors = repmat(colors,100,1);

low_x = -max(LB_x)/5;
high_x = 1.1*max(LB_x);
low_y = -10*min(LB_y(:));
high_y = 1.1*max(LB_y(:));

for i = 1:s2_cons
    
% Here we use Matlab 'fit' function with the 'robust' option as encoded in 
% the local function 'Robust_Linear_Least_Squares'.    

[FITRESULT,GOF] = Robust_Linear_Least_Squares(LB_x,LB_y(:,i)','Bisquare');

% Check the fit quality:
fit_params = coeffvalues(FITRESULT);
a(i) = fit_params(1);
b(i) = fit_params(2);
R2(i) = GOF.rsquare;

Vmax1(i) = 1/b(i);   % 1/fit_params(2);
Km1(i) = a(i)/b(i);  % fit_params(1)/fit_params(2);
slope1(i) = a(i);
end

% for i = 1:2:s2_cons
for i = 1:s2_cons
x=[low_x,high_x];
y=[low_x*a(i)+b(i),high_x*a(i)+b(i)];

plot(LB_x,LB_y(:,i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
hold on
line(x,y,'Color',colors(i,:));
end

xlim([low_x,high_x]);
ylim([low_y,high_y]);
% ylim([-0.15E3,high_y]);

vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});

% String_mat2 = cell(s2_cons+1,1);
% vec2 = [1 3:2:s2_cons];
% String_mat2(1:2:s2_cons) = String_mat(vec2);
% 
% for i = 2:2:s2_cons+1
% String_mat2{i,1} = 'Fit';
% end

String_mat1 = cell(s2_cons*2,1);
String_mat1(1:2:(s2_cons*2-1)) = String_mat_s2;

for i = 2:2:s2_cons*2
String_mat1{i,1} = 'Fit';
end

legend(String_mat1,'Location','Best')

xlabel('1/[S1]');
ylabel('1/V');


% We can get a similar reciprocal plot changing the substrate in
% the inner and outer loops.

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

Vmax2 = zeros(s1_cons,1);
Km2 = zeros(s1_cons,1);
slope2 = zeros(s1_cons,1); 
R2 = zeros(s1_cons,1);
LB_x = 1./s2_conc_vec;
% We take the transpose of the 'Initial_Rate' matrix
LB_y = 1./Initial_Rate';
a = zeros(s1_cons,1);
b = zeros(s1_cons,1);
        
colors = get(gca,'ColorOrder');
colors = repmat(colors,100,1);

low_x = -max(LB_x)/5;
high_x = 1.1*max(LB_x);
low_y = -10*min(LB_y(:));
high_y = 1.1*max(LB_y(:));

for i = 1:s1_cons
    
% Here we use Matlab 'fit' function with the 'robust' option as encoded in 
% the local function 'Robust_Linear_Least_Squares'.    

[FITRESULT,GOF] = Robust_Linear_Least_Squares(LB_x,LB_y(:,i)','Bisquare');

% Check the fit quality:
fit_params = coeffvalues(FITRESULT);
a(i) = fit_params(1);
b(i) = fit_params(2);
R2(i) = GOF.rsquare;

Vmax2(i) = 1/b(i);   % 1/fit_params(2);
Km2(i) = a(i)/b(i);  % fit_params(1)/fit_params(2);
slope2(i) = a(i);
end

% for i = 1:2:s2_cons
for i = 1:s1_cons
x=[low_x,high_x];
y=[low_x*a(i)+b(i),high_x*a(i)+b(i)];

plot(LB_x,LB_y(:,i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
hold on
line(x,y,'Color',colors(i,:));
end

xlim([low_x,high_x]);
ylim([low_y,high_y]);
% ylim([-0.15E3,high_y]);

vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});

% String_mat2 = cell(s2_cons+1,1);
% vec2 = [1 3:2:s2_cons];
% String_mat2(1:2:s2_cons) = String_mat(vec2);
% 
% for i = 2:2:s2_cons+1
% String_mat2{i,1} = 'Fit';
% end

String_mat2 = cell(s1_cons*2,1);
String_mat2(1:2:(s1_cons*2-1)) = String_mat_s1;

for i = 2:2:s1_cons*2
String_mat2{i,1} = 'Fit';
end

legend(String_mat2,'Location','Best')

xlabel('1/[S2]');
ylabel('1/V');


%%
TWO_SUBSTRATES_MM_2 = figure;
    set(TWO_SUBSTRATES_MM_2,'Units','normalized','Position',[0.5 0.5 0.5 0.5],...
        'Name','Two Substrtates MM Kinetics: linear regression');clf

% Here we derive Vmax, Km_S1 and Km_S2 values from 4 different replots of 
% the first two reciprocal plots.

    subplot1 = subplot(2,2,1,'Parent',figure(gcf));
    box(subplot1,'on');
    grid(subplot1,'on');
    hold(subplot1,'all');    
    
[FITRESULT,~] = Robust_Linear_Least_Squares(1./s2_conc_vec,1./Vmax1','bisquare');

% Check the fit quality:
fit_params = coeffvalues(FITRESULT);
a_i = fit_params(1);
b_i = fit_params(2);
% R2_i = GOF.rsquare;

Vmax_S2 = 1/b_i;
Km_S2 = a_i/b_i;
Km_S2_1 = Km_S2;

% low_x = min(1./s2_conc_vec);
low_x = 2*(-1/Km_S2);
high_x = 1.1*max(1./s2_conc_vec);
low_y = low_x*a_i+b_i;
high_y = high_x*a_i+b_i;

x=[low_x,high_x];
y=[low_y,high_y];


    for i = 1:s2_cons
    plot(1./s2_conc_vec(i),1./Vmax1(i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
         hold on
    end

line(x,y,'Color','cyan','LineStyle','-');

vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});

xlim([-2/Km_S2,1.1*high_x]);
% ylim([low_y,1.1*high_y]);
ylim([-1/Vmax_S2,1.1*high_y]);

% legend({String_mat{1:2:s2_cons}},'Location','Best')
% Axes label
xlabel('1/[S2] (µM)');
ylabel('1/Vmax (µM/s)');

% text
string1 = ['Vmax = ' num2str(Vmax_S2,'%6.2e'),' (µM/sec)  '];
string2 = ['KmS2 = ' num2str(Km_S2,'%6.2e'),' (µM)      '];
Result = [string1;string2];

annotation(gcf,'textbox',...
    [0.3 0.5 0.2 0.2],...
    'String',Result,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

clear low_x high_x low_y high_y

%
    subplot3 = subplot(2,2,3,'Parent',figure(gcf));
    box(subplot3,'on');
    grid(subplot3,'on');
    hold(subplot3,'all');    
    
[FITRESULT,~] = Robust_Linear_Least_Squares(1./s1_conc_vec,1./Vmax2','bisquare');

% Check the fit quality:
fit_params = coeffvalues(FITRESULT);
a_i = fit_params(1);
b_i = fit_params(2);
% R2_i = GOF.rsquare;

Vmax_S1 = 1/b_i;
Km_S1 = a_i/b_i;
Km_S1_1 = Km_S1;

% low_x = min(1./s1_conc_vec);
low_x = -1/Km_S1;
high_x = 1.1*max(1./s1_conc_vec);
low_y = low_x*a_i+b_i;
high_y = high_x*a_i+b_i;

x=[low_x,high_x];
y=[low_y,high_y];


    for i = 1:s1_cons
    plot(1./s1_conc_vec(i),1./Vmax2(i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
         hold on
    end

line(x,y,'Color','cyan','LineStyle','-');

vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});

xlim([-1.2/Km_S1,1.1*high_x]);
% ylim([low_y,1.1*high_y]);
ylim([-0.2/Vmax_S1,1.3*high_y]);

% legend({String_mat{1:2:s2_cons}},'Location','Best')
% Axes label
xlabel('1/[S1] (µM)');
ylabel('1/Vmax (µM/s)');

% text
string1 = ['Vmax = ' num2str(Vmax_S1,'%6.2e'),' (µM/sec)  '];
string2 = ['KmS1 = ' num2str(Km_S1,'%6.2e'),' (µM)      '];
Result = [string1;string2];

annotation(gcf,'textbox',...
    [0.15 0.2 0.2 0.2],...
    'String',Result,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

clear low_x high_x low_y high_y

%
    subplot2 = subplot(2,2,2,'Parent',figure(gcf));
    box(subplot2,'on');
    grid(subplot2,'on');
    hold(subplot2,'all');    
    
[FITRESULT,~] = Robust_Linear_Least_Squares(1./s2_conc_vec,slope1','bisquare');

% Check the fit quality:
fit_params = coeffvalues(FITRESULT);
a_i = fit_params(1);
b_i = fit_params(2);
% R2_i = GOF.rsquare;

Km_S1oVmax = b_i;
KiS1xKmS2oKmS1 = a_i/b_i;
Ki_S1_1 = KiS1xKmS2oKmS1*Km_S1/Km_S2;
Km_S1_2 = Km_S1oVmax*Vmax_S2;

% low_x = min(1./s2_conc_vec);
low_x = -1/KiS1xKmS2oKmS1;
high_x = 1.1*max(1./s2_conc_vec);
low_y = low_x*a_i+b_i;
high_y = high_x*a_i+b_i;

x=[low_x,high_x];
y=[low_y,high_y];


    for i = 1:s2_cons
    plot(1./s2_conc_vec(i),slope1(i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
         hold on
    end

line(x,y,'Color','cyan','LineStyle','-');

vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});

xlim([-10/KiS1xKmS2oKmS1,1.1*high_x]);
% ylim([low_y,1.1*high_y]);
ylim([-2*b_i,1.1*high_y]);

% legend({String_mat{1:2:s2_cons}},'Location','Best')
% Axes label
xlabel('1/S2 (µM)');
ylabel('slope\_S1');

% text
string1 = ['KmS1/Vmax =' num2str(Km_S1oVmax,'%6.2e'),' (µM/sec)'];
string2 = ['KiS1 =   ' num2str(KiS1xKmS2oKmS1*Km_S1/Km_S2,'%6.2e'),' (µM)      '];
string3 = ['KmS1 =   ' num2str(Km_S1oVmax*Vmax_S2,'%6.2e'),' (µM)      '];
Result = [string1;string2;string3];

annotation(gcf,'textbox',...
    [0.7 0.5 0.2 0.2],...
    'String',Result,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

clear low_x high_x low_y high_y

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

[FITRESULT,~] = Robust_Linear_Least_Squares(1./s1_conc_vec,slope2','bisquare');

% Check the fit quality:
fit_params = coeffvalues(FITRESULT);
a_i = fit_params(1);
b_i = fit_params(2);
R2_i = GOF.rsquare;

Km_S2oVmax = b_i;
Ki_S1 = a_i/b_i;
Ki_S1_2 = Ki_S1;
Km_S2_2 = Km_S2oVmax*Vmax_S1;

%low_x = min(1./s1_conc_vec);
low_x = -1/Ki_S1;
high_x = 1.1*max(1./s1_conc_vec);
low_y = low_x*a_i+b_i;
high_y = high_x*a_i+b_i;

x=[low_x,high_x];
y=[low_y,high_y];


    for i = 1:s1_cons
    plot(1./s1_conc_vec(i),slope2(i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
         hold on
    end

    line(x,y,'Color','cyan','LineStyle','-');

vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});

xlim([-2/Ki_S1,1.1*high_x]);
% ylim([low_y,1.1*high_y]);
ylim([-b_i,1.1*high_y]);

% legend({String_mat{1:2:s2_cons}},'Location','Best')
% Axes label
xlabel('1/S1 (µM)');
ylabel('slope\_S2');

% text
string1 = ['KmS2/Vmax =' num2str(Km_S2oVmax,'%6.2e'),' (µM/sec)'];
string2 = ['KiS1 =   ' num2str(Ki_S1,'%6.2e'),' (µM)      '];
string3 = ['KmS2 =   ' num2str(Km_S2oVmax*Vmax_S1,'%6.2e'),' (µM)      '];
Result = [string1;string2;string3];

annotation(gcf,'textbox',...
    [0.6 0.2 0.2 0.2],...
    'String',Result,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

%%
% Finally we derive approximate values for some of the rate constants.
kcatf.Value = mean([Vmax_S1 Vmax_S2])
koff1.Value = mean([Ki_S1_1 Ki_S1_2])*kon1.Value - kcatf.Value
koff2.Value = mean([Km_S2_1 Km_S2_2])*kon2.Value - kcatf.Value
% For the other variable we will guess some reasonable values:
kcatr.Value = 0.5*kcatf.Value
koff3.Value = 2*koff2.Value
koff4.Value = 2*koff1.Value
% 
%% Global fit
% Here we fit globally all the progress curves with a single set of 
% parameters. However, this time we can take take only the relevant part 
% of the data. 

RelTime_mat = nan(ntimes,s1_cons,s2_cons);
RelProduct_mat = nan(ntimes,s1_cons,s2_cons);

for k = 1:s2_cons
    for i = 1:s1_cons
        X = Time_mat(:,i);
        Y = Product_mat(:,i,k);
        %
        first_der = gradient(Y,X);
        second_der = gradient(first_der,X);
        smooth_2der = smooth(-second_der,10);
        range_ind = find(abs(smooth_2der) >= 0.05*max(smooth_2der));
        % plot(Time_mat(range_ind,i),Product_mat(range_ind,i,k),'-');
        % hold on
        RelTime_mat(range_ind,i,k) = Time_mat(range_ind,i);
        RelProduct_mat(range_ind,i,k) = Product_mat(range_ind,i,k);
        clear X Y
    end
end

%% Final check of the data before the fit
% If we don't like the trimming of the curves we can redo it, for example 
% increasing or decreasing the fractionn of the curves we keep, or we can 
% reset everything to the original values:

%     for k = 1:s2_cons
%         RelTime_mat(:,:,k) = Time_mat(:,:);
%         RelProduct_mat(:,:,k) = Product_mat(:,:,k);
%     end

%% Here we plot the progress curves within their modified ranges.
    
for k = 1:s2_cons
plot(RelTime_mat(:,:,k),RelProduct_mat(:,:,k),'.-');
hold on
end
% legend(gca,String_mat,'Location','BestOutside')
xlabel('Time (sec)');
ylabel('[P2] (µM)');


%% 
% Here we convert all the progress curves into a single consecutive
% progress curve.

global Time_vector_ind Product_vector_ind

Product_vector = RelProduct_mat(:);
Time_vector = RelTime_mat(:);
Time_vector_ind = ~isnan(Time_vector);
Product_vector_ind = ~isnan(Product_vector);
Time_vector = Time_vector(Time_vector_ind);
Product_vector = Product_vector(Product_vector_ind);
plot(Product_vector,'-r')
% Check:
Comb_vector = [Time_vector Product_vector];
openvar Comb_vector

%% Choice of parameters to refine in the 1st refinement round
% We need to set some initial value for the refinement of parameters: we
% will refine koff1, kcat, and koff2. For example, we can take the initial 
% values from the previous fit
p1 = koff3.Value;
p2 = koff4.Value;
p3 = kcatr.Value;
pin = [p1;p2;p3];
% We reinitialize the initial value for the enzyme
% concentration and the cell volume.
E.InitialAmount = E_init;
C.capacity = C_init;

%% Non-linear global fit of all the progress curves.
nobs = length(Time_vector);
nvars = length(pin);

% options = ...
% optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-8,'TolX',1e-8);
% [FP1,sos1,R1,~,~,~,J1] = ...
% lsqcurvefit(@bi_bi_fit_1,pin,Time_vector,Product_vector,...
% [],[],options); 
% 
% MSE1 = sos1/(nobs-nvars);
% RMSE1 = sqrt(MSE1)
% 
% % We calculate the covariance matrix of the variables from the last value
% % of the Jacobian: J'*J is the 'precision' matrix, its inverse is the
% % 'covariance' matrix (LECTURE 9).
% 
% J1 = full(J1);
% Cov1 = inv(J1'*J1).*MSE1;
% 
% % The covariance matrix can be converted into a correlation matrix; also
% % the diagonal of the covariance matrix gives the variances and therefore
% % the sigmas.
% 
% % [Corr,sigma] = corrcov(Cov2);
% % FP2ci = [FP2-1.96*sigma FP2+1.96*sigma] 
% 
% % The following works for both nlinfit and lsqcurvefit.
% [fcurve1,delta1] = nlpredci('bi_bi_fit_1',Time_vector,FP1,R1,'Covar',Cov1);
% FP1ci = nlparci(FP1,R1,'covar',Cov1);

% 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-8,'Display','iter','RobustWgtFun','',...
    'UseParallel',true,'MaxIter',200);
[FP1,R1,J1,Cov1,MSE1] = ...
    nlinfit(Time_vector,Product_vector,'bi_bi_fit_1',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
[fcurve1,delta1] = nlpredci('bi_bi_fit_1',Time_vector,FP1,R1,'Covar',Cov1);
FP1ci = nlparci(FP1,R1,'covar',Cov1);
% The following alternative syntax gives the same result as the previous 
% two lines.
% [fcurve1,delta1] = nlpredci('fit5fitM',lt,FP1,R1,'jacobian',J1);
% FPci = nlparci(FP1,R1,'jacobian',J1);
Comb2_vector = [Comb_vector fcurve1];
openvar Comb2_vector
plot(Comb2_vector(:,1),Comb2_vector(:,2),'.r',Comb2_vector(:,1),Comb2_vector(:,3),'-b')
% plot(Comb2_vector(:,2),'.r');hold on
% plot(Comb2_vector(:,3),'-b')

%% Choice of parameters to refine in the 2nd refinement round
% We need to set some initial value for the refinement of parameters: we
% will refine koff1, kcat, and koff2. For example, we can take the initial 
% values from the previous fit
p1 = koff1.Value;
p2 = koff2.Value;
p3 = koff3.Value;
p4 = koff4.Value;
p5 = kcatf.Value;
p6 = kcatr.Value;
pin = [p1;p2;p3;p4;p5;p6];

% We reinitialize the initial value for the enzyme
% concentration and the cell volume.
E.InitialAmount = E_init;
C.capacity = C_init;

%% Non-linear global fit of all the progress curves.
% nobs = length(Time_vector);
% nvars = length(pin);
% 
% options = ...
% optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-8,'TolX',1e-8);
% [FP1,sos1,R1,~,~,~,J1] = ...
% lsqcurvefit(@bi_bi_fit,pin,Time_vector,Product_vector,...
% [],[],options); 
% 
% MSE1 = sos1/(nobs-nvars);
% RMSE1 = sqrt(MSE1)
% 
% % We calculate the covariance matrix of the variables from the last value
% % of the Jacobian: J'*J is the 'precision' matrix, its inverse is the
% % 'covariance' matrix (LECTURE 9).
% 
% J1 = full(J1);
% Cov1 = inv(J1'*J1).*MSE1;
% 
% % The covariance matrix can be converted into a correlation matrix; also
% % the diagonal of the covariance matrix gives the variances and therefore
% % the sigmas.
% 
% % [Corr,sigma] = corrcov(Cov2);
% % FP2ci = [FP2-1.96*sigma FP2+1.96*sigma] 
% 
% % The following works for both nlinfit and lsqcurvefit.
% [fcurve1,delta1] = nlpredci('bi_bi_fit',Time_vector,FP1,R1,'Covar',Cov1);
% FP1ci = nlparci(FP1,R1,'covar',Cov1);

% 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-8,'Display','iter','RobustWgtFun','',...
    'UseParallel',true,'MaxIter',200);
[FP1,R1,J1,Cov1,MSE1] = ...
    nlinfit(Time_vector,Product_vector,'bi_bi_fit_2',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
[fcurve1,delta1] = nlpredci('bi_bi_fit_2',Time_vector,FP1,R1,'Covar',Cov1);
FP1ci = nlparci(FP1,R1,'covar',Cov1);
% The following alternative syntax gives the same result as the previous 
% two lines.
% [fcurve1,delta1] = nlpredci('fit5fitM',lt,FP1,R1,'jacobian',J1);
% FPci = nlparci(FP1,R1,'jacobian',J1);
Comb2_vector = [Comb_vector fcurve1];
openvar Comb2_vector
plot(Comb2_vector(:,1),Comb2_vector(:,2),'.r',Comb2_vector(:,1),Comb2_vector(:,3),'-b')
% plot(Comb2_vector(:,2),'.r');hold on
% plot(Comb2_vector(:,3),'-b')

%% Regeneration of the calculated individual progress curves.
koff1.Value = FP1(1);
koff2.Value = FP1(2);
koff3.Value = FP1(3);
koff4.Value = FP1(4);
kcatf.Value = FP1(5);
kcatr.Value = FP1(6);

calc_product_mat = nan(ntimes,s1_cons,s2_cons);

for k = 1:s2_cons
    for i = 1:s1_cons
        set(C,'Capacity',C_init);
        set(E,'InitialAmount',E_init);
        set(S1,'InitialAmount',s1_conc_vec(i));
        set(S2,'InitialAmount',s2_conc_vec(k));
        
        bi_bi_kinetics = sbiosimulate(m1);
        % Find the P2 product (8th field in the data array) concentration 
        % at 1 second.
        Product = bi_bi_kinetics.Data(2,8);
    	% Time = bi_bi_kinetics.Time(2);
        % Fill the 2D array with times and product values for each run.
        calc_product_mat(:,i,k) = bi_bi_kinetics.Data(:,8);
        % calc_initial_rate(i,k) = Product/Time;
    end
end

%%
% Here we plot some of the layer of product_mat and calc_product_mat for comparison:
for k = 1:s2_cons
plot(Time_mat,Product_mat(:,:,k),'.',Time_mat,calc_product_mat(:,:,k),'-');
hold on
end
hold off

% or just a single one:
% plot(Time_mat,Product_mat(:,:,6),'.',Time_mat,calc_product_mat(:,:,6),'-');

%% Backward check
% Recalc_Product_vector = zeros(ntimes*s1_cons*s2_cons,1);
% % Recalc_Time_vector = zeros(ntimes*s1_cons*s2_cons,1);
% for i = 1:s2_cons
%     for k = 1:s1_cons
%         for j = 1:ntimes
%         Recalc_Product_vector(j+(k-1)*ntimes+(i-1)*s1_cons*ntimes,1) = calc_product_mat(j,k,i);
%         end;
%     end;
% end
Recalc_Product_vector = calc_product_mat(:);
Recalc_Product_vector = Recalc_Product_vector(Product_vector_ind);
% plot(Time_vector,Product_vector,'.r')
Comb3_vector = [Comb2_vector Recalc_Product_vector];
openvar Comb3_vector

%% Regeneration of the calculated individual progress curves (alternative way).
% calc_product_mat = zeros(ntimes,s1_cons,s2_cons);
% for i = 1:s2_cons
%     for k=1:s1_cons
%         for j=1:ntimes
%         calc_product_mat(j,k,i)=fcurve1(j+(k-1)*ntimes+(i-1)*s1_cons*ntimes,1);
%         end
%     end
% end
% 
% % Here we derive the calculated initial velocities
% calc_initial_rate = zeros(s1_cons,s2_cons);
% for k = 1:s2_cons
% calc_initial_rate(:,k) = (calc_product_mat(2,:,k)./Time_mat(2,:))';
% end

% Check a few curves.
% plot(Time_mat,Product_mat(:,:,5),'+');
% hold on
% colors = get(gca,'ColorOrder');
% plot(Time_mat,calc_product_mat(:,:,5));

%% Recalculation of all initial rates

obs_initial_rate = zeros(s1_cons,s2_cons);
calc_initial_rate = zeros(s1_cons,s2_cons);
for k = 1:s2_cons
    obs_initial_rate(:,k) = (Product_mat(2,:,k)./Time_mat(2,:))';
    calc_initial_rate(:,k) = (calc_product_mat(2,:,k)./Time_mat(2,:))';
end
    
%%
TWO_SUBSTRATES_MM_3 = figure;
    set(TWO_SUBSTRATES_MM_3,'Units','normalized','Position',[0.4 0.4 0.6 0.6],...
        'Name','Two Substrates MM Kinetics: non-linear global fit');clf
    
%     subplot1 = subplot(1,1,1,'Parent',figure(gcf));
%     box(subplot1,'on');
%     grid(subplot1,'on');
%     hold(subplot1,'all');    

colors = get(gca,'ColorOrder');
colors = repmat(colors,100,1);

for i = 1:s2_cons
plot(s1_conc_vec,obs_initial_rate(:,i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
hold on
plot(s1_conc_vec,calc_initial_rate(:,i),'--',...
             'Color',colors(i,:));
         
end

    xlim([-10,1.1*max(s1_conc_vec)]);
    ylim([-0.001,1.1*max(obs_initial_rate(:))]);
    xlabel('[S1] (µM)');
    ylabel('Vo (µM/sec)');

String_mat2 = cell(s2_cons*2,1);
String_mat2(1:2:(s2_cons*2-1)) = String_mat_s2;

for i = 2:2:s2_cons*2
String_mat2{i,1} = 'Fit';
end

legend(String_mat2,'Location','BestOutside')

hline(0,{'g-', 'LineWidth', 0.5});

% text
koff1_fit = FP1(1);
koff2_fit = FP1(2);
koff3_fit = FP1(3);
koff4_fit = FP1(4);
% kcatf_fit = FP1(5);
% Vmax = kcatf.Value - kcatr_fit;
Vmax = kcatf.Value;
Kd_S1 = koff1_fit/kon1.Value;
Kd_S2 = koff2_fit/kon2.Value;
Kd_P1 = koff3_fit/kon3.Value;
Kd_P2 = koff4_fit/kon4.Value;
string1 = ['Vmax =   ' num2str(Vmax,'%6.2e'),' (µM/sec)'];
string2 = ['KdS1 =   ' num2str(Kd_S1,'%6.2e'),' (µM)    '];
string3 = ['KdS2 =   ' num2str(Kd_S2,'%6.2e'),' (µM)    '];
string4 = ['KdP1 =   ' num2str(Kd_P1,'%6.2e'),' (µM)    '];
string5 = ['KdP2 =   ' num2str(Kd_P2,'%6.2e'),' (µM)    '];
Result = [string1;string2;string3;string4;string5];

annotation(gcf,'textbox',...
    [0.14 0.71 0.2 0.2],...
    'String',Result,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

%% Plot correlation between parameters
PARAMETERS_CORRELATION_1 = figure;
    set(PARAMETERS_CORRELATION_1,'Units','normalized','Position',[0.4 0.6 0.6 0.4],...
        'Name','Correlation and Partial Correlation between Parameters');clf
    
% We can also calculate the partial correlation matrix
    N = length(pin);
    rho = zeros(N,N);
    cov_inv = inv(Cov1);

    for i = 1:N
        for j = i:N
            rho(i,j) = -cov_inv(i,j)/sqrt(cov_inv(i,i)*cov_inv(j,j));
            rho(j,i) = rho(i,j);
        end
    end
    
subplot1 = subplot(1,2,1,'Parent',figure(gcf));
% box(subplot1,'on');
% grid(subplot1,'on');
% hold(subplot1,'all');

imagesc(Corr1);
set(gca,'YDir','normal');
xlabel('Variable Number');
ylabel('Variable Number');
colorbar    
title('Correlation');
    
subplot2 = subplot(1,2,2,'Parent',figure(gcf));

imagesc(rho)
set(gca,'YDir','normal');
xlabel('Variable Number');
ylabel('Variable Number');
colorbar
title('Partial Correlation');

% vline(3.5,{'m-', 'LineWidth', 0.5});
% vline(23.5,{'m-', 'LineWidth', 0.5});
% vline(43.5,{'m-', 'LineWidth', 0.5});
% hline(3.5,{'m-', 'LineWidth', 0.5});
% hline(23.5,{'m-', 'LineWidth', 0.5});
% hline(43.5,{'m-', 'LineWidth', 0.5});


