## Contents

- Two Substrates Bi-Uni Michaelis Menten Kinetics: Simulation and Analysis.
- Global fit
- Choice of parameters to refine.
- Non-linear global fit of all the progress curves.
- Regeneration of the calculated individual progress curves.
- Plot correlation between parameters
- Parameter determination with a subset of 10 curves.
- Global fit
- Choice of parameters to refine.
- Non-linear global fit of all the progress curves.

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

       SimBiology Compartment - Cell

       Compartment Components:
         Capacity:          1
         CapacityUnits:     milliliter
         Compartments:      0
         ConstantCapacity:  true
         Owner:
         Species:           6

       SimBiology Species Array

       Index:    Compartment:    Name:    InitialAmount:    InitialAmountUnits:
       1         Cell            S1       100               micromolarity
       2         Cell            E        1                 micromolarity
       3         Cell            ES1      0                 micromolarity
       4         Cell            P        0                 micromolarity
       5         Cell            S2       50                micromolarity
       6         Cell            ES2      0                 micromolarity

       SimBiology Parameter Array

       Index:    Name:    Value:    ValueUnits:
       1         koff1    200       1/second
       2         kon1     10        1/(micromolarity*second)
       3         kcat     0.05      1/second
       4         kon2     10        1/(micromolarity*second)
       5         koff2    50        1/second

       SimBiology Reaction Array

       Index:    Reaction:
       1         S1 + E <-> ES1
       2         ES2 -> P + E
       3         S2 + ES1 <-> ES2

    equations =

    ODEs:
    d(S1)/dt = 1/Cell*(-ReactionFlux1)
    d(E)/dt = 1/Cell*(-ReactionFlux1 + ReactionFlux2)
    d(ES1)/dt = 1/Cell*(ReactionFlux1 - ReactionFlux3)
    d(P)/dt = 1/Cell*(ReactionFlux2)
    d(S2)/dt = 1/Cell*(-ReactionFlux3)
    d(ES2)/dt = 1/Cell*(-ReactionFlux2 + ReactionFlux3)

    Fluxes:
    ReactionFlux1 = (kon1*S1*E)*Cell-(koff1*ES1)*Cell
    ReactionFlux2 = (kcat*ES2)*Cell
    ReactionFlux3 = (kon2*S2*ES1)*Cell-(koff2*ES2)*Cell

    Parameter Values:
    koff1 = 200 1/second
    kon1 = 10 1/(micromolarity*second)
    kcat = 0.05 1/second
    kon2 = 10 1/(micromolarity*second)
    koff2 = 50 1/second
    Cell = 1 milliliter

    Initial Conditions:
    S1 = 100 micromolarity
    E = 1 micromolarity
    ES1 = 0 micromolarity
    P = 0 micromolarity
    S2 = 50 micromolarity
    ES2 = 0 micromolarity

                            Active: 1
                    CompileOptions: [1x1 SimBiology.CompileOptions]
                              Name: 'default'
                             Notes: ''
                    RuntimeOptions: [1x1 SimBiology.RuntimeOptions]
        SensitivityAnalysisOptions: [1x1 SimBiology.SensitivityAnalysisOptions]
                     SolverOptions: [1x1 SimBiology.ODESolverOptions]
                        SolverType: 'ode15s'
                          StopTime: 10000
               MaximumNumberOfLogs: Inf
                  MaximumWallClock: Inf
                         TimeUnits: 'second'
                              Type: 'configset'

