% Copyright (c) 2013, 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.
%
%%
clear, clc, close all

%% Michaelis Menten Kinetics: Competitive Inhibition Simulation and Analysis.

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

global m1
global Time_points ntimes
global ncons conc_vec i_ncons i_conc_vec
global C C_init E E_init S I ES EI P kon koff kcat ki_on ki_off

% Load the project
sbioloadproject('Comp_Inhib_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.
Comp_Inhib_MM_Kinetics = sbiosimulate(m1);

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

% We can plot the simulation
% sbioplot(Comp_Inhib_MM_Kinetics)

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

%%
close all

% Then we extend the simulation time, we set the solver and request a
% smaller tolerance
set(configset_m1, 'StopTime', 13000);
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_points = [0 1 10 100:100:13000]';

ntimes = size(Time_points,1);
set(configset_m1.SolverOptions, 'OutputTimes',Time_points);

get(configset_m1)

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

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

% Here we change the value of ki_off for the inhibitor.
ki_off.Value = 50;

% In order to examine the effect of different concentrations of the
% substrate and the inhibitpr 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 ~20
% times the Kd and Ki values in order to completely saturate the enzyme with
% substrate and/or inhibitor:

ncons = 20;
i_ncons = 10;
conc_vec = logspace(log10(1),log10(400),ncons);
i_conc_vec = logspace(log10(1),log10(400),i_ncons);
i_conc_vec =[0 i_conc_vec];
i_ncons = length(i_conc_vec);
% openvar conc_vec
Initial_Rate = zeros(ncons,i_ncons);
Product_mat = nan(ntimes,ncons,i_ncons);
Time_mat = nan(ntimes,ncons);
String_mat = cell(i_ncons,1);
Init_mat = nan(ncons,4,i_ncons);

for k = 1:i_ncons
    % Here we just get the legend entries.
    String_mat{k,1} = [num2str(i_conc_vec(k),'%6.2f') ' µM']; 
    for i = 1:ncons
        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));
        set(S,'InitialAmount',conc_vec(i)*error(3));
        set(I,'InitialAmount',i_conc_vec(k)*error(4));
        Init_mat(i,1,k) = C.Capacity;
        Init_mat(i,2,k) = E.InitialAmount;
        Init_mat(i,3,k) = S.InitialAmount;
        Init_mat(i,4,k) = I.InitialAmount;
        
        Comp_Inhib_MM_Kinetics = sbiosimulate(m1);
        % Find the product concentration at 1 second.
        Product = Comp_Inhib_MM_Kinetics.Data(2,4);
    	Time = Comp_Inhib_MM_Kinetics.Time(2);
        % Fill the 2D array with times and product values for each run.
        Product_mat(:,i,k) = Comp_Inhib_MM_Kinetics.Data(:,4);
        Time_mat(:,i) = Comp_Inhib_MM_Kinetics.Time(:,1);
        Initial_Rate(i,k) = Product/Time;
    end
end

