by Ben Noe
University of Colorado at Boulder
Updated April 29, 2012
In collaboration with Dr. David Bortz and other members of the Mathematical Biology Undergraduate Research Team, my
Change in Enzyme Concentration
research has focused on accurately modeling the kinetics of any simple substrate – enzyme reaction (in the example below, I used the hydrolysis of benzoyl-L-arginine ethyl ester by
trypsin). More specifically, I created a computer program with Matlab that demonstrates how the concentrations of each chemical species changes over time, shows the reaction rate over time, and determines when the reaction is in the standard quasi-steady-state. For the sake of completeness, the documented code for my program will follow a description of the principles and ideas involved.
We’ll start with a basic discussion of simple substrate-enzyme kinetics. In 1913, Leonor Michaelis and Maud Menten proposed the following chemical reaction to relate the substrate, enzyme, substrate-enzyme complex, and product.
and correspond to the arrows above, which represent the rates of each step in the reaction.
Using the Law of Mass Action, where the rate of a reaction is directly proportional to the product of the concentrations of
Change in Substrate Concentration
the reactants, we can derive four differential equations to model how the concentration of each chemical species changes.
In 1925, G.E. Briggs and J.B.S. Haldane defined the standard quasi-steady-state (sQSSA) as the part of the reaction where the concentration of the intermediate substrate-enzyme complex does not change appreciably. In mathematical terms, Failed to parse (lexing error): \frac{d[S]}{dt} ≈ 0
. From this approximation, we can use the enzyme conservation law,
and substitute to derive the velocity, or rate of product formation, of the substrate-enzyme reaction.
Where and
Armed with the above knowledge, I proceeded to create a computer program in Matlab to easily study the progression of the reaction over time. By utilizing a numerical solver such as ode15s, one can easily obtain solutions for the system of four differential equations described above by inputting the appropriate reaction rates and initial concentrations.
Change in Substrate-Enzyme Complex Concentration
Unfortunately, values for and are not readily available in many documented substrate-enzyme reactions. However, values for the Michaelis_Menten constant, , and the rate-limiting turnover number,, have been widely studied. Therefore, the next step in creating a useful simulation was to calculate and in terms of and so that the user could input readily available values.
When calculating the values of and , I derived the following solutions:
However, regarding the differential equation solver ode15s in Matlab, the solutions simply had too many variables with the concentrations of the chemical species. Therefore, it extended beyond the scope of the function. Thus, we decided to approximate
and take
as zero. We successfully tested many reactions and found that this
Change in Product Concentration
approximation was valid with nearly identical solutions since
is usually negligible in comparison to the other rates.
In addition, I also created another similar function that can use and values if those values are available or if the above assumption does not hold for a given reaction, thus providing a slightly more accurate solution to the given problem.
The final step in creating a useful program was to develop a test to see when the reaction can be considered in the quasi-steady-state. In previous work on Michaelis-Menten kinetics, mathematicians and scientists have focused on under what conditions the steady state assumption would be valid, but have not to our knowledge developed a test with a specific
epsilon value to see when the reaction was indeed in the steady state. But first we’ll look at a review of what others have worked on.
In 1989, Lee and Slemrod 1 proposed the conditions that with . They also proposed that the substrate-enzyme complex could be modeled by the equation during the QSSA nullcline. Furthermore, the stated the condition that , where is the duration of the initial fast transient and is the time for significant change in [S] in the post transient period. Lastly, they calculated under what conditions the errors in the quasi-steady state approximation would be small (for example, ).
In 1996, Borghans et al. 2 proposed the condition and also introduced the following three steady states:
sQSSA Failed to parse (lexing error): \frac{d[C]}{dt} ≈ 0
rQSSA Failed to parse (lexing error): \frac{d[S]}{dt} ≈ 0
tQSSA Failed to parse (lexing error): \frac{d[S_{t}]}{dt} ≈ 0
where ,math> S_{t} = S + C + P </math>
In 2003, A.R. Tzafriri 2 proposed that the value , where and f(r), or the Van Slyke Cullen constant, equals .
Lastly, Banghe Li, Bo Li, and Yuefeng Shen 3 proposed two quasi-steady state laws. First, they stated, in much more rigorous mathematical language with appropriate proofs than stated here, that and . Second, they said the also has to be relatively equal to zero. However, when it comes to a particular epsilon value to determine the steady state, they say, “Epsilon can be any positive number that depends on the requirement of the experiments”.
Thus, former work has focused largely on under what reaction conditions will the quasi-steady state hold, but not particularly when such reactions are in the QSS. Thus, we attempt to do so here. Keeping in mind the work of former authors, our method of determining when the reaction is in the QSS is as follows:
We take the maximum value of the substrate-enzyme complex, which occurs within the steady state, and put the epsilon value at 1% of that value. Thus, anything within 1% of the maximum concentration is an indication of the solution remaining relatively constant and may be considered within the steady state.
Furthermore, we find the value of the derivative of the complex at a point very near the maximum concentration (for our example, we chose a time step of 0.1 seconds away). The reason for choosing this value is that the derivative at the maximum is zero. Then using that value, we say that if the derivative is within 1% of that value, thus indicating that the derivative is very near zero, then the reaction may be in the steady state.
In addition, if one takes the limit as the substrate approaches infinity, the concentration of the complex approaches the initial enzyme concentration . Thus, our third QSS criteria is that the concentration of the substrate enzyme complex is within 0.75 % of the initial enzyme concentration.
When all three conditions are met, the Matlab program returns the time values for when the reaction in question is in the quasi-steady state.
(In our example reaction, the quasi-steady-state was determined to last from 0.1 - 5.7 seconds)
After having established the criteria for determining when the reaction is in the QSS, we decided to explore adding forcing to the substrate differential equation. To start, we added the constant forcing alpha:
Cognizant of the fact that , I found an expression for alpha:
Note: This result is based on previous approximations that .
Further analysis revealed that as one increases the value of alpha above the calculated value provided by the preceding equation, the amount of time that the reaction is in the QSS grows exponentially.
Armed with this tool, my current work is focused on sinusoidal forcing in place of alpha, thus modeling a heartbeat or periodic eating.
The code for both programs, one with a user input of and and the other with a user input of and , follows:
% This code is written for the hydrolysis of benzoyl-L-arginine ethyl ester by
% trypsin. However, it can be modified to fit any single substrate-enzyme reaction.
% s = substrate
% e = enzyme
% se = substrate enzyme complex
% p = product
% conc = concentration
% k = reaction rates
global k_m k_cat; % Allowing variables to be used in other functions
% Reaction rates
k_m = 10e-6; % Michaelis-Menten Constant
k_cat = 15; % Rate-limiting turnover number (equal to k2)
% Initial concentrations (mol)
conce_init = 10e-5; % Usually nanomolar
concs_init = 10e-3; % Usually 5-6 orders of magnitude larger than conce_init
concse_init = 0; % Standard reaction analysis assumes the condition that
concp_init = 0; % the complex and product concentrations are initially zero
t_period = .1; % Time period between function evaluation points
t_end = 10; % Total time evaluated
% Calculating the solutions
[t,y] = ode15s ('system_template1',0:t_period:t_end,[concs_init,conce_init,concse_init,concp_init]);
concs = y(:,1);
conce = y(:,2);
concse = y(:,3); % Assigning values to a column of the solution array
concp = y(:,4);
% Setting up time for graphs (in seconds) - must be same as time evaluation
% in ode15s above
T = 0:t_period:t_end;
% Quasi-steady-state Approximation (Velocity of the Reaction)
V_qss =(k_cat*conce_init*concs)/(k_m + concs);
%Graphing
figure(1);
plot(T, conce);
ylabel('Enzyme Concentration (mol)');
xlabel('Time (s)');
figure(2);
plot(T, concs);
ylabel('Substrate Concentration (mol)');
xlabel('Time (s)');
figure(3);
plot(T, concse);
ylabel('Complex Concentration (mol)');
xlabel('Time (s)');
figure(4);
plot(T, concp);
ylabel('Product Concentration (mol)');
xlabel('Time (s)');
figure(6);
plot(T, V_qss);
ylabel('Reaction rate');
xlabel('Time');
Numerical = conce_init - concse;
Percentage = (Numerical/conce_init)*100;
% Test for when reaction is in quasi-steady-state
k1 = k_cat/k_m;
k_neg1 = k1*k_m - k_cat; % Equals 0
for i = 1:t_end/t_period + 1
concseprime(i) = k1*conce(i)*concs(i)-k_neg1*concse(i)-k_cat*concse(i); % Calculating dC/dt
end
max_concse = max(concse); % Obtaining the maximum value of the complex
interval = max_concse * 0.01; % Setting epsilon to 1% of max
for i = 1:t_end/t_period + 1
difference(i) = max_concse - concse(i); % Calculating the difference between the maximum value of the complex and
%each value of the function
if difference(i) < interval % Sets condition for QSS
[C, I] = max(concse);
value = concseprime(I-1) * 0.01; % Going one element to the left since the derivative at the max should be 0
differenceprime = concseprime(I) - concseprime(i);
if differenceprime < value
if Percentage(i) < .75
disp(i*t_period - t_period) % Returns exact time values for when reaction is in QSS
end
end
end
% This function is used in MichaelisMenten.m. It assigns differential
% equations for each chemical species to a column in a solution array,
% which is filled in MichaelisMenten by an ode solver.
function yprime = system_template1(t,y)
% Import global parameters
global k_m k_cat; % The Michaelis-Menten constant and rate-limiting turnover number
% Assigning the concentrations of each chemical species to a column in the array
concs = y(1); % Substrate
conce = y(2); % Enzyme
concse = y(3); % Substrate-Enzyme complex
concp = y(4); % Product
k1 = k_cat/k_m;
k_neg1 = k1*k_m - k_cat; % Equals 0
% Michaelis-Menten Equations
yprime(1,1) = -k1*conce*concs+k_neg1*concse;
yprime(2,1) = -k1*conce*concs+k_neg1*concse+k_cat*concse;
yprime(3,1) = k1*conce*concs-k_neg1*concse-k_cat*concse;
yprime(4,1) = k_cat*concse;
end
% Note: k_cat = k2 for a simple substrate-enzyme reaction as explained by
% Michaelis and Menten (1913)
% This code is written for the hydrolysis of benzoyl-L-arginine ethyl ester by
% trypsin. However, it can be modified to fit any single substrate-enzyme reaction.
% s = substrate
% e = enzyme
% se = substrate enzyme complex
% p = product
% conc = concentration
% k = reaction rates
global k_1 k_neg1 k_2; % Allowing variables to be used in other functions
% Reaction rates
k_1 = 4e6;
k_neg1 = 25;
k_2 = 15;
% Initial concentrations (mol)
conce_init = 10e-5; % Usually nanomolar
concs_init = 10e-3; % Usually 5-6 orders of magnitude larger than conce_init
concse_init = 0; % Standard reaction analysis assumes the condition that
concp_init = 0; % the complex and product concentrations are initially zero
t_period = .1; % Time period between function evaluation points
t_end = 10; % Total time evaluated
% Calculating the solutions
[t,y] = ode15s ('system_template2',0:t_period:t_end,[concs_init,conce_init,concse_init,concp_init]);
concs = y(:,1);
conce = y(:,2);
concse = y(:,3); % Assigning values to a column of the solution array
concp = y(:,4);
% Setting up time for graphs (in seconds) - must be same as time evaluation
% in ode15s above
T = 0:t_period:t_end;
k_m = (k_neg1 + k_2)/k_1; % Michaelis-Menten constant
% Quasi-steady-state Approximation (Velocity of the Reaction)
V_qss =(k_2*conce_init*concs)/(k_m + concs);
%Graphing
figure(1);
plot(T, conce);
ylabel('Enzyme Concentration (mol)');
xlabel('Time (s)');
figure(2);
plot(T, concs);
ylabel('Substrate Concentration (mol)');
xlabel('Time (s)');
figure(3);
plot(T, concse);
ylabel('Complex Concentration (mol)');
xlabel('Time (s)');
figure(4);
plot(T, concp);
ylabel('Product Concentration (mol)');
xlabel('Time (s)');
figure(6);
plot(T, V_qss);
ylabel('Reaction rate');
xlabel('Time');
Numerical = conce_init - concse;
Percentage = (Numerical/conce_init)*100;
% Test for when reaction is in quasi-steady-state
for i = 1:t_end/t_period + 1
concseprime(i) = k_1*conce(i)*concs(i)-k_neg1*concse(i)-k_2*concse(i); % Calculating dC/dt
end
max_concse = max(concse); % Obtaining the maximum value of the complex
interval = max_concse * 0.01; % Setting epsilon to 1% of max
for i = 1:t_end/t_period + 1
difference(i) = max_concse - concse(i); % Calculating the difference between the maximum value of the complex and %each value of the function
if difference(i) < interval % Sets condition for QSS
[C, I] = max(concse);
value = concseprime(I-1) * 0.01; % Going one element to the left since the derivative at the max should be 0
differenceprime = concseprime(I) - concseprime(i);
if differenceprime < value
if Percentage(i) < .75
disp(i*t_period - t_period) % Returns exact time values for when reaction is in QSS
end
end
end
% This function is used in MichaelisMenten1.m. It assigns differential
% equations for each chemical species to a column in a solution array,
% which is filled in MichaelisMenten1 by an ode solver.
function yprime = system_template2(t,y)
% Import global parameters
global k_1 k_neg1 k_2; % The reaction rates
% Assigning the concentrations of each chemical species to a column in the array
concs = y(1); % Substrate
conce = y(2); % Enzyme
concse = y(3); % Substrate-Enzyme complex
concp = y(4); % Product
% Michaelis-Menten Equations
yprime(1,1) = -k_1*conce*concs+k_neg1*concse;
yprime(2,1) = -k_1*conce*concs+k_neg1*concse+k_2*concse;
yprime(3,1) = k_1*conce*concs-k_neg1*concse-k_2*concse;
yprime(4,1) = k_2*concse;
end