Now we run a small simulation with the existing parameters.

    Two_Substr_MM_Kinetics = sbiosimulate(m1);

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

    % We can plot the simulation
    sbioplot(Two_Substr_MM_Kinetics)

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

             Data: [39x6 double]
        DataCount: [1x1 struct]
         DataInfo: {6x1 cell}
        DataNames: {6x1 cell}
        ModelName: 'Two_Substrates_Michaelis_Menten'
             Name: ''
            Notes: ''
          RunInfo: [1x1 struct]
             Time: [39x1 double]
        TimeUnits: 'second'
         UserData: []

    close all

    % Then we set the simulation time and the solver tolerance
    stop_time = 6000;

    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.
    Two_Substr_MM_Kinetics = sbiosimulate(m1);

    % We can plot the  new simulation
    %
    sbioplot(Two_Substr_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');
    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;

    % Check the variable values
    % S
    % E
    % I
    % kon1
    % koff1
    % kon2
    % koff2
    % kcat

    % Let's see what happens if we change the value of koff2 for the second
    % substrate and its initial concentration.
    % ES2.InitialAmount = 50;
    % koff2.Value = 50;

    % We rerun the simulation:
    % Two_Substr_MM_Kinetics = sbiosimulate(m1);

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

                            Active: 1
                    CompileOptions: [1x1 SimBiology.CompileOptions]
                              Name: 'default'
                             Notes: ''
                    RuntimeOptions: [1x1 SimBiology.RuntimeOptions]
        SensitivityAnalysisOptions: [1x1 SimBiology.SensitivityAnalysisOptions]
                     SolverOptions: [1x1 SimBiology.ODESolverOptions]
                        SolverType: 'ode15s'
                          StopTime: 6000
               MaximumNumberOfLogs: Inf
                  MaximumWallClock: Inf
                         TimeUnits: 'second'
                              Type: 'configset'

    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;

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

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

    kcat.Value = mean([Vmax_S1 Vmax_S2])
    koff1.Value = Ki_S1*kon1.Value
    koff2.Value = Km_S2*kon2.Value - kcat.Value

       SimBiology Parameter Array

       Index:    Name:    Value:       ValueUnits:
       1         kcat     0.0520255    1/second

       SimBiology Parameter Array

       Index:    Name:    Value:    ValueUnits:
       1         koff1    182.66    1/second

       SimBiology Parameter Array

       Index:    Name:    Value:     ValueUnits:
       1         koff2    63.2039    1/second

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

    Product_vector = Product_mat(:);
    Time_vector = repmat(Time_mat(:),s2_cons,1);
    % 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 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];
    % 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=statset('TolX',1e-8,'Display','iter','RobustWgtFun','fair');

    options=statset('TolX',1e-8,'Display','iter','RobustWgtFun','');
    [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);

                                         Norm of         Norm of
       Iteration             SSE        Gradient           Step
      -----------------------------------------------------------
               0         163.896
               1         131.035         45.2721         13.0714
               2         130.987         9235.49         1.35223
               3         130.986         219.656        0.100145
               4         130.955         71.5803        0.689964
               5         130.951         33.3225        0.703919
               6         130.949         113.857        0.188867
               7         130.946         140.447        0.158554
               8         130.946         145.304       0.0181233
               9         130.945         146.917       0.0178173
              10         130.945         147.857       0.0180977
              11          130.92         189.172       0.0187523
              12         130.919         199.126       0.0148818
              13         130.918          194.53       0.0150349
              14         130.918         192.877       0.0140211
              15         130.918         192.075      0.00148284
              16         130.918         191.773      0.00145003

    Error using SBCompiler.SimulationObject/simulate

    Error in sbiosimulate (line 143)
    [t x] = simobj.simulate(mobj, cs, variants, doses);

    Error in bi_uni_fit (line 46)
            Sim_mat = sbiosimulate(m1);

    Error in nlinfit>@(b,x)w.*model(b,x) (line 274)
            modelw = @(b,x) w.*model(b,x);

    Error in nlinfit>LMfit (line 578)
                yfit = model(beta,X);

    Error in nlinfit (line 275)
            [beta,J,~,cause,fullr] = LMfit(X,yw, modelw,beta,options,verbose,maxiter);

    Error in Bi_Uni_Michaelis_Menten_Global_Fit (line 706)
    [FP1,R1,J1,Cov1,MSE1] = ...

## Regeneration of the calculated individual progress curves.

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

    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.

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

## Parameter determination with a subset of 10 curves.

Here we select randomly 10 curves

    global Small_conc_mat ncurves
    close all
    ncurves = 10;
    pick_S1S2 = randi(ncurves,100,2);
    pick_S1S2 = unique(pick_S1S2,'rows');
    pick_S1 = pick_S1S2(:,1);
    pick_S2 = pick_S1S2(:,2);

    select_ind_1 = pick_S1 ~= pick_S2;
    pick_S1 = pick_S1(select_ind_1);
    pick_S2 = pick_S2(select_ind_1);
    select_ind_2 = randi(sum(select_ind_1),ncurves,1);
    pick_S1 = pick_S1(select_ind_2);
    pick_S2 = pick_S2(select_ind_2);
    Small_product_mat = zeros(ntimes,ncurves);
    for i = 1:ncurves
    Small_product_mat(:,i) = Product_mat(:,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),1,ncurves);
    figure;plot(Small_time_mat,Small_product_mat,'-o')

## Global fit

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

## Choice of parameters to refine.

    p1 = 50;
    p2 = 0.1;
    p3 = 200;
    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','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
    [fcurve2,delta2] = nlpredci('small_bi_uni_fit',Small_time_vector,FP2,R2,'Covar',Cov2);
    FP2ci = nlparci(FP2,R2,'covar',Cov2);

    % Regeneration of the calculated individual progress curves.

    small_calc_product_mat = reshape(fcurve2,ntimes,ncurves);

Check a few curves.

    figure;plot(Small_time_mat,Small_product_mat,'+');
    hold on
    colors = get(gca,'ColorOrder');
    plot(Small_time_mat,small_calc_product_mat);

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

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