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

%% Michaelis Menten Kinetics: Simulation and Analysis.

% Add to the path the enhanced scripts for drawing lines at specified
% coordinates.

global m1

% Load the project
sbioloadproject('Simple_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);

% Now we run a small simulation with the existing parameters. First, let's 
% get the configuration of the simulation.
configset_m1 = getconfigset(m1);
get(configset_m1)

Simple_MM_Kinetics = sbiosimulate(m1);

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

% We can plot the simulation
sbioplot(Simple_MM_Kinetics)

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

%%
close all

% Now we can change the simulation time.
set(configset_m1, 'StopTime', 4000);
set(configset_m1,'SolverType','ode15s');
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.0e-9);

% We can also configure the solver so that states are logged only at
% specific times.
global Time_points ntimes
Time_points = [0 1 10 100:100:4000]';
ntimes = size(Time_points,1);
set(configset_m1.SolverOptions, 'OutputTimes',Time_points);
get(configset_m1)

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

% We can plot the  new simulation
% 
sbioplot(Simple_MM_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.
global C E S ES P kon koff kcat
C = sbioselect(m1,'Name','Cell');
S = sbioselect(m1,'Name','S');
E = sbioselect(m1,'Name','E');
ES = sbioselect(m1,'Name','ES');
P = sbioselect(m1,'Name','P');
kon = sbioselect(m1,'Name','kon');
koff = sbioselect(m1,'Name','koff');
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;

% Check the variable values
% S
% E
% kon
% koff
% kcat

%%
close all

% Now we can set different values. We notice that Kd = koff/kon = 200/10 = 20
% microM, so we can try a value of the substrate equal to 1/10 the Kd and a
% value equal to twice the Kd. It is important to remember here that Km is
% different from Kd as Km = (koff + kcat)/kon = (200+0.05)/10 = 20.005. 
%
% In order to examine the effect of different concentrations of the
% substrate it is convenient to write a 'loop'. First we define a vector
% with all the concentrations we want to explore spaced logarithmically. It
% is a good idea to cover a range from 1/20 to ~50 times the Kd value
% in order to completely saturate the enzyme with substrate; we will go
% from 1 µM to 1 mM.

global ncons conc_vec
ncons = 20;
conc_vec = logspace(log10(1),log10(1000),ncons);
openvar conc_vec
Initial_Rate = zeros(ncons,1);
Product_mat = nan(ntimes,ncons);
Time_mat = nan(ntimes,ncons);
String_mat = cell(ncons,1);

for i = 1:ncons
    set(S,'InitialAmount',conc_vec(i));
    Simple_MM_Kinetics = sbiosimulate(m1);
    % Find the product concentration at 1 second.
    Product = Simple_MM_Kinetics.Data(2,4);
    Time = Simple_MM_Kinetics.Time(2);
    % Fill the 2D array with times and product values for each run.
    Product_mat(:,i) = Simple_MM_Kinetics.Data(:,4);
    Time_mat(:,i) = Simple_MM_Kinetics.Time(:,1);
    Initial_Rate(i) = Product/Time;
    String_mat{i,1} = [num2str(conc_vec(i),'%6.2f') ' µM']; 
end

%%
SIMPLE_MM = figure;
    set(SIMPLE_MM,'Units','normalized','Position',[0.4 0.0 0.6 1.0],...
        'Name','Simple Michaelis Menten Kinetics');clf
    subplot1 = subplot(3,2,1,'Parent',figure(gcf));
    box(subplot1,'on');
    grid(subplot1,'on');
    hold(subplot1,'all');
    
plot(Time_mat(:,2:2:20),Product_mat(:,2:2:20));
legend(gca,String_mat(2:2:20),'Location','BestOutside')
xlabel('Time (sec)');
ylabel('[P] (µM)');


% Now we can choose various ways of analyzing the data. We start by
% determining the initial velocity of the reaction for each initial
% concentration of the substrate. To measure this initial velocity we look
% at the reaction rate only during the 1st second.

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

% For example we can look at the reaction in the presence of 26.4 µM substrate.    
plot(Time_mat(:,10),Product_mat(:,10));
line_y = Initial_Rate(10)*1000;
line([0 1000],[0 line_y],'LineStyle','--','Color','r');
legend({String_mat{10} 'Vo (µM/sec)'},'Location','Best')
xlabel('Time (sec)');
ylabel('[P] (µM)');

% Then we plot all the initial velocities vs the corresponding concentrations 
% of substrate. 

    subplot3 = subplot(3,2,3,'Parent',figure(gcf));
    box(subplot3,'on');
    grid(subplot3,'on');
    hold(subplot3,'all');
    
plot(conc_vec,Initial_Rate,'-ob',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','c',...
             'MarkerSize',6);
%legend(gca,'show','Location','BestOutside')
xlabel('[S] (µM)');
ylabel('Vo (µM/sec)');

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

LB_x = 1./conc_vec;
LB_y = 1./Initial_Rate;
plot(LB_x,LB_y,'ob',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','c',...
             'MarkerSize',6);

xlim([-max(LB_x)/5,max(LB_x)]);
hold on

f = fittype('a*x + b');
[LB,GOF] = fit(LB_x',LB_y,f);

% Check the fit quality:
Rsquare2 = GOF.rsquare;

plot(LB,'-r');
legend('LB','Fit','Location','Best');
xlabel('1/[S]');
ylabel('1/V');
vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});
LB_params = coeffvalues(LB);

