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