% 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-Uni Michaelis Menten Kinetics: Simulation and Analysis.

% Add to the path the enhanced scripts for drawing lines at specified
% coordinates.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
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 P kon1 koff1 kcat kon2 koff2

% Load the project
sbioloadproject('../TOOLBOXES/ENZYME_KINETICS/Bi_Uni_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_uni_kinetics = sbiosimulate(m1);
% 
% % We can get some information about the simulation.
% get(bi_uni_kinetics)
% 
% % We can plot the simulation
% sbioplot(bi_uni_kinetics)
% 
% % If we want to see the actual numbers, the time points are in the array
% % bi_uni_kinetics.Time and the the time-course of the simulation is in the 
% % bi_uni_kinetics.Data array
% % openvar bi_uni_kinetics.Time
% % openvar bi_uni_kinetics.Data

%%
close all

% Then we set the simulation time and the solver 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 1 10 100:100:stop_time]';
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_uni_kinetics = sbiosimulate(m1);

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

% 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');
S1 = sbioselect(m1,'Name','S1');
E = sbioselect(m1,'Name','E');
S2 = sbioselect(m1,'Name','S2');
ES1 = sbioselect(m1,'Name','ES1');
ES2 = sbioselect(m1,'Name','ES2');
P = sbioselect(m1,'Name','P');
kon1 = sbioselect(m1,'Name','kon1');
koff1 = sbioselect(m1,'Name','koff1');
kon2 = sbioselect(m1,'Name','kon2');
koff2 = sbioselect(m1,'Name','koff2');
kcat = sbioselect(m1,'Name','kcat');

% 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_conc_vec =[0 s2_conc_vec];
s2_cons = length(s2_conc_vec);
% openvar s1_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;
        
        bi_uni_kinetics = sbiosimulate(m1);
        % Find the product concentration at 1 second.
        Product = bi_uni_kinetics.Data(2,4);
    	Time = bi_uni_kinetics.Time(2);
        % Fill the 2D array with times and product values for each run.
        Product_mat(:,i,k) = bi_uni_kinetics.Data(:,4);
        Time_mat(:,i) = bi_uni_kinetics.Time(:,1);
        Initial_Rate(i,k) = Product/Time;
    end
end

for k = 1:s2_cons
    for i = 1:s1_cons
        Substrate_mat(i,1,k) = s1_conc_vec(i);
        Substrate_mat(i,2,k) = s2_conc_vec(k);        
    end
end

dlmwrite('Bi_Uni_Michaelis_Menten_Time_mat.txt',Time_mat,'delimiter','\t');
dlmwrite('Bi_Uni_Michaelis_Menten_Product_mat.txt',Product_mat,'delimiter','\t');
dlmwrite('Bi_Uni_Michaelis_Menten_Substrate_mat.txt',Substrate_mat,'delimiter','\t');

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

% 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

%
    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;
Ki_S2 = a_i/b_i;