Vmax = 1/LB_params(2);
Km = LB_params(1)/LB_params(2);
% text
R2string = ['Vmax = ' num2str(Vmax,'%6.4f'),' (µM/sec) ';...
            'Km = ' num2str(Km,'%6.2f'),    ' (µM)        '];
annotation(gcf,'textbox',...
    [0.64 0.4 0.2 0.2],...
    'String',R2string,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

%%
% Of course this is what we would obtain with perfect data. We can also 
% simulate the presence of errors in the set-up, which would result in the 
% presence of 'noise' in the data.

% close all

ncons = 20;
conc_vec = logspace(log10(1),log10(1000),ncons);
% openvar conc_vec
Initial_Rate = zeros(ncons,1);
Product_mat = nan(ntimes,ncons);
Time_mat = nan(ntimes,ncons);
Init_mat = nan(ncons,3);

for i = 1:ncons
    
    error = (rand(3,1)-0.5);
    error(1) = 1 + error(1)*0.10;
    error(2) = 1 + error(2)*0.05;
    error(3) = 1 + error(3)*0.05;
    set(C,'Capacity',C_init*error(1));
    set(E,'InitialAmount',E_init*error(2));
    set(S,'InitialAmount',conc_vec(i)*error(3));
    Init_mat(i,1) = C.Capacity;
    Init_mat(i,2) = E.InitialAmount;
    Init_mat(i,3) = S.InitialAmount;
    
    Simple_MM_Kinetics = sbiosimulate(m1);
    % Find the product concentration at 1 second.
    Product = Simple_MM_Kinetics.Data(2,4);
    Time = Simple_MM_Kinetics.Time(2);
    % Fill the 2D array with times and product values for each run.
    Product_mat(:,i) = Simple_MM_Kinetics.Data(:,4);
    Time_mat(:,i) = Simple_MM_Kinetics.Time(:,1);
    Initial_Rate(i) = Product/Time;
    String_mat{i,1} = [num2str(conc_vec(i),'%6.2f') ' µM']; 
end

% Then we plot all the initial velocities vs the corresponding 
% concentrations of substrate. 

    subplot5 = subplot(3,2,5,'Parent',figure(gcf));
    box(subplot5,'on');
    grid(subplot5,'on');
    hold(subplot5,'all');
    
plot(conc_vec,Initial_Rate,'-ob',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','c',...
             'MarkerSize',6);
%legend(gca,'show','Location','BestOutside')
xlabel('[S] (µM)');
ylabel('Vo (µM/sec)');

% figure
% semilogx(conc_vec,Initial_Rate);

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

LB_x = 1./conc_vec;
LB_y = 1./Initial_Rate;
% plot(LB_x,LB_y,'ob');
plot(LB_x,LB_y,'ob',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor','c',...
             'MarkerSize',6);

xlim([-max(LB_x)/5,max(LB_x)]);
hold on

f = fittype('a*x + b');
[LB,GOF] = fit(LB_x',LB_y,f);

% Check the fit quality:
Rsquare = GOF.rsquare;

plot(LB,'-r');
legend('LB','Fit','Location','Best');
xlabel('1/[S]');
ylabel('1/V');
vline(0,{'g-', 'LineWidth', 0.5});
hline(0,{'g-', 'LineWidth', 0.5});
LB_params = coeffvalues(LB);

Vmax = 1/LB_params(2);
Km = LB_params(1)/LB_params(2);
% text
R2string = ['Vmax = ' num2str(Vmax,'%6.4f'),' (µM/sec) ';...
            'Km = ' num2str(Km,'%6.2f'),    ' (µM)        '];
annotation(gcf,'textbox',...
    [0.64 0.1 0.2 0.2],...
    'String',R2string,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

%%
% However, rather than using a LB plot we could fit directly the
% plot of initial rate vs the amount of substrate (L). This plot is a
% hyperbola of the type y = a*x/(b+x) given by the Briggs-Haldane 
% relationship: v = Vmax[S]/(Km + [S])

X = conc_vec';
Y = Initial_Rate;

SIMPLE_MM_2 = figure;
    set(SIMPLE_MM_2,'Units','normalized','Position',[0.4 0.2 0.6 0.8],...
        'Name','Simple Michaelis Menten Kinetics');clf

subplot1 = subplot(2,2,1,'Parent',figure(gcf));
box(subplot1,'on');
grid(subplot1,'on');
hold(subplot1,'all');
    
plot(X,Y,'o',...
             'MarkerEdgeColor','b',...
             'MarkerFaceColor','c',...
             'MarkerSize',5);
f = fittype('a*x/(b + x)');
[Hyperb,GOF] = fit(X,Y,f,'StartPoint',[0.05 20]);
plot(Hyperb,'-r');
legend('Vo','Fit','Location','Best');
xlabel('[S]');
ylabel('Vo');

% Retrieve the fit parameters.
% Hyperb
% GOF
% or simply:
Hyperb_params = coeffvalues(Hyperb);
Vmax = Hyperb_params(1);
Km = Hyperb_params(2);

string1 = ['Vmax = ' num2str(Vmax,'%6.4f\n') ' (µM/sec)'];
string2 = ['Km = ' num2str(Km,'%6.4f\n') ' (µM)     '];
string = [string1;string2];
% string2 = num2str(Kd,5);
% Create textbox
annotation(gcf,'textbox',...
    [0.3 0.6 0.1 0.25],...
    'String',string,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);

%% Global fit
% 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. The vectors 
% 'Time_vector_ind' and 'Product_vector_ind' are not used here, but are
% derived anyway for consistent use of the function 'uni_uni_fit_all'

global Time_vector_ind Product_vector_ind

Product_vector = zeros(ntimes*ncons,1);
Time_vector = zeros(ntimes*ncons,1);
for k = 1:ncons
    for j = 1:ntimes
    Product_vector(j+(k-1)*ntimes,1) = Product_mat(j,k);
    Time_vector(j+(k-1)*ntimes,1) = Time_mat(j,k);
    end;
end;
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);

%% Choice of parameters to refine (part I)
% We need to set some initial value for the refinement of parameters: we
% will refine kon, koff, and kcat. For example, we can take the initial 
% values from the previous fit
p1 = 10; % Standard value for 2nd order rate constant
p2 = Hyperb_params(2)*p1;
p3 = Hyperb_params(1);
pin=[p1;p2;p3];
% MIN=[1;20;0.005];
% MAX=[100;2000;0.5];

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

%--------------------------------------------------------------------------
% Here we use the nlinfit function from the Statistics Toolbox. This is the
% best choice if the toolbox is available.
options=optimset('TolX',1e-10);
[FP1,R1,J1,Cov1,MSE1] = nlinfit(Time_vector,Product_vector,'uni_uni_fit',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
N = length(Product_mat(:));
N2 = length(pin);
MSE2 = sos1/(N-N2);

% Cov1 is the covariance matrix of the predicted parameters (FP1): on the
% diagonal there are the squared errors of the 3 paramaters and their
% products represent the off-diagonal terms. Check:
test = sqrt(diag(Cov1));
stest = test*test';

% The Jacobian is the matrix of partial derivatives of the function of 3 
% variables the represents the 3 unknown parameters, in the immediate 
% neighbourood of its minimum. The covariance matrix is actually calculated 
% as the inverse of the J1'*J1 product, scaled by the mean squared error.

J1product = J1'*J1;

% Here we calculate the scaled covariance matrix from the Jacobian
Cov2 = inv(J1product).*MSE2;

[Corr,sigma] = corrcov(full(Cov2));

% We can also calculate the partial correlation matrix from the inverse of
% the covariance matrix:
% cov_inv = inv(Cov2);
% or
cov_inv = J1product./MSE2;

    rho = zeros(N2,N2);

    for i = 1:N2
        for j = i:N2
            rho(i,j) = cov_inv(i,j)/sqrt(cov_inv(i,i)*cov_inv(j,j));
            rho(j,i) = rho(i,j);
        end
    end
    
% Here we calculate the predicted progress curves and the margin of
% confidence.
[fcurve1,delta1] = nlpredci('uni_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);
%--------------------------------------------------------------------------

theo = zeros(ntimes,ncons);
for k=1:ncons
    for j=1:ntimes
    theo(j,k)=fcurve1(j+(k-1)*ntimes,1);
    end;
end;

%%
Vmax = FP1(3);
Km = (FP1(2)+FP1(3))/FP1(1);

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

plot(Time_mat(1:ntimes,:),Product_mat(1:ntimes,:),'+');
get(gca,'ColorOrder');hold on
plot(Time_mat(1:ntimes,:),theo(1:ntimes,:),'-');
xlabel('Time (sec)');
ylabel('[P] (µM)');

string1 = ['Vmax = ' num2str(Vmax,'%6.4f\n') ' (µM/sec)'];
string2 = ['Km = ' num2str(Km,'%6.4f\n') ' (µM)     '];
string = [string1;string2];
% string2 = num2str(Kd,5);
% Create textbox
annotation(gcf,'textbox',...
    [0.6 0.6 0.1 0.25],...
    'String',string,...
    'FitBoxToText','on','EdgeColor','r','LineWidth',0.2);


%% Plot partial correlation between parameters

subplot3 = subplot(2,2,3,'Parent',figure(gcf));
box(subplot3,'on');
grid(subplot3,'on');
hold(subplot3,'all');
    
imagesc(Corr)
set(gca,'YDir','normal');
colorbar
vline(1.5,{'c-', 'LineWidth', 0.5});
vline(2.5,{'c-', 'LineWidth', 0.5});
hline(1.5,{'c-', 'LineWidth', 0.5});
hline(2.5,{'c-', 'LineWidth', 0.5});
xlabel('Variable Number');
ylabel('Variable Number');
title('Correlation');

subplot4 = subplot(2,2,4,'Parent',figure(gcf));
box(subplot4,'on');
grid(subplot4,'on');
hold(subplot4,'all');
    
imagesc(rho)
set(gca,'YDir','normal');
colorbar
vline(1.5,{'c-', 'LineWidth', 0.5});
vline(2.5,{'c-', 'LineWidth', 0.5});
hline(1.5,{'c-', 'LineWidth', 0.5});
hline(2.5,{'c-', 'LineWidth', 0.5});
xlabel('Variable Number');
ylabel('Variable Number');
title('Partial Correlation');

%%