%%
COMP_INHIB_MM_1 = figure;
    set(COMP_INHIB_MM_1,'Units','normalized','Position',[0.0 0.4 1.0 0.6],...
        'Name','Competitive Inhibition 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(conc_vec,Initial_Rate(:,1:2:i_ncons));
%     legend({String_mat{1:2:i_ncons}},'Location','Best')
    plot(conc_vec,Initial_Rate);
%    legend({String_mat{1:2:i_ncons}},'Location','Best')
    legend(String_mat,'Location','Best')
    xlim([-10,1.1*max(conc_vec)]);
    ylim([0,1.1*max(Initial_Rate(:))]);
    xlabel('[S] (µ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');

Vmax = zeros(i_ncons,1);
Km = zeros(i_ncons,1);
R2 = zeros(i_ncons,1);
LB_x = 1./conc_vec;
LB_y = 1./Initial_Rate;
a = zeros(i_ncons,1);
b = zeros(i_ncons,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:i_ncons
    
% 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;

% We can also calculate directly using the closed form solution.
% [a(i),b(i),R2(i)] = linreg(LB_x,LB_y(:,i)');

Vmax(i) = 1/b(i);   % 1/fit_params(2);
Km(i) = a(i)/b(i);  % fit_params(1)/fit_params(2);

end

% for i = 1:2:i_ncons
for i = 1:i_ncons
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]);

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

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

String_mat2 = cell(i_ncons*2,1);
String_mat2(1:2:(i_ncons*2-1)) = String_mat;

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

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

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


% Here we derive the Km and Ki values from a replot of the apparent Km's at
% different concentrations of the inhibitor.

    subplot3 = subplot(1,3,3,'Parent',figure(gcf));
    box(subplot3,'on');
    grid(subplot3,'on');
    hold(subplot3,'all');    
    
[FITRESULT,GOF] = Robust_Linear_Least_Squares(i_conc_vec',Km,'bisquare');

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

low_x = -30;
high_x = 1.1*max(i_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:i_ncons
    plot(i_conc_vec(i),Km(i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
    end

% plot(i_conc_vec,Km,'o',...
%              'MarkerEdgeColor','k',...
%              'MarkerFaceColor','cyan',...
%              'MarkerSize',6);

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

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

xlim([low_x,1.1*high_x]);
% ylim([low_y,1.1*high_y]);
ylim([low_y,1.1*max(Km)]);

% legend({String_mat{1:2:i_ncons}},'Location','Best')
% Axes label
xlabel('[I] (µM)');
ylabel('Km (µM)');

% text
Vmax_mean = mean(Vmax);
Vmax_std = std(Vmax);
string1 = ['Vmax = ' num2str(Vmax_mean,'%6.2e'),' (µM/sec)  '];
string2 = ['Km =   ' num2str(Km_linreg,'%6.2e'),' (µM)      '];
string3 = ['Ki =   ' num2str(Ki_linreg,'%6.2e'),' (µM)      '];
Result = [string1;string2;string3];

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

%--------------------------------------------------------------------------
% Would we get the same result if we used our own least-square?
% By now we also know how to do a linear least-square with weighting
% accounting for the spread of the data. Previously we have developed the 
% least-square code to fit a hyperbola. It should quite easy to fit a 
% straight line to the data. As usual, first we write the
% equations:
% b + ax(1) = y(1)
% b + ax(2) = y(2)
% b + ax(3) = y(3)
% ...
% b + ax(m) = y(m)
% 
% % Then we represent them as a single matrix equation:
% [ 1 x(1)] [b]   [y(1)]
% | 1 x(2)|*[ ] = |y(2)]
% | 1 x(3)| [a]   |y(3)]
% |...    |       |... |
% [ 1 x(m)]       [y(m)]
% 
% In this case 'i_cons_vec' is 'y' and 'Km' is 'x'. Thus our matrix equation 
% is M*U = Km:
% 
% [ones(i_ncons,1) i_conc_vec'] * U = Km 
%  
% from which we derive the 'normal equation'
% M'*M*U = M'*Km, which we solve as:

M = [ones(i_ncons,1) i_conc_vec'];
U = inv(M'*M)*M'*Km;
U
% U value are quite different from those obtained using the robust
% function. However we know how to carry out a 'weighted least-square'. First
% we determine the covariance matrix for the values in 'Km':
cKm = Km - mean(Km);
sigma = cKm*cKm';

% COV = sigma;
% lambda = 0.00001;
%     imatrix = eye(i_ncons);
%     fmatrix = mean(diag(COV))*imatrix;
%     
%     sCOV = COV;
%     while lambda <= 1.0
%         try
%             sCOV = lambda*fmatrix + (1-lambda)*COV;    
%             test_chol = chol(sCOV);
%         end
%         if exist('test_chol','var')
%             break
%         end
%         lambda = lambda + delta;
%     end
% sigma = sCOV;

% We take only the diagonal
sigma = sigma .* eye(i_ncons);
inv_sigma = inv(sigma);
U = inv(M'*inv_sigma*M)*M'*inv_sigma*Km;
U

b_i2 = U(1);
a_i2 = U(2);

low_y = low_x*a_i2+b_i2;
high_y = high_x*a_i2+b_i2;

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

line(x,y,'Color','magenta','LineStyle','-');
%--------------------------------------------------------------------------

% Finally we derive values for some of the rate constants.
koff.Value = Km_linreg*kon.Value - Vmax_mean;
kcat.Value = Vmax_mean;
ki_off.Value = Ki_linreg*ki_on.Value;

%% 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 could
% be derived anyway for consistent use of the function 'sbiosim_fit_all'

% global Time_vector_ind Product_vector_ind

Product_vector = zeros(ntimes*ncons*i_ncons,1);
Time_vector = zeros(ntimes*ncons*i_ncons,1);
for i = 1:i_ncons
    for k = 1:ncons
        for j = 1:ntimes
        Product_vector(j+(k-1)*ntimes+(i-1)*ncons*ntimes,1) = Product_mat(j,k,i);
        Time_vector(j+(k-1)*ntimes+(i-1)*ncons*ntimes,1) = Time_mat(j,k);
        end;
    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.
% We need to set some initial value for the refinement of parameters: we
% will refine koff, kcat, and ki_off. For example, we can take the initial 
% values from the previous fit
p1 = koff.Value;
p2 = kcat.Value;
p3 = ki_off.Value;
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;

%% 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=optimset('TolX',1e-10);
[FP1,R1,J1,Cov1,MSE1] = ...
    nlinfit(Time_vector,Product_vector,'comp_inhib_fit',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
[fcurve1,delta1] = nlpredci('comp_inhib_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);

%% Regeneration of the calculated individual progress curves.
calc_product_mat = zeros(ntimes,ncons,i_ncons);
for i = 1:i_ncons
    for k=1:ncons
        for j=1:ntimes
        calc_product_mat(j,k,i)=fcurve1(j+(k-1)*ntimes+(i-1)*ncons*ntimes,1);
        end
    end
end

% Here we derive the calculated initial velocities
calc_initial_rate = zeros(ncons,i_ncons);
for k = 1:i_ncons
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));

%%
COMP_INHIB_MM_2 = figure;
    set(COMP_INHIB_MM_2,'Units','normalized','Position',[0.4 0.4 0.6 0.6],...
        'Name','Competitive Inhibition 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:i_ncons
plot(conc_vec,Initial_Rate(:,i),'o',...
             'MarkerEdgeColor','k',...
             'MarkerFaceColor',colors(i,:),...
             'MarkerSize',6);
hold on
plot(conc_vec,calc_initial_rate(:,i),'--',...
             'Color',colors(i,:));
         
end

    xlim([-10,1.1*max(conc_vec)]);
    ylim([-0.001,1.1*max(Initial_Rate(:))]);
    xlabel('[S] (µM)');
    ylabel('Vo (µM/sec)');

String_mat2 = cell(i_ncons*2,1);
String_mat2(1:2:(i_ncons*2-1)) = String_mat;

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

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

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

% text
koff_fit = FP1(1);
kcat_fit = FP1(2);
ki_off_fit = FP1(3);
Vmax = kcat_fit;
Km = (koff_fit + kcat_fit)/kon.Value;
Ki = ki_off_fit/ki_on.Value;
string1 = ['Vmax = ' num2str(Vmax,'%6.2e'),' (µM/sec)'];
string2 = ['Km =   ' num2str(Km,'%6.2e'),' (µM)    '];
string3 = ['Ki =   ' num2str(Ki,'%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});