% low_x = min(1./s2_conc_vec);
low_x = -1/Ki_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),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/Ki_S2,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 = ['KiS2 =   ' num2str(Ki_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

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

% 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

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

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

%%
% text
% Finally we derive approximate values for some of the rate constants.
Km_S1 = mean([Km_S1 Km_S1oVmax*Vmax_S2])
Km_S2 = mean([Km_S2 Km_S2oVmax*Vmax_S1])
kcat.Value = mean([Vmax_S1 Vmax_S2])
koff1.Value = Ki_S1*kon1.Value
koff2.Value = Km_S2*kon2.Value - kcat.Value

kcat_bk = mean([Vmax_S1 Vmax_S2]);
koff1_bk = Ki_S1*kon1.Value;
koff2_bk = Km_S2*kon2.Value - kcat.Value;

%% Global fit
close all

% Alternatively we could fit globally all the progress curves with a single
% set of parameters. For this purpose we will first convert all the
% progress curves into a single consecutive progress curve. 

Product_vector = Product_mat(:);
Time_vector = repmat(Time_mat(:),s2_cons,1);

% Check a few curves.
k = 5;
figure;plot(Time_mat,Product_mat(:,:,k),'--+');
xlabel('Time (s)  ')
ylabel('[Product]  ')
hline(0,'-y')

figure;plot(Product_vector,'-');
xlim([-50 3650]);ylim([-10 110])
%% Choice of parameters to refine.
% 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 = kcat.Value;
p3 = koff2.Value;
pin=[p1;p2;p3];

% We also need to 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.
% Here we use the nlinfit function from the Statistics Toolbox. This is the
% best choice if the toolbox is available.

% options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','fair');
options=statset('TolX',1e-6,'TolFun',1e-6,'Display','iter','RobustWgtFun','','UseParallel',true);
[FP1,R1,J1,Cov1,MSE1] = ...
    nlinfit(Time_vector,Product_vector,'bi_uni_fit',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
[fcurve1,delta1] = nlpredci('bi_uni_fit',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);

FP1

%% Regeneration of the calculated individual progress curves.

calc_product_mat = reshape(fcurve1,ntimes,s1_cons,s2_cons);

% 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.
close all
k = 10
figure;plot(Time_mat,Product_mat(:,:,k),'+');
hold on
colors = get(gca,'ColorOrder');
plot(Time_mat,calc_product_mat(:,:,k));

%%
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,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(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);
kcat_fit = FP1(2);
koff2_fit = FP1(3);
Vmax = kcat_fit;
Kd_S1 = koff1_fit/kon1.Value;
Kd_S2 = koff2_fit/kon2.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)    '];
Result = [string1;string2;string3];

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

%% Optional trimming of time points.
% ntimes = ntimes_bk;
% Time_points = Time_points_bk;
% Here we reset the time vector if we want to discard some time points.
ntimes_bk = ntimes;
Time_points_bk = Time_points;
skip_times = 0;
ntimes = ntimes - skip_times;
Time_points = Time_points(1:ntimes);
set(configset_m1.SolverOptions, 'OutputTimes',Time_points);

%% Parameter determination with a random subset of curves
global Small_conc_mat ncurves

% all_S1S2 = nchoosek([1:10],2);
% all_S1S2 = [all_S1S2;all_S1S2(:,[2 1])];
% npicks = size(all_S1S2,1);

%% Random selection
% close all
% select_ind = randi(npicks,3*ncurves,1);
% select_ind = unique(select_ind,'stable');
% pick_S1S2 = all_S1S2(select_ind(1:ncurves),:)
% pick_S1S2 = sortrows(pick_S1S2)
% pick_S1 = pick_S1S2(:,1);
% pick_S2 = pick_S1S2(:,2);

%% Rational selection
close all
ncurves = 10;

% pick_S1 = [2 7 1 1 2 2 3 3 4  4 10 5 6  6 8  7  9 10]';
% pick_S2 = [1 2 3 7 4 8 5 9 6 10  7 8 8 10 9 10 10 10]';
pick_S1 = [1 2 3  4 10 6  6 8  7  9 ]';
pick_S2 = [3 8 9 10  7 8 10 9 10 10 ]';

pick_S1S2 = [pick_S1 pick_S2]

Small_product_mat = zeros(ntimes,ncurves);
for i = 1:ncurves
Small_product_mat(:,i) = Product_mat(1:ntimes,pick_S1(i),pick_S2(i));
end

Small_conc_mat = [s1_conc_vec(pick_S1);s2_conc_vec(pick_S2)];
Small_time_mat = repmat(Time_mat(1:ntimes,1),1,ncurves);
figure;plot(Small_time_mat,Small_product_mat,'-o')

Small_product_vector = Small_product_mat(:);
Small_time_vector = Small_time_mat(:);

figure;plot(Small_product_vector,'--r');hold on
plot(Small_product_vector,'xb')

% Choice of parameters to refine for the global fit.
p1 = 50;
p2 = 0.1;
p3 = 200;
pin=[p1;p2;p3];

nobs = length(Small_time_vector);
nvars = length(pin);

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

% Solution with lsqcurvefit
options = ...
optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-6,'TolX',1e-6);
[FP2,sos2,R2,~,~,~,J2] = ...
lsqcurvefit(@small_bi_uni_fit,pin,Small_time_vector,Small_product_vector,...
[],[],options); 

MSE2 = sos2/(nobs-nvars);
RMSE2 = sqrt(MSE2)

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

J2 = full(J2);
Cov2 = inv(J2'*J2).*MSE2;

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

% Solution with nlinfit
% options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','bisquare');
% options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','');
% [FP2,R2,J2,Cov2,MSE2] = ...
%     nlinfit(Small_time_vector,Small_product_vector,'small_bi_uni_fit',pin,options);
% [Corr2,sigma2] = corrcov(Cov2);
% sos2 = R2'*R2

% The following works for both nlinfit and lsqcurvefit.
[fcurve2,delta2] = nlpredci('small_bi_uni_fit',Small_time_vector,FP2,R2,'Covar',Cov2);
FP2ci = nlparci(FP2,R2,'covar',Cov2);
FP2

%% Refinement with rational selection including substrates concentrations
close all
ncurves = 10;

% pick_S1 = [2 7 1 1 2 2 3 3 4  4 10 5 6  6 8  7  9 10]';
% pick_S2 = [1 2 3 7 4 8 5 9 6 10  7 8 8 10 9 10 10 10]';
pick_S1 = [1 2 3  4 10 6  6 8  7  9 ]';
pick_S2 = [3 8 9 10  7 8 10 9 10 10 ]';

pick_S1S2 = [pick_S1 pick_S2]

Small_product_mat = zeros(ntimes,ncurves);
for i = 1:ncurves
Small_product_mat(:,i) = Product_mat(1:ntimes,pick_S1(i),pick_S2(i));
end

Small_conc_mat = [s1_conc_vec(pick_S1);s2_conc_vec(pick_S2)];
Small_time_mat = repmat(Time_mat(1:ntimes,1),1,ncurves);
figure;plot(Small_time_mat,Small_product_mat,'-o')

Small_product_vector = Small_product_mat(:);
Small_time_vector = Small_time_mat(:);

figure;plot(Small_product_vector,'--r');hold on
plot(Small_product_vector,'xb')

% Choice of parameters to refine for the global fit.
% p1 = 50;
% p2 = 0.1;
% p3 = 200;
% pin=[p1;p2;p3];
% pin=[pin ; s1_conc_vec(pick_S1)' ; s2_conc_vec(pick_S2)'];
pin=[FP2 ; s1_conc_vec(pick_S1)' ; s2_conc_vec(pick_S2)'];
lpin = pin - 0.2*pin;
upin = pin + 0.2*pin;

nobs = length(Small_time_vector);
nvars = length(pin);

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

% Solution with lsqcurvefit
options = ...
optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-6,'TolX',1e-6);
% [FP3,sos3,R3,~,~,~,J3] = ...
% lsqcurvefit(@small_bi_uni_fit_2,pin,Small_time_vector,Small_product_vector,...
% [],[],options); 
[FP3,sos3,R3,~,~,~,J3] = ...
lsqcurvefit(@small_bi_uni_fit_2,pin,Small_time_vector,Small_product_vector,...
lpin,upin,options); 

MSE3 = sos3/(nobs-nvars);
RMSE3 = sqrt(MSE3)

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

J3 = full(J3);
Cov3 = inv(J3'*J3).*MSE3;

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

% Solution with nlinfit
% options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','');
% [FP3,R3,J3,Cov3,MSE3] = ...
%     nlinfit(Small_time_vector,Small_product_vector,'small_bi_uni_fit_2',pin,options);
% [Corr3,sigma3] = corrcov(Cov3);
% sos3 = R3'*R3

% The following works for both nlinfit and lsqcurvefit.
[fcurve3,delta3] = nlpredci('small_bi_uni_fit_2',Small_time_vector,FP3,R3,'Covar',Cov3);
FP3ci = nlparci(FP3,R3,'covar',Cov3)
FP3

%% Regeneration of the calculated individual progress curves.
small_calc_product_mat = reshape(fcurve3,ntimes,ncurves);

%%
% Check individual curves.
figure;plot(Small_time_mat,Small_product_mat,'+');
hold on
colors = get(gca,'ColorOrder');
plot(Small_time_mat,small_calc_product_mat);
xlabel('Time (s) ')
ylabel('[Product] ')

% Check concatenated vector
figure;plot(Small_product_vector,'+');
hold on
plot(fcurve3,'-r');
xlabel('TimeVector ')
ylabel('[Product] ')
% xlim([-10 370])
% ylim([-0.1*max(fcurve2) 1.1*max(fcurve2)])


%% Parameter determination with a distribution of subsets of 10 curves.
% % global Small_conc_mat ncurves
% % ncurves = 10;
% % ntrial = 400;
% % % FP2_distr = zeros(3,ntrial);
% % pick_S1S2_distr = zeros(ncurves,2,ntrial);
% % 
% % all_S1S2 = nchoosek([1:10],2);
% % all_S1S2 = [all_S1S2;all_S1S2(:,[2 1])]
% % npicks = size(all_S1S2,1);
% % 
% % %% Determine the distribution of the solution with 'ntrial' runs:
% % 
% % for idistr = 1:ntrial
% %     
% % % Here we select randomly 10 curves
% % close all
% % select_ind = randi(npicks,3*ncurves,1);
% % select_ind = unique(select_ind,'stable');
% % pick_S1S2 = all_S1S2(select_ind(1:ncurves),:)
% % pick_S1S2_distr(:,:,idistr) = pick_S1S2;
% % pick_S1 = pick_S1S2(:,1);
% % pick_S2 = pick_S1S2(:,2);
% % 
% % Small_product_mat = zeros(ntimes,ncurves);
% % for i = 1:ncurves
% % Small_product_mat(:,i) = Product_mat(1:ntimes,pick_S1(i),pick_S2(i));
% % end
% % 
% % Small_conc_mat = [s1_conc_vec(pick_S1);s2_conc_vec(pick_S2)];
% % Small_time_mat = repmat(Time_mat(1:ntimes,1),1,ncurves);
% % % figure;plot(Small_time_mat,Small_product_mat,'-o')
% % 
% % Small_product_vector = Small_product_mat(:);
% % Small_time_vector = Small_time_mat(:);
% % 
% % figure;plot(Small_product_vector,'--r');hold on
% % plot(Small_product_vector,'xb')
% % 
% % % Choice of parameters to refine for the global fit.
% % p1 = 50;
% % p2 = 0.1;
% % p3 = 200;
% % pin=[p1;p2;p3];
% % 
% % nobs = length(Small_time_vector);
% % nvars = length(pin);
% % 
% % % We reinitialize the initial value for the enzyme concentration and the
% % % cell volume.
% % E.InitialAmount = E_init;
% % C.capacity = C_init;
% % 
% % % Solution with lsqcurvefit
% % options = ...
% % optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-8,'TolX',1e-8);
% % [FP2,sos2,R2,~,~,~,J2] = ...
% % lsqcurvefit(@small_bi_uni_fit,pin,Small_time_vector,Small_product_vector,...
% % [],[],options); 
% % 
% % MSE2 = sos2/(nobs-nvars);
% % RMSE2 = sqrt(MSE2)
% % 
% % % 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).
% % 
% % J2 = full(J2);
% % cov2 = inv(J2'*J2).*MSE2;
% % 
% % % 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] 
% % 
% % % Solution with nlinfit
% % % options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','bisquare');
% % % options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','');
% % % [FP2,R2,J2,Cov2,MSE2] = ...
% % %     nlinfit(Small_time_vector,Small_product_vector,'small_bi_uni_fit',pin,options);
% % % [Corr2,sigma2] = corrcov(Cov2);
% % % sos2 = R2'*R2
% % 
% % % The following works for both nlinfit and lsqcurvefit.
% % [fcurve2,delta2] = nlpredci('small_bi_uni_fit',Small_time_vector,FP2,R2,'Covar',Cov2);
% % FP2ci = nlparci(FP2,R2,'covar',Cov2);
% % FP2
% % FP2_distr(:,idistr) = FP2;
% % 
% % end
% % 
% % %% Analyze the distribution
% % figure; hist(FP2_distr(1,:),30)
% % figure; hist(FP2_distr(2,:),30)
% % figure; hist(FP2_distr(3,:),30)
% % 
% % muld = 0.03;
% % p1d = muld*200;
% % p2d = muld*0.05;
% % p3d = muld*50;
% % 
% % p1_l = 200 - p1d;
% % p1_u = 200 + p1d;
% % p2_l = 0.05 - p2d;
% % p2_u = 0.05 + p2d;
% % p3_l = 50 - p3d;
% % p3_u = 50 + p3d;
% % good_FP2 = find(FP2_distr(1,:)>p1_l & FP2_distr(1,:)<p1_u & ...
% %                 FP2_distr(2,:)>p2_l & FP2_distr(2,:)<p2_u & ...
% %                 FP2_distr(3,:)>p3_l & FP2_distr(3,:)<p3_u);
% % 
% % for i = 1:length(good_FP2)
% %     for j = 1:length(good_FP2)
% %     good_FP2_score(i,j) = size(intersect(pick_S1S2_distr(:,:,i),...
% %         pick_S1S2_distr(:,:,j),'rows'),1);
% %     end
% % end
% % 
% % FP2_best = FP2_distr(:,good_FP2);
% % pick_S1S2_best = pick_S1S2_distr(:,:,good_FP2);
% % 
% % %% Repeat fit with best combination of S1 and S2
% % best_pick = 22;
% % pick_S1S2 = pick_S1S2_best(:,:,best_pick)
% % 
% % %%
% % % Alternatively load the combinations from a previous distribution search
% % load('TOOLBOXES/ENZYME_KINETICS/Bi_uni_Michaelis_Menten_Global_Fit','pick_S1S2')
% % 
% % pick_S1 = pick_S1S2(:,1);
% % pick_S2 = pick_S1S2(:,2);
% % 
% % Small_product_mat = zeros(ntimes,ncurves);
% % for i = 1:ncurves
% % Small_product_mat(:,i) = Product_mat(1:ntimes,pick_S1(i),pick_S2(i));
% % end
% % 
% % Small_conc_mat = [s1_conc_vec(pick_S1);s2_conc_vec(pick_S2)];
% % Small_time_mat = repmat(Time_mat(1:ntimes,1),1,ncurves);
% % figure;plot(Small_time_mat,Small_product_mat,'-o')
% % 
% % Small_product_vector = Small_product_mat(:);
% % Small_time_vector = Small_time_mat(:);
% % 
% % figure;plot(Small_product_vector,'--r');hold on
% % plot(Small_product_vector,'xb')
% % 
% % p1 = 50;
% % p2 = 0.1;
% % p3 = 200;
% % pin=[p1;p2;p3];
% % 
% % nobs = length(Small_time_vector);
% % nvars = length(pin);
% % 
% % E.InitialAmount = E_init;
% % C.capacity = C_init;
% % 
% % % Solution with lsqcurvefit
% % options = ...
% % optimoptions('lsqcurvefit','Display','iter','FinDiffType','central','TolFun',1e-8,'TolX',1e-8);
% % [FP2,sos2,R2,~,~,~,J2] = ...
% % lsqcurvefit(@small_bi_uni_fit,pin,Small_time_vector,Small_product_vector,...
% % [],[],options); 
% % 
% % MSE2 = sos2/(nobs-nvars);
% % RMSE2 = sqrt(MSE2)
% % J2 = full(J2);
% % cov2 = inv(J2'*J2).*MSE2;
% % 
% % % [Corr,sigma] = corrcov(cov2);
% % % FP2ci = [FP2-1.96*sigma FP2+1.96*sigma] 
% % 
% % % Further refinement with nlinfit
% % pin = FP2;
% % options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','');
% % [FP2,R2,J2,Cov2,MSE2] = ...
% %     nlinfit(Small_time_vector,Small_product_vector,'small_bi_uni_fit',pin,options);
% % [Corr2,sigma2] = corrcov(Cov2);
% % sos2 = R2'*R2
% % 
% % % The following works for both nlinfit and lsqcurvefit.
% % [fcurve2,delta2] = nlpredci('small_bi_uni_fit',Small_time_vector,FP2,R2,'Covar',Cov2);
% % FP2ci = nlparci(FP2,R2,'covar',Cov2);
% % FP2
% % 
% % %% Regeneration of the calculated individual progress curves.
% % small_calc_product_mat = reshape(fcurve2,ntimes,ncurves);
% % 
% % %%
% % % Check individual curves.
% % figure;plot(Small_time_mat,Small_product_mat,'+');
% % hold on
% % colors = get(gca,'ColorOrder');
% % plot(Small_time_mat,small_calc_product_mat);
% % xlabel('Time (s) ')
% % ylabel('[Product] ')
% % 
% % % Check concatenated vector
% % figure;plot(Small_product_vector,'+');
% % hold on
% % plot(fcurve2,'-r');
% % xlabel('TimeVector ')
% % ylabel('[Product] ')
% % xlim([-10 370])
% % ylim([-0.1*max(fcurve2) 1.1*max(fcurve2)])
% % 
% % %%
% % % We can rerun with increased weights for the regions that are not fit well. 
% % W = ones(360,1);
% % W([206:216 279:288]) = 10*W([206:216 279:288]);
% % [FP3,R3,J3,Cov3,MSE3] = ...
% %     nlinfit(Small_time_vector,Small_product_vector,'small_bi_uni_fit',pin,options,...
% %     'Weights',W);
% % [Corr3,sigma3] = corrcov(Cov3);
% % sos3 = R3'*R3
% % [fcurve3,delta3] = nlpredci('small_bi_uni_fit',Small_time_vector,FP3,R3,'Covar',Cov3);
% % FP3ci = nlparci(FP3,R3,'covar',Cov3);
% 
% 
% figure;plot(Small_product_vector,'+');
% hold on
% plot(fcurve2,'-r');
% plot(fcurve3,'-c');
% plot(W,'--g');
