- Email: [email protected]
Metabolic rate constants: Some computational aspects
Available online at www.sciencedirect.com
ScienceDirect Mathematics and Computers in Simulation (
)
– www.elsevier.com/locate/matcom
Original articles
Metabolic rate constants: Some computational aspects Stanko Dimitrov a,∗ , Svetoslav Markov b a University of Sofia, Faculty of Mathematics and Informatics, Bulgaria b Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Bulgaria
Received 31 December 2014; received in revised form 22 September 2015; accepted 7 November 2015
Abstract In this work we pose the question how reliable the Michaelis constant is as an enzyme kinetic parameter in situations when the Michaelis–Menten equation is not a good approximation of the true substrate dynamics as it may be in the case of metabolic processes in living cells. We compare the Michaelis–Menten substrate–product kinetics with the complete substrate–enzyme–product kinetics induced by the reaction scheme originally proposed by V. Henri. The Henri reaction scheme involves four concentrations and three rate constants and via the law of mass action is translated to a system of four ordinary differential equations (denoted as HMM-system). We propose a method for the computation of the three HMM rate constants, which can be applied in any situation whenever time course measurement data are available. The proposed method has been tested on the case study of acetylcholine hydrolysis. Our approach provides for the validation of the HMM-system by taking into account the uncertainties in the measurement data, and focuses on the use of contemporary computational tools. c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Enzyme kinetics; Reaction equations; Michaelis constant; Biomass–substrate–product dynamics; Computation of rate constants
1. Introduction Since its discovery around 1900 the Henri enzyme kinetic law has been intensively applied in many fields of life sciences. Two types (or stages) of its application can be distinguished: a first type of applications in biotechnological (fermentation) processes where the enzymes are often produced by bacterial cells, and a second type of applications in metabolic processes in living cells. The first type of applications are often denoted as “in vitro” whereas the second one—“in vivo”. During “in vitro” processes one is usually interested in the transformation of as much as possible quantity of substrate into product using thereby as little as possible quantity of enzymes. Thus the ratio enzyme/substrate “in vitro” is usually very small, whereas “in vivo” this may not be the case. For the “in vitro” case the approximate Michaelis–Menten model has been proposed together with a simple protocol for the calculation of the two parameters involved in the model (Vmax and K m ). The Michaelis–Menten model has proved to be extremely useful and has become very popular amongst biological scientists. Plausible biochemical mechanisms (known as Briggs–Haldane interpretations) have been found also known as (standard) quasi-steady-state approximation (sQSSA) criteria for the validity of the Michaelis–Menten model. With the necessity of considering “in vivo” processes, two approaches can be possibly followed. One is to extend the applicability (validity) of the Michaelis–Menten model by finding more general criteria than the standard enzyme/substrate ratio criterion (the sQSSA condition), such as the ∗ Corresponding author.
E-mail address: [email protected] (S. Dimitrov). http://dx.doi.org/10.1016/j.matcom.2015.11.003 c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.
2
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
so-called rQSSA and the tQSSA criteria. A disadvantage of this approach is that these criteria are more complex to check and involve the computation of new parameters, thereby still producing often very rough results. Another approach could be to make use of other approximate models together with corresponding validity criteria—such as the first order decay model. As demonstrated in Section 3 in certain situations the decay model may produce better approximation than the Michaelis–Menten model. A third approach would be to make use of the original (exact) Henri enzyme kinetic law and the induced system of ODEs. The proponents of the approximate models underline that the Henri law induced system of ODEs is too complex to use in real applications and cumbersome to handle. In our opinion this may have been the case some decades ago, however there exist powerful contemporary mathematical and computational tools to make the exact Henri model easily applicable. In particular, familiar computer algebra systems (CAS) such as Mathematica and MATLAB can be used to handle the ODE system induced by the Henri law. In order to apply efficiently the original Henri law it is necessary to compute the three Henri’s rate constants from time-course kinetic measurements using computer simulations. However, such time-course measurement data are usually not explicitly available in the literature, in the best case some graphs are published together with the computed Michaelis constant K m which is of little use for the computation of the three Henri rate constants. We remind that the Michaelis–Menten protocol does not make use of time-course measurement data; this probably is a possible reason for the popularity of the Michaelis–Menten model amongst experimental biochemists and biologists. Here we study mathematically and computationally experimental data paying special attention to the measurement errors involved. We describe and motivate our modeling approach applying V. Henri’s biochemical reaction scheme of the simple enzyme–substrate dynamics where two fractions of the enzyme (free and bound) are involved. Our approach is tested on a set of simulated time course measurements for the dynamics of acetylcholine hydrolysis by acetylcholinesterase [34]. Henri’s reaction scheme induces (via the mass action law) a system of ODEs further briefly denoted as HMMsystem. It is well-known that the Michaelis–Menten equation (MM-equation) for the substrate is an approximation of the HMM-system and this approximation is a good one only under certain conditions, for recent reviews see [6,28]. In this work we confirm that the MM-equation may produce rather rough approximations in certain situations. We focus on some contemporary computational tools that are available for dealing directly with the HMM-system, so that there is no need to make use of approximate models such as the MM-equation. Section 2 contains a review of the original Henri–Michaelis–Menten reaction scheme (HMM-scheme) and the approximate Michaelis–Menten model (MM-model). The HMM reaction scheme induces (via the mass action law) a system of ordinary differential equations (HMM-system) while the approximate MM-model is derived from the HMM-system using well known assumptions [4,3,5,18,16,21,22]. We briefly comment on the criteria for the “validity” of the MM-model. Section 3 is devoted to computational experiments that reveal some details on the conditions under which the approximate Michaelis–Menten model is adequate. We base our analysis on a well defined criterion for the validity of the model as well as general analysis of its dynamics based on the ODE system. Cases where the biomass solutions of the two models are close, but the suggested dynamics differ significantly are also presented. The use of the proposed validation criterion and analysis for the approximate model are encouraged before any analysis of the process dynamics based solely on the values of constants such as the Michaelis constant K m . Section 4 considers the computation of the HMM-system rate parameters using time course data. In particular, we make use of some Mathematica and MATLAB tools allowing for the estimation of the three rate parameters of the HMM-system by means of appropriate adjusting simulations of the system solutions to match available time course measurement data. 2. Enzyme kinetics basic models 2.1. Henri’s reaction scheme The following reaction scheme of a simple enzyme–substrate dynamics where two fractions of the enzyme (free and bound) are involved has been proposed by Victor Henri [7,13,14]: k1 k2 −→ S+E C −→ P + E. ←− k−1
(1)
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
3
Henri’s reaction scheme (1) describes the reaction mechanism between an enzyme E with a single active site and a substrate S, forming reversibly an enzyme–substrate complex C, which then yields irreversibly a product P. Reaction scheme (1) says that during the transition of the substrate S into product P the enzyme E bounds the substrate into a complex C having different properties than the free enzyme and thus being necessarily considered as a separate substance. Denoting the concentrations s = [S], e = [E], c = [C], p = [P] and applying the Mass Action Law to Henri’s reaction scheme (1) we obtain the following system of ODEs: ds/dt = −k1 es + k−1 c, de/dt = −k1 es + (k−1 + k2 )c, dc/dt = k1 es − (k−1 + k2 )c,
(2)
d p/dt = k2 c, to be further called the HMM-system—in tribute to V. Henri, as well to L. Michaelis and M. Menten [21]. If the three rate constants k1 , k−1 , k2 are known, system (2) can be treated as a Cauchy ODE problem with initial conditions s(0) = s0 > 0, e(0) = e0 > 0, c(0) = 0, p(0) = 0. In practice the HMM-system rate constants k1 , k−1 , k2 are often not known and have to be determined for every specific pair enzyme–substrate. The contemporary approach to this task is to consider the rate constants as parameters in the dynamic HMM-system (2) and to compute them by fitting the solutions of the system to time course experimentally measured data. This problem will be considered in Section 4. 2.2. Michaelis–Menten equation Applying the quasi-steady-state assumption (QSSA) to HMM-system (2) one derives the following approximate equations for the substrate rate ds/dt and the complex c [5,21,22]: e0 s . Km + s ds Vmax s =− = −µ(s) dt Km + s
c=
(3) (4)
where the right-hand side function µ(s) = Vmax s/(K m + s) is known as “specific growth function” in cell growth modeling [19]. The solution sm of Eq. (4) is an approximation of the solution s of (2). The QSSA is reasonably applied e.g. when the ratio [E]/[S] is small so that fermentation continues for a considerably long time interval during which the concentration c = [C] of the bound enzyme is nearly constant. The latter can be achieved e.g. when the condition e0 ≪ s0 holds. Let us also remind that under the validity of QSSA the complex concentration c = [E S] is expressed in terms of the MM-substrate concentration via (3). In their paper [21] Michaelis and Menten discuss in detail Henri’s reaction scheme and Eq. (4) known as the Michaelis–Menten equation to be further denoted as MM-model. In addition Michaelis and Menten proposed a protocol for the practical calculation of the constant Vmax and K m in (4). The constant K m is known as Michaelis constant [7,30]. The MM-equation (4) can be written in the form ds Vmax s Vmax =− =− , dt Km + s K m /s + 1 showing that for large values of s the right-hand side is close to the constant −Vmax ; hence the substrate uptake rate |ds/dt| is nearly equal to the constant Vmax . Eq. (4) produces a good approximation sm for s under certain conditions (sometimes called validity criteria) [4,5,10,12,18,28,25,26]. Thus the condition e0 ≪ s0 assures good approximation and is ubiquitous for many fermentation and biotechnological processes, but may not be present, e.g. in living cells [1,29]. In the sequel we discuss certain methods and tools for the computation of the Michaelis constant using both the MM-model (4) and the HMM-system (2). In particular, we are interested in computational methods for the evaluation of the three rate parameters of the HMM-system.
4
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
1.0
0.8
0.6
0.4
0.2
1
2
3
4
Fig. 1. Graphics of the substrate dynamics according to the HMM-system (2) and MM-model (4). The rate constants of (2) are k1 = 2.62, k−1 = 0.1, k2 = 1.25, initial conditions: s0 = 1, e0 = 1.5 (the system is dimensionless); the parameters of (4) are K m = (k−1 + k2 )/k1 = 0.51526, Vmax = k2 e0 = 1.875.
2.3. Relations between the rate parameters In Fig. 1 the substrate uptake s is visualized in two different ways. The two graphics present the approximate solution sm to the MM-model (4) as well as the original solution s of the HMM-system (2). In this example the kinetic parameters are chosen so that the validity conditions of the MM-model do not hold and the difference of the solutions for the substrate of the two systems is visible. We would also want to compare the two solutions in case the validity conditions of the approximate model do hold. In order to correctly compare the solutions s and sm one has to establish certain consistency relations between the parameters in the MM-equation and the HMM-system. From the derivation of the MM-equation, see e.g. [22], we know that the parameters in (4) can be expressed by the parameters in (2) by means of the relations: Vmax = k2 e0 ,
K m = (k−1 + k2 )/k1 .
(5)
On the other side the three parameters in (2) cannot be uniquely determined by the two parameters in (4). The following computational examples show how the approximate substrate concentration solution sm to the MM-equation may look like depending on the initial values of the substrate s0 and the enzyme e0 . According to the MM-model we have ds/dt = Vmax s/(K m + s). The two constants Vmax and K m can be experimentally determined according to the MM-protocol. If the MM-model is valid we have Vmax ≈ k2 e0 . From this expression one can calculate the value kcat = Vmax /e0 . Hence kcat ≈ k2 if the validity conditions for the MMmodel take place. Assuming that k2 is a universal rate constant for a given specific reaction, then it should not depend on whether the validity conditions take place or not (e.g. if k2 is computed in an experiment with e0 ≪ s0 , it should be the same constant in another experiment with e0 > s0 .) Thus we may assume kcat ≈ k2 . This gives a method for evaluating the HMM-rate constant k2 . However, it is not possible to obtain any information about the HMM-rate parameters k1 , k−1 on the basis of the MM-model, resp. the MM-protocol. 3. Computational examples: model comparison In this section we present the results of several numerical studies aiming the comparison of the substrate dynamics of the two models (2), (4). The values of the dimensionless parameters K = K m /s0 and ε = e0 /s0 play an important role in our analysis because they can significantly influence the behavior of the system. The values of the two system characteristics K , ε are varied in the numerical study in order to present results for the cases when K ≤ 1 and ε ≪ 1, ε > 1, as well as K ≫ 1 and ε ≪ 1, ε > 1. The values of the HMM rate constants k1 , k−1 , k2 are kept the
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
5
same for all cases as follows: k1 = 0.1, k−1 = 0.01, k2 = 10. The Michaelis constant consistent with these values is according to (5): K m = (k−1 + k2 )/k1 = 100.1. 3.1. Criterion for good approximation An important requirement for the validity of the MM-model (4) is the following condition: k2 e0 e0 1 ≪ 1. = s0 + K m 1 + (k−1 /k2 ) + (s0 k1 /k2 ) k1 (s0 + K m )2
(6)
Condition (6) can be derived from the relation between the estimates of two timescales tc and ts , see [18], or [22, ch. 6.2]. The timescale tc corresponds to the “faster” reaction, involving the enzyme and the complex. Complex formation is the most significant process during this timescale. The longer timescale, ts , corresponds to the product synthesis and substrate utilization. It is natural to assume tc ≪ ts which in turn leads to (6). We see that e0 ≪ s0 implies (6) but in case when e0 ≪ s0 does not hold, the MM-model could still be a good approximation if K m is large enough (K m ≫ s0 and also k−1 ≪ k2 ). In order to provide a more comprehensive explanation of this statement we should consider the non-dimensionalized system proposed in [22], as well as the form of the MM-model under the considered conditions: τ = k1 e0 t, λ = k2 /(k1 s0 ),
u(τ ) = s(t)/s0 , K = K m /s0 ,
v(τ ) = c(t)/e0 , ε = e0 /s0 ,
du/dτ = −u + (u + K − λ)v, εdv/dτ = u − (u + K )v,
u(0) = 1, v(0) = 0.
(7) (8)
First, the dynamics of v can be summarized as follows: v increases from t = 0 until it reaches its maximum at v = u/(u + K ), after that it decreases back to 0. An important point is that if K is “big” then the maximum of v will be smaller. Furthermore, larger K indicates either larger k2 and/or larger k−1 values in comparison to s0 k1 . The most common ranges of the values of k2 and k1 are considered to be k2 ∈ [103 , 106 ]s −1 , k1 ∈ [107 , 1010 ]M−1 s −1 from many sources (such as [2,31,33]). Some sources [31,24] also state that for “superefficient enzymes” (or enzymes with perfect kinetics) the rate-limiting step is the substrate–enzyme association step. Such enzymes are said to have reached kinetic perfection or catalytic perfection [32] and that k2 ≫ k−1 holds for them [31]. Additionally, in certain rare cases the concentration of substrate is very small in comparison to the available enzyme [1]. Such cases may include formation of cellular micro-compartments, enzymatic activity in mammalian muscle tissue and others [1]. We may assume that in these cases the concentration of substrate can even reach values s0 ≪ k2 /k1 . We will assume the following hold for some of the cases we study: (1) the substrate concentration is sufficiently low, so that k2 ≫ s0 k1 , (2) the formed enzyme–substrate complex breaks down to product and enzyme faster than it dissociates to substrate and enzyme, thus k2 ≫ k−1 . For the enzymes under consideration we can summarize that the rate-limiting step of their reactions is the substrate–enzyme association step and due to their catalytic perfection a small amount of complex is present during the whole process. It should be noted that for the case studies in the next section we have used values of the kinetic parameters different than the ones derived from real experiments for convenience and simplicity. However, the conditions we rely on (k2 ≫ s0 k1 , k2 ≫ k−1 ) are kept valid when examining the case K ≫ 1 and we argue that parameters close to the ranges observed in real experiments exist, based on the aforementioned statements. We could also rewrite the Michaelis–Menten model (4) under the assumptions K ≫ 1, k−1 ≪ k2 in the following form: Vmax s k2 e0 s k 1 k 2 e0 s ds =− = − k +k =− ∼ −k1 e0 s. −1 2 dt Km + s k−1 + k2 + k1 s +s k1
The obtained form suggests a much simpler reaction and a negligible role of the intermediate complex. Thus, the product formation is only limited by the kinetic constant k1 . The derived equation is also identical to the exact substrate ODE from (2) if we assume e = e0 , c ∼ 0, which would be natural given the kinetic parameter assumptions we have made.
6
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
3.2. The exponential decay model We may argue that under the conditions discussed above the second terms in the non-dimensionalized system (7)–(8) are almost 0. System (7)–(8) can be approximated as follows: du/dτ = −u,
dv/dτ =
u ; ε
u(0) = 1, v(0) = 0.
(9)
The solution u of this dynamical system is identical to the solution sd of the substrate dynamics induced by the simple reaction scheme: k
S + E −→ P + E,
(10)
wherein S is the substrate, E is the enzyme and P is the product. Reaction scheme (10) says, differently to (1), that during the transition of the substrate S into product P the enzyme E does not bound the substrate into a complex. Applying the mass action law, assuming e = [E] constant, the kinetic scheme (10) leads to the following “exponential decay” differential equation for the substrate concentration sd = [S]: dsd = −kesd . dt
(11)
Solving the Cauchy problem related to Eq. (11), having set e = e0 , we obtain k = k1 . Indeed, the solution of Eq. (11) is sd (t) = s0 e−e0 kt . If we then solve the Cauchy problem (9), we would obtain s(t) = s0 e−e0 k1 t after converting it back to the initial variables. Our analysis shows that model (11) deserves to be included in the comparison of the two models (2), (4). We shall further show that under the conditions K ≫ 1, k−1 ≪ k2 the three models (2), (4) and (11) do indeed provide almost identical results and that the substrate solution of the MM-model may provide a good approximation even though e0 ≪ s0 does not necessarily hold. 3.3. Computational examples: three cases Case 1. K ≈ 1 (s0 ≈ K m ) (also valid for the range K ≤ 1 (s0 ≥ K m )). Computational example 1.1. In this example the MM–HMM-solutions for the substrate are close when ε is “small” as seen in Fig. 2. The values of the initial conditions used in the models are as follows: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0. The values of the parameters in (4) are consistent with the initial conditions and rate parameters in (2) and are calculated using (5) as Vmax = k2 e0 , K m = (k−1 + k2 )/k1 . The parameter k in the “exponential decay” differential equation (11) is set to be equal to k1 as derived in the previous section. In all the examples that follow, only s0 and e0 will vary, all the other model parameters will have the same values as in this example or they will be calculated using the corresponding s0 , e0 . This numerical example demonstrates that when the condition ε ≪ 1 holds then the MM-model (4) can be a very good approximation of the original HMM-system (2). Note that in this example the Michaelis constant K m used for the computation of the approximate MM-solution is derived from the HMM-rate parameters k1 , k−1 , k2 . This means that both models describe the process dynamics using equivalent rate constants and, since we consider the HMM-system to be true, this implies that the approximate model is also valid under the assumption ε ≪ 1. Our next numerical examples aim to examine what happens whenever the assumption ε ≪ 1 does not hold. Computational example 1.2. Our second numerical example shows how the solutions start to deviate when s0 and e0 are close to each other (Fig. 3). We can observe that with ε ∼ 1 the transient phase, related to the timescale tc , takes much longer and the assumption that tc ≪ ts is violated, thus invalidating the use of the MM-model. The goodness of the approximation of the complex c can be a better indicator of the applicability of the MM-model than the approximation of the substrate s, as c undergoes a boundary layer effect in the beginning of the reaction. The values of the initial conditions used are as follows: s0 = 100, e0 = 50, c0 = 0, p0 = 0. As before, we use (5) to calculate Vmax = k2 e0 , K m = (k−1 + k2 )/k1 .
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
7
Fig. 2. The solutions of the substrate (for models (2), (4), (11)) and the complex (for models (2), (4)) in the case ε = 10−3 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 +k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0. Under these conditions the MM-model clearly serves as a good approximation of the exact HMM-system.
Computational example 1.3. In this numerical example ε > 1. The HMM-solution for the substrate is even closer to the exponential decay solution (Fig. 4). As observed from the approximations of the complex in the same figure, with ε > 1, the MM-model is inapplicable due to the much longer timescale of the complex buildup. The values of the initial conditions used for the given solution are as follows: s0 = 100, e0 = 400, c0 = 0, p0 = 0. Here, and everywhere in the sequel Vmax = k2 e0 , K m = (k−1 + k2 )/k1 using (5). A figure of the specific growth rate function µ(s) = Vmax s/(K m + s) is included for reference, see Fig. 5. Case 2. K > 1. Computational example 2.1. It is clear that for values of K larger than 1 (K ≈ 2 in this example) but not K ≫ 1, the relations between the substrate solutions of the three models (2), (4) and (11) are preserved (Fig. 6). The values of the initial conditions used in the models are as follows: s0 = 50, e0 = 0.1, c0 = 0, p0 = 0. Computational example 2.2. In the second numerical example of this case the substrate solutions deviate when s0 and e0 are close to each other (Fig. 7), just as in example 1.2 and we can draw the same conclusion from the behavior of the complex, namely the assumption that tc ≪ ts has been violated, invalidating the use of the MM-model. All these follow from our choice of s0 and e0 (ε = 1). The values of the initial conditions used for the given solution are as follows: s0 = 50, e0 = 50, c0 = 0, p0 = 0. Computational example 2.3. In this numerical example ε > 1 (Fig. 8). The solutions are identical to the solutions of example 1.3. The approximation of the complex c in Fig. 8 can lead us to the conclusion that the MM-model is
8
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Fig. 3. The solutions of the substrate (for models (2), (4), (11)) and the complex (for models (2), (4)) in the case ε = 0.5 ∼ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 500; initial conditions: s0 = 100, e0 = 50, c0 = 0, p0 = 0.
inapplicable due to the much longer timescale of the complex buildup. This again follows from our choice of initial conditions s0 = 50, e0 = 150, meaning ε = 3 > 1. The values of the initial conditions used for the given solution are as follows: s0 = 50, e0 = 150, c0 = 0, p0 = 0. A figure of the specific growth rate function µ(s) = Vmax s/(K m + s) is included for reference, see Fig. 9. Case 3. K ≫ 1. Computational example 3.1. In this numerical example ε ≪ 1 and the solutions of the three models (2), (4) and (11) are identical (Fig. 10). The values of the initial conditions used for the given solution are as follows: s0 = 0.1, e0 = 0.001, c0 = 0, p0 = 0. Computational example 3.2. Although we have set the initial conditions so that ε ≫ 1, the substrate solutions of the three models (2), (4) and (11) are still very close to each other (Fig. 11). The values of the initial conditions used in the models are as follows: s0 = 0.1, e0 = 2000, c0 = 0, p0 = 0. A figure of the specific growth rate function µ(s) = Vmax s/(K m + s) is included for reference, cf. Fig. 12. As described above, in case 3: K ≫ 1 (also keeping in mind the assumption k−1 ≪ k2 ), the MM-model can serve as a very good approximation to the Henri–Michaelis–Menten law in terms of substrate dynamics, regardless of the range of initial conditions. As seen from the solutions of the complex though (Fig. 11), the MM-model completely fails to provide a good approximation of the complex dynamics when ε ≫ 1. This is expected since the approximate model only provides the outer solutions of the system in terms of singular perturbation analysis and in this case the
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
9
Fig. 4. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 4 > 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 4000; initial conditions: s0 = 100, e0 = 400, c0 = 0, p0 = 0.
Fig. 5. The specific growth rate function µ(s) = Vmax s/(K m + s) from (4) calculated in the range s from 0 to s0 = 100. Case—ε = 10−3 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0.
10
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Fig. 6. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 0.002 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 50, e0 = 0.1, c0 = 0, p0 = 0.
inner solution has to be accounted for due to the longer timescale of the complex buildup process. According to the MM-model, the solution of the enzyme–substrate complex can be derived as c(t) = e0 s(t)/(K m + s(t)), which under the proposed conditions (K ≫ 1, k−1 ≪ k2 and ε ≫ 1) cannot serve as a good approximation of the exact complex behavior. If we consider the non-dimensionalized system (9) under the same conditions, we can conclude that dv/dτ ∼ 0 hence v = const. The specific growth rate function visualized on Fig. 12 is almost linear in this case (again, regardless of whether ε ≪ 1 or ε ≫ 1) in contrast to cases 1 and 2. It is thus natural to expect the approximate model to predict substrate solutions similar to exponential decay in this case. Finally, the criterion (6) is obviously violated for computational example 3.2. (in contrast to all other examples where it is evaluated to less than 1): k2 e0 = 19.9. k1 (s0 + K m )2 Nevertheless, since K ≫ 1, k−1 ≪ k2 , the proposed explanations in the beginning of the section are indeed manifested and the substrate solutions of all three considered models are identical. As we also noted, this does not necessarily mean that the solutions of the other molecular species are correct, e.g. the MM-model solutions for the complex are incorrect in case ε ≫ 1.
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
11
Complex
Fig. 7. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 500; initial conditions: s0 = 50, e0 = 50, c0 = 0, p0 = 0.
Conclusion. There exist theoretical cases when the decay model produces better approximation to the exact HMMscheme than the MM-model. Additionally, in certain cases in which the MM-model can provide good substrate approximations, it may completely fail to give correct estimates of the complex or enzyme behavior. 3.4. Is the Michaelis constant a reliable kinetic parameter? Examples 1.2., 1.3., 2.2. and 2.3. demonstrate that the approximate model’s solutions are far from those of the exact HMM-system despite of the consistency of the parameters used in (2) and (4) even though we derive the Vmax and K m constants from the solutions of the exact model. This, of course, follows from the fact that under the initial conditions we have set the approximations embedded in the MM-model do not hold. However, if we free the MMmodel parameters from the conditions binding them to the exact scheme (5), we may still be able to find a solution closer to the exact one (compare Figs. 13(a) and 14 to Figures resp. 4 and 8). On the other hand, even if such a solution exists, it would correspond to completely different dynamics as we could easily determine from the solutions of the complex c in Fig. 13(b). The duration of the transient phase, as observed in Fig. 13(b), is related to the parameter ε (ε > 1 in the considered case) and it indicates that the MM-model may not be suitable under the given conditions since the approximations embedded in it are in contradiction to the behavior of the system. As a consequence any parameter estimation would suggest values of the kinetic parameters much different than the real ones. Although other sources [26] have also revealed the flaws of the inappropriate application of such approximate models under conditions which do not hold, some authors still prefer to fit their experimental data with them. As in Fig. 13(a) the solutions may appear to be close to those of the exact scheme but the dynamics that the approximate
12
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
complex
Fig. 8. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 3 > 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1500; initial conditions: s0 = 50, e0 = 150, c0 = 0, p0 = 0.
Fig. 9. The specific growth rate function µ(s) = Vmax s/(K m + s) from (4) calculated in the range s from 0 to s0 = 50. Case — ε = 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 500; initial conditions: s0 = 50, e0 = 50, c0 = 0, p0 = 0.
model suggests could be completely different than the real one (i.e. the MM-model parameters used in Figs. 13(a), 14 are respectively (K m = 200, Vmax = 9000) for Case 1 and (K m = 85, Vmax = 1500) for Case 2; in both cases K m can be evaluated to K m = 100.1 while Vmax = 4000 in Case 1 and Vmax = 1500 in Case 2).
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
13
Complex 1.X10-6
3.X10-7
5.X10-7
4.X10-7
2.X10-7
10000
20000
30000
40000
50000
Time
Fig. 10. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 0.01 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 0.01; initial conditions: s0 = 0.1, e0 = 0.001, c0 = 0, p0 = 0.
As explained in the previous section it is always beneficial to verify the criterion for the validity of the approximate model before applying it. Some general analysis of the behavior of the ODE system, based on the initial conditions, could also be valuable. We can similarly find relatively good fits of the decay model (11) for Examples 1.1. and 2.1. if we unbind the rate constant k. Figs. 15 and 16 show that in some cases the exponential decay solution may still approximate the system behavior to a certain degree, although the rate constant k of the decay model may not be strictly connected to the parameters of the exact HMM-model. Of course, if ε is small enough so that the substrate uptake is almost linear, no solution of (10) would be close to the exact solutions produced by (2). Conclusions. The presented computational examples show that the concept of “reliability” of the approximate models (the MM-model and the decay model) is not well defined. This concept obtains different meanings depending on whether one keeps the parameters of the approximate models bound to the HMM-system parameters, or not. 4. Computation of the HMM-system rate parameters using time course data 4.1. General settings The dynamics of the following biochemical process is studied using the HMM-system. The evaluation of the models is based on the goodness of fit achieved with them, keeping in mind the measurement errors contained in the data. The correctness of the dynamics suggested by the model is assessed by deriving the Michaelis constant for
14
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Time Complex
Time
Fig. 11. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 2 ∗ 104 ≫ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 2 ∗ 104 ; initial conditions: s0 = 0.1, e0 = 2000, c0 = 0, p0 = 0.
Fig. 12. The specific growth rate function µ(s) = Vmax s/(K m +s) from (4) calculated in the range s from 0 to s0 = 0.1. Case — ε = 2∗104 ≫ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 2 ∗ 104 ; initial conditions: s0 = 0.1, e0 = 2000, c0 = 0, p0 = 0.
each experiment and comparing it to known values from the literature corresponding to the same biochemical process. The conclusion is that approximate models should be avoided as much as possible. In the next section we propose
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
15
Substrate concentration
Time
(a) A substrate solution of the MM-model closer to the exact one may exist if we free the MM-model parameters from the conditions binding them to the exact scheme (5). Complex
120 100 80 60 40 20
0.02
0.04
0.06
0.08
0.10
0.12
0.14
Time
(b) The poor approximation of the complex indicates that the MM-model may not be applicable in the presented case. Fig. 13. The solutions are as in Example 1.3 (ε = 4 > 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 4000; initial conditions: s0 = 100, e0 = 400, c0 = 0, p0 = 0), with the only difference in (4) being: K m = 200, Vmax = 9000.
contemporary tools allowing the computation of the HMM-system parameters directly from time course measurement data. The process we follow is similar to the parameter fitting procedures and reverse engineering described in other papers, such as [26,27,23]. 4.2. Case study: fitting the HMM-system against experimental data of acetylcholine hydrolysis We fit the HMM model to experimental data of acetylcholine hydrolysis by acetylcholinesterase that has been obtained using the derived rate constants from the HXA method of the experiments in [34]. The rate constants were used to compute the solutions of (2) and evaluate them for a number of time points ti . Finally, a certain amount of noise was added to the data according to the standard deviations given in [34]. We take into account the measurement errors in the data and they are displayed as concentration intervals for each point in the plots. Our computational problem can be formulated as follows. Given time course experimental data (together with measurement errors) for the substrate, enzyme and product concentrations find values for the parameters k−1 , k1 , k2 and initial values E 0 , s0 such that the solution of the HMM system (2) fit well against the experimental data and
16
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Time
Fig. 14. The solutions are as in Case 2, Example 3 (ε = 3 > 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1500; initial conditions: s0 = 50, e0 = 150, c0 = 0, p0 = 0), with the only difference in (4) being: K m = 85. Substrate concentration
Time
Fig. 15. The solutions are as in Case 1, Example 1 (ε = 10−3 ≪ 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m =
(k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0), with the only difference in (11) being: k = 0.07, rather than k = 0.1.
possibly fit into the measurement intervals. Our procedure for solving this problem passes in two stages: first we find an initial rough “guess” for the parameter values, then we iteratively improve the parameter set, resp. the solutions, using some optimization routines of the numerical computing environment MATLAB. More precisely, we have used MATLAB’s lsqnonlin procedure (optimization algorithm defaults to Trust-Region-Reflective Algorithm), where we aim to minimize the difference between the experimental data and the computed enzyme kinetics model solutions (using ode23 or ode23s solvers) for a given set of rate constants and initial conditions. At the first step, the solution of the ODEs HMM-system (2) with initial “guess” parameters (Fig. 17) is shown in Fig. 18. The values of the parameters are as follows: k1 = 25000M−1 s −1 , k−1 = 0s −1 , k2 = 10s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M. The results of the iterative optimization procedure we applied using the initial parameters can be seen in Fig. 19. The values of the parameters used for the given solution are as follows: k1 = 16847M−1 s −1 , k−1 = 7s −1 , k2 = 12s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M. Using the parameter set we have found to solve the HMM-system we get solutions that fit nicely the experimental data we started with. Using K m = (k−1 +k2 )/k1 we can derive the Michaelis constant from the optimal parameters in the HMM-system, K m = 0.00112M. This value is close to known values of the Michaelis constant from other sources [9,11], hence we can conclude that the model estimates correctly the behavior of the underlying biochemical process.
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
17
Substrate concentration 50
40
30
20
10
0
0
100
200
300 Time
400
500
600
Fig. 16. The solutions are as in Case 2, Example 1 (ε = 0.002 ≪ 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 50, e0 = 0.1, c0 = 0, p0 = 0), with the only difference in (11) being: k = 0.08, rather than k = 0.1.
Fig. 17. An example of the usage of dynamical numerical solution of the ODE system and their graphical representation using the Mani pulate and N DSolve functions of CAS Mathematica. An initial rough fit of the experimental data is obtained in this example.
The MM-model (4) on the other hand provides very different results for the same experimental data (Fig. 20). Although the obtained parameters (Vmax = 7.12 ∗ 10−7 Ms −1 ; K m = 2.193 ∗ 10−3 M) are relatively close to the real ones, the solutions do not pass through the intervals of the data. As in Example 2 and Example 3 of the first two Cases, the approximate model’s solutions deviate from the real ones even though its parameters may be close to the exact model’s rate constants. This can easily be explained by the fact that the complex c has more complicated dynamics (see Fig. 19 for the dynamics of c and e) than what the MM-model can estimate, i.e. data points from the transient phase cannot be accounted for. The obtained parameters are thus unreliable as they do not correspond to a good approximation of the complex c and the enzyme e. A rather good fit using the same approximate model (4) can be achieved by omitting the first observation data point of the enzyme (Fig. 21). The values of the model parameters (Vmax = 6.6 ∗ 10−7 Ms −1 ; K m = 1.26 ∗ 10−3 M) are even closer to those found using the HMM model. In general, we cannot expect data points related to the “faster” reaction (taking place during the tc timescale) to be approximated well with the MM-model because it is entirely focused on the “slower” reaction and the solutions we get cannot possibly follow the dynamics during both timescales.
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation ( Enzyme Kinetics Model
Concentration (substrate,product [mM])
x10-3
)
–
x10-7 substrate - s product - p experimental data - s experimental data - p enzyme e complex c experimental data - e errors
3
2 0.6
0.4 1
concentration - enzyme, complex [mM]
18
0.2
0
0
2000
4000
6000
8000
10000
12000
0 14000
Fig. 18. The numerical solution of the HMM-system with an initial parameter “guess” set of the unknown parameters k1 = 25000M−1 s −1 , k−1 = 0s −1 , k2 = 10s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M. Enzyme Kinetics Model
x10-7 substrate - s product - p experimental data - s experimental data - p enzyme e complex c experimental data - e errors
2.5
2 0.6 1.5
0.5 0.4
1 0.3
concentration - enzyme, complex [mM]
Concentration (substrate,product [mM])
x10-3
0.2
0.5
0.1 0
0
2000
4000
6000
8000
10000
12000
14000
Fig. 19. The numerical solution of the HMM-system (2) after the fitting of the model. The parameters and initial values used for the given solution are k1 = 16847M−1 s −1 , k−1 = 7s −1 , k2 = 12s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M.
5. Concluding remarks In this work we demonstrate that the value of the Michaelis constant may strongly depend on the way one makes use of the experimental data and the method of calculation. On the other side the rate constants in the enzyme kinetic Henri’s reaction scheme are well defined and can be efficiently computed from time course experimental data. Moreover, the modeling of metabolic pathways makes use of multiple reactions and can be made more rigorous when based on Henri’s reactions instead of Michael–Menten approximate reactions, being especially imprecise when applied to living cells. Such arguments confirm the need of refining the general methods for the determination of the
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
Enzyme Kinetics Model
x10-3
19
– x10-7
experimental data - s experimental data - p enzyme e complex c errors experimental data - e 2.5
2
0.7 0.6
1.5
0.5 0.4
1 0.3 0.2
0.5
0
0
2000
4000
6000
8000
10000
12000
concentration - enzyme, complex [mM]
Concentration (substrate,product [mM])
substrate - s product - p
14000
Fig. 20. The numerical solution of the MM-model after fitting it to the experimental data. Model parameters obtained from the optimization procedure: Vmax = 7.12 ∗ 10−7 Ms −1 ; K m = 2.193 ∗ 10−3 M. Enzyme Kinetics Model
x10-3
x10-7
experimental data - p enzyme e complex c errors experimental data - e 2.5
2 0.6 1.5
0.5 0.4
1 0.3
concentration - enzyme, complex [mM]
Concentration (substrate,product [mM])
substrate - s product - p experimental data - s
0.2
0.5
0.1 0 0
2000
4000
6000
8000
10000
12000
14000
Fig. 21. The numerical solution of the MM-model after fitting it to the experimental data. The first experimental data point for the enzyme has been omitted in the optimization. Model parameters obtained from the optimization procedure are now Vmax = 6.6 ∗ 10−7 Ms −1 ; K m = 1.26 ∗ 10−3 M.
Henri’s reaction scheme rate constants. Such general methods using contemporary software tools are suggested and discussed in the present work. In addition we wish to point out the need of validation methods in combination with the availability of experimental data involving measurement errors [17]; a more detailed introduction of the use of ⁀ validated interval methods is presented in [8]. In this work we make use of CAS Mathematica as computational tool, however other suitable software tools can be used, such as MATLAB and COPASI [15,20]. Ensuring the validity of the Michaelis–Menten model has been identified as a significant step in its application which has led to the comprehensive studies of validation criteria [28].
20
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
When modeling kinetic reactions we suggest to start with a description of the process by a biochemical reaction equations (schemes) and then pass to mathematical differential equations via the mass action principle. Such an approach is useful in keeping close to the biochemical mechanism of the modeling process; we recommend the use of this approach not only in fields like enzyme kinetic and metabolic networks but also in fields like cell growth and population dynamics [19]. Acknowledgments The authors are very grateful to the anonymous reviewers whose valuable reports contributed much to the quality of the paper. References [1] K.R. Albe, M.H. Butler, B.E. Wright, Cellular concentrations of enzymes and their substrates, J. Theoret. Biol. 143 (1990) 163–195. [2] M.N. Berberan-Santos, A general treatment of Henri–Michaelis–Menten enzyme kinetics: exact series solution and approximate analytical solutions, MATCH Commun. Math. Comput. Chem. 63 (2010) 283–318. [3] A.M. Bersani, E. Bersani, G. Dell’Acqua, M.G. Pedersen, New trends and perspectives in nonlinear intracellular dynamics: one century from Michaelis–Menten paper, Contin. Mech. Thermodyn. 27 (2014) 659–684. [4] A.M. Bersani, G. Dell’Acqua, Is there anything left to say on enzyme kinetic constants and QSSA? J. Math. Chem. 50 (2012) 335–344. [5] G.E. Briggs, J.B.S. Haldane, A note on the kinetic of enzyme action, Biochem. J. 19 (1925) 338–339. [6] W.W. Chen, M. Niepel, P.K. Sorger, Classic and contemporary approaches to modeling biochemical reactions, Genes. Dev. 24 (2010) 1861–1875. [7] U. Deichmann, S. Schuster, J.-P. Mazat, A. Cornish-Bowden, Commemorating the 1913 Michaelis–Menten paper Die Kinetik der Invertinwirkung: three perspectives, FEBS J. 281 (2014) 435–463. J.-P. Mazat, Part 3: pp.452–463. [8] S. Dimitrov, G. Velikova, V. Beschkov, S. Markov, On the numerical computation of enzyme kinetic parameters, Biomath Commun. 1 (2) (2014). [9] A. Forsberg, Kinetics for the inhibition of acetylcholinesterase from the electric eel by some organophosphates and carbamates, Eur. J. Biochem. 140 (1984) 153–156. [10] R. Grima, N.G. Walter, S. Schnell, Single molecule enzymology a` la Michaelis–Menten, FEBS J. 281 (2014) 518–530. [11] A. Hai, et al., Acetylcholinesterase-ISFET based system for the detection of acetylcholine and acetylcholinesterase inhibitors, Biosens. Bioelectron. 22 (2006) 605–612. [12] S.M. Hanson, S. Schnell, The reactant stationary approximation in enzyme kinetics, J. Phys. Chem. A 112 (2008) 8654–8658. [13] V. Henri, Recherches sur la loi de l’action de la sucrase, C. R. Hebd. Acad. Sci. 133 (1901) 891–899. [14] V. Henri, Lois g´en´erales de l’action des diastases, Librairie Scientifique A. Hermann, Paris, 1903. [15] S. Hoops, et al., COPASI—a complex pathway simulator, Bioinformatics 22 (24) (2006) 3067–3074. [16] K.A. Johnson, R.S. Goody, The original Michaelis constant: translation of the 1913 Michaelis–Menten paper, Biochemistry 50 (39) (2011) 8264–8269. [17] W. Kraemer, J. Wolff, v. Gudenberg (Eds.), Scientific Computing, Validated Numerics, Interval Methods, in: Proc. of SCAN-2000/Interval2000, Kluwer/Plenum, 2001. [18] K.J. Laidler, Theory of the transient phase in kinetics, with special reference to enzyme systems, Can. J. Chem. 33 (1955) 1614–1624. [19] S. Markov, Cell growth models using reaction schemes: batch cultivation, Biomath. 2 (2) (2013) 1312301–9. [20] P. Mendes, S. Hoops, S. Sahle, R. Gauges, J. Dada, U. Kummer, Computational modeling of biochemical networks using COPASI, Methods Mol. Biol. 500 (2009) 17–59. [21] L. Michaelis, M.L. Menten, Die Kinetik der Invertinwirkung, Biochem. Z. 49 (1913) 333–369. [22] J.D. Murray, Mathematical Biology: I. An Introduction, third ed., Springer, 2002. [23] N. Nikolova, K. Tenekedjiev, K. Kolev, Uses and misuses of progress curve analysis in enzyme kinetics, Cent. Eur. J. Biol. 3 (4) (2008) 345–350. [24] M.H.M. Olsson, P.E.M. Siegbahn, A. Warshel, Simulations of the large kinetic isotope effect and the temperature dependence of the hydrogen atom transfer in lipoxygenase, J. Am. Chem. Soc. 126 (9) (2004) 2820–2828. [25] M.G. Pedersen, A.M. Bersani, E. Bersani, The quasi steady-state approximation for fully competitive enzyme reactions, Bull. Math. Biol. 69 (2007) 433–457. [26] M.G. Pedersen, A.M. Bersani, E. Bersani, G. Cortese, The total quasi-steady-state approximation for complex enzyme reactions — a word of caution, Math. Comput. Simul. 79 (2008) 1010–1019. [27] T.D. Pollard, E.M. De La Cruz, Take advantage of time in your experiments: a guide to simple, informative kinetics assays, Mol. Biol. Cell. 24 (8) (2013) 1103–1110. [28] S. Schnell, Validity of the Michaelis–Menten equation — steady-state, or reactant stationary assumption: that is the question, FEBS J. 281 (2014) 464–472. [29] S. Schnell, P.K. Maini, Enzyme kinetics at high enzyme concentration, Bull. Math. Biol. 62 (2000) 483–499. [30] S. Schnell, P.K. Maini, A century of enzyme kinetics: Reliability of the K m and vmax estimates, Comments Theor. Biol. 8 (2003) 169–187. [31] M.E. Stroppolo, M. Falconi, A.M. Caccuri, A. Desideri, Superefficient enzymes, Cell. Mol. Life Sci. 58 (2001) 1451–1460. [32] L. Stryer, J.M. Berg, J.L. Tymoczko, Biochemistry, fifth ed., W H Freeman, New York, 2002. [33] I.B. Wilson, M.A. Harrison, Turnover number of acetylcholinesterase, J. Biol. Chem. 236 (8) (1961) 2292–2295. [34] P. Zdrazilova, et al., Kinetics of total enzymatic hydrolysis of acetylcholine and acetylthiocholine, Z. Nat.forsch. C. 61 (3–4) (2006) 289–294.
ScienceDirect Mathematics and Computers in Simulation (
)
– www.elsevier.com/locate/matcom
Original articles
Metabolic rate constants: Some computational aspects Stanko Dimitrov a,∗ , Svetoslav Markov b a University of Sofia, Faculty of Mathematics and Informatics, Bulgaria b Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Bulgaria
Received 31 December 2014; received in revised form 22 September 2015; accepted 7 November 2015
Abstract In this work we pose the question how reliable the Michaelis constant is as an enzyme kinetic parameter in situations when the Michaelis–Menten equation is not a good approximation of the true substrate dynamics as it may be in the case of metabolic processes in living cells. We compare the Michaelis–Menten substrate–product kinetics with the complete substrate–enzyme–product kinetics induced by the reaction scheme originally proposed by V. Henri. The Henri reaction scheme involves four concentrations and three rate constants and via the law of mass action is translated to a system of four ordinary differential equations (denoted as HMM-system). We propose a method for the computation of the three HMM rate constants, which can be applied in any situation whenever time course measurement data are available. The proposed method has been tested on the case study of acetylcholine hydrolysis. Our approach provides for the validation of the HMM-system by taking into account the uncertainties in the measurement data, and focuses on the use of contemporary computational tools. c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Enzyme kinetics; Reaction equations; Michaelis constant; Biomass–substrate–product dynamics; Computation of rate constants
1. Introduction Since its discovery around 1900 the Henri enzyme kinetic law has been intensively applied in many fields of life sciences. Two types (or stages) of its application can be distinguished: a first type of applications in biotechnological (fermentation) processes where the enzymes are often produced by bacterial cells, and a second type of applications in metabolic processes in living cells. The first type of applications are often denoted as “in vitro” whereas the second one—“in vivo”. During “in vitro” processes one is usually interested in the transformation of as much as possible quantity of substrate into product using thereby as little as possible quantity of enzymes. Thus the ratio enzyme/substrate “in vitro” is usually very small, whereas “in vivo” this may not be the case. For the “in vitro” case the approximate Michaelis–Menten model has been proposed together with a simple protocol for the calculation of the two parameters involved in the model (Vmax and K m ). The Michaelis–Menten model has proved to be extremely useful and has become very popular amongst biological scientists. Plausible biochemical mechanisms (known as Briggs–Haldane interpretations) have been found also known as (standard) quasi-steady-state approximation (sQSSA) criteria for the validity of the Michaelis–Menten model. With the necessity of considering “in vivo” processes, two approaches can be possibly followed. One is to extend the applicability (validity) of the Michaelis–Menten model by finding more general criteria than the standard enzyme/substrate ratio criterion (the sQSSA condition), such as the ∗ Corresponding author.
E-mail address: [email protected] (S. Dimitrov). http://dx.doi.org/10.1016/j.matcom.2015.11.003 c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.
2
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
so-called rQSSA and the tQSSA criteria. A disadvantage of this approach is that these criteria are more complex to check and involve the computation of new parameters, thereby still producing often very rough results. Another approach could be to make use of other approximate models together with corresponding validity criteria—such as the first order decay model. As demonstrated in Section 3 in certain situations the decay model may produce better approximation than the Michaelis–Menten model. A third approach would be to make use of the original (exact) Henri enzyme kinetic law and the induced system of ODEs. The proponents of the approximate models underline that the Henri law induced system of ODEs is too complex to use in real applications and cumbersome to handle. In our opinion this may have been the case some decades ago, however there exist powerful contemporary mathematical and computational tools to make the exact Henri model easily applicable. In particular, familiar computer algebra systems (CAS) such as Mathematica and MATLAB can be used to handle the ODE system induced by the Henri law. In order to apply efficiently the original Henri law it is necessary to compute the three Henri’s rate constants from time-course kinetic measurements using computer simulations. However, such time-course measurement data are usually not explicitly available in the literature, in the best case some graphs are published together with the computed Michaelis constant K m which is of little use for the computation of the three Henri rate constants. We remind that the Michaelis–Menten protocol does not make use of time-course measurement data; this probably is a possible reason for the popularity of the Michaelis–Menten model amongst experimental biochemists and biologists. Here we study mathematically and computationally experimental data paying special attention to the measurement errors involved. We describe and motivate our modeling approach applying V. Henri’s biochemical reaction scheme of the simple enzyme–substrate dynamics where two fractions of the enzyme (free and bound) are involved. Our approach is tested on a set of simulated time course measurements for the dynamics of acetylcholine hydrolysis by acetylcholinesterase [34]. Henri’s reaction scheme induces (via the mass action law) a system of ODEs further briefly denoted as HMMsystem. It is well-known that the Michaelis–Menten equation (MM-equation) for the substrate is an approximation of the HMM-system and this approximation is a good one only under certain conditions, for recent reviews see [6,28]. In this work we confirm that the MM-equation may produce rather rough approximations in certain situations. We focus on some contemporary computational tools that are available for dealing directly with the HMM-system, so that there is no need to make use of approximate models such as the MM-equation. Section 2 contains a review of the original Henri–Michaelis–Menten reaction scheme (HMM-scheme) and the approximate Michaelis–Menten model (MM-model). The HMM reaction scheme induces (via the mass action law) a system of ordinary differential equations (HMM-system) while the approximate MM-model is derived from the HMM-system using well known assumptions [4,3,5,18,16,21,22]. We briefly comment on the criteria for the “validity” of the MM-model. Section 3 is devoted to computational experiments that reveal some details on the conditions under which the approximate Michaelis–Menten model is adequate. We base our analysis on a well defined criterion for the validity of the model as well as general analysis of its dynamics based on the ODE system. Cases where the biomass solutions of the two models are close, but the suggested dynamics differ significantly are also presented. The use of the proposed validation criterion and analysis for the approximate model are encouraged before any analysis of the process dynamics based solely on the values of constants such as the Michaelis constant K m . Section 4 considers the computation of the HMM-system rate parameters using time course data. In particular, we make use of some Mathematica and MATLAB tools allowing for the estimation of the three rate parameters of the HMM-system by means of appropriate adjusting simulations of the system solutions to match available time course measurement data. 2. Enzyme kinetics basic models 2.1. Henri’s reaction scheme The following reaction scheme of a simple enzyme–substrate dynamics where two fractions of the enzyme (free and bound) are involved has been proposed by Victor Henri [7,13,14]: k1 k2 −→ S+E C −→ P + E. ←− k−1
(1)
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
3
Henri’s reaction scheme (1) describes the reaction mechanism between an enzyme E with a single active site and a substrate S, forming reversibly an enzyme–substrate complex C, which then yields irreversibly a product P. Reaction scheme (1) says that during the transition of the substrate S into product P the enzyme E bounds the substrate into a complex C having different properties than the free enzyme and thus being necessarily considered as a separate substance. Denoting the concentrations s = [S], e = [E], c = [C], p = [P] and applying the Mass Action Law to Henri’s reaction scheme (1) we obtain the following system of ODEs: ds/dt = −k1 es + k−1 c, de/dt = −k1 es + (k−1 + k2 )c, dc/dt = k1 es − (k−1 + k2 )c,
(2)
d p/dt = k2 c, to be further called the HMM-system—in tribute to V. Henri, as well to L. Michaelis and M. Menten [21]. If the three rate constants k1 , k−1 , k2 are known, system (2) can be treated as a Cauchy ODE problem with initial conditions s(0) = s0 > 0, e(0) = e0 > 0, c(0) = 0, p(0) = 0. In practice the HMM-system rate constants k1 , k−1 , k2 are often not known and have to be determined for every specific pair enzyme–substrate. The contemporary approach to this task is to consider the rate constants as parameters in the dynamic HMM-system (2) and to compute them by fitting the solutions of the system to time course experimentally measured data. This problem will be considered in Section 4. 2.2. Michaelis–Menten equation Applying the quasi-steady-state assumption (QSSA) to HMM-system (2) one derives the following approximate equations for the substrate rate ds/dt and the complex c [5,21,22]: e0 s . Km + s ds Vmax s =− = −µ(s) dt Km + s
c=
(3) (4)
where the right-hand side function µ(s) = Vmax s/(K m + s) is known as “specific growth function” in cell growth modeling [19]. The solution sm of Eq. (4) is an approximation of the solution s of (2). The QSSA is reasonably applied e.g. when the ratio [E]/[S] is small so that fermentation continues for a considerably long time interval during which the concentration c = [C] of the bound enzyme is nearly constant. The latter can be achieved e.g. when the condition e0 ≪ s0 holds. Let us also remind that under the validity of QSSA the complex concentration c = [E S] is expressed in terms of the MM-substrate concentration via (3). In their paper [21] Michaelis and Menten discuss in detail Henri’s reaction scheme and Eq. (4) known as the Michaelis–Menten equation to be further denoted as MM-model. In addition Michaelis and Menten proposed a protocol for the practical calculation of the constant Vmax and K m in (4). The constant K m is known as Michaelis constant [7,30]. The MM-equation (4) can be written in the form ds Vmax s Vmax =− =− , dt Km + s K m /s + 1 showing that for large values of s the right-hand side is close to the constant −Vmax ; hence the substrate uptake rate |ds/dt| is nearly equal to the constant Vmax . Eq. (4) produces a good approximation sm for s under certain conditions (sometimes called validity criteria) [4,5,10,12,18,28,25,26]. Thus the condition e0 ≪ s0 assures good approximation and is ubiquitous for many fermentation and biotechnological processes, but may not be present, e.g. in living cells [1,29]. In the sequel we discuss certain methods and tools for the computation of the Michaelis constant using both the MM-model (4) and the HMM-system (2). In particular, we are interested in computational methods for the evaluation of the three rate parameters of the HMM-system.
4
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
1.0
0.8
0.6
0.4
0.2
1
2
3
4
Fig. 1. Graphics of the substrate dynamics according to the HMM-system (2) and MM-model (4). The rate constants of (2) are k1 = 2.62, k−1 = 0.1, k2 = 1.25, initial conditions: s0 = 1, e0 = 1.5 (the system is dimensionless); the parameters of (4) are K m = (k−1 + k2 )/k1 = 0.51526, Vmax = k2 e0 = 1.875.
2.3. Relations between the rate parameters In Fig. 1 the substrate uptake s is visualized in two different ways. The two graphics present the approximate solution sm to the MM-model (4) as well as the original solution s of the HMM-system (2). In this example the kinetic parameters are chosen so that the validity conditions of the MM-model do not hold and the difference of the solutions for the substrate of the two systems is visible. We would also want to compare the two solutions in case the validity conditions of the approximate model do hold. In order to correctly compare the solutions s and sm one has to establish certain consistency relations between the parameters in the MM-equation and the HMM-system. From the derivation of the MM-equation, see e.g. [22], we know that the parameters in (4) can be expressed by the parameters in (2) by means of the relations: Vmax = k2 e0 ,
K m = (k−1 + k2 )/k1 .
(5)
On the other side the three parameters in (2) cannot be uniquely determined by the two parameters in (4). The following computational examples show how the approximate substrate concentration solution sm to the MM-equation may look like depending on the initial values of the substrate s0 and the enzyme e0 . According to the MM-model we have ds/dt = Vmax s/(K m + s). The two constants Vmax and K m can be experimentally determined according to the MM-protocol. If the MM-model is valid we have Vmax ≈ k2 e0 . From this expression one can calculate the value kcat = Vmax /e0 . Hence kcat ≈ k2 if the validity conditions for the MMmodel take place. Assuming that k2 is a universal rate constant for a given specific reaction, then it should not depend on whether the validity conditions take place or not (e.g. if k2 is computed in an experiment with e0 ≪ s0 , it should be the same constant in another experiment with e0 > s0 .) Thus we may assume kcat ≈ k2 . This gives a method for evaluating the HMM-rate constant k2 . However, it is not possible to obtain any information about the HMM-rate parameters k1 , k−1 on the basis of the MM-model, resp. the MM-protocol. 3. Computational examples: model comparison In this section we present the results of several numerical studies aiming the comparison of the substrate dynamics of the two models (2), (4). The values of the dimensionless parameters K = K m /s0 and ε = e0 /s0 play an important role in our analysis because they can significantly influence the behavior of the system. The values of the two system characteristics K , ε are varied in the numerical study in order to present results for the cases when K ≤ 1 and ε ≪ 1, ε > 1, as well as K ≫ 1 and ε ≪ 1, ε > 1. The values of the HMM rate constants k1 , k−1 , k2 are kept the
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
5
same for all cases as follows: k1 = 0.1, k−1 = 0.01, k2 = 10. The Michaelis constant consistent with these values is according to (5): K m = (k−1 + k2 )/k1 = 100.1. 3.1. Criterion for good approximation An important requirement for the validity of the MM-model (4) is the following condition: k2 e0 e0 1 ≪ 1. = s0 + K m 1 + (k−1 /k2 ) + (s0 k1 /k2 ) k1 (s0 + K m )2
(6)
Condition (6) can be derived from the relation between the estimates of two timescales tc and ts , see [18], or [22, ch. 6.2]. The timescale tc corresponds to the “faster” reaction, involving the enzyme and the complex. Complex formation is the most significant process during this timescale. The longer timescale, ts , corresponds to the product synthesis and substrate utilization. It is natural to assume tc ≪ ts which in turn leads to (6). We see that e0 ≪ s0 implies (6) but in case when e0 ≪ s0 does not hold, the MM-model could still be a good approximation if K m is large enough (K m ≫ s0 and also k−1 ≪ k2 ). In order to provide a more comprehensive explanation of this statement we should consider the non-dimensionalized system proposed in [22], as well as the form of the MM-model under the considered conditions: τ = k1 e0 t, λ = k2 /(k1 s0 ),
u(τ ) = s(t)/s0 , K = K m /s0 ,
v(τ ) = c(t)/e0 , ε = e0 /s0 ,
du/dτ = −u + (u + K − λ)v, εdv/dτ = u − (u + K )v,
u(0) = 1, v(0) = 0.
(7) (8)
First, the dynamics of v can be summarized as follows: v increases from t = 0 until it reaches its maximum at v = u/(u + K ), after that it decreases back to 0. An important point is that if K is “big” then the maximum of v will be smaller. Furthermore, larger K indicates either larger k2 and/or larger k−1 values in comparison to s0 k1 . The most common ranges of the values of k2 and k1 are considered to be k2 ∈ [103 , 106 ]s −1 , k1 ∈ [107 , 1010 ]M−1 s −1 from many sources (such as [2,31,33]). Some sources [31,24] also state that for “superefficient enzymes” (or enzymes with perfect kinetics) the rate-limiting step is the substrate–enzyme association step. Such enzymes are said to have reached kinetic perfection or catalytic perfection [32] and that k2 ≫ k−1 holds for them [31]. Additionally, in certain rare cases the concentration of substrate is very small in comparison to the available enzyme [1]. Such cases may include formation of cellular micro-compartments, enzymatic activity in mammalian muscle tissue and others [1]. We may assume that in these cases the concentration of substrate can even reach values s0 ≪ k2 /k1 . We will assume the following hold for some of the cases we study: (1) the substrate concentration is sufficiently low, so that k2 ≫ s0 k1 , (2) the formed enzyme–substrate complex breaks down to product and enzyme faster than it dissociates to substrate and enzyme, thus k2 ≫ k−1 . For the enzymes under consideration we can summarize that the rate-limiting step of their reactions is the substrate–enzyme association step and due to their catalytic perfection a small amount of complex is present during the whole process. It should be noted that for the case studies in the next section we have used values of the kinetic parameters different than the ones derived from real experiments for convenience and simplicity. However, the conditions we rely on (k2 ≫ s0 k1 , k2 ≫ k−1 ) are kept valid when examining the case K ≫ 1 and we argue that parameters close to the ranges observed in real experiments exist, based on the aforementioned statements. We could also rewrite the Michaelis–Menten model (4) under the assumptions K ≫ 1, k−1 ≪ k2 in the following form: Vmax s k2 e0 s k 1 k 2 e0 s ds =− = − k +k =− ∼ −k1 e0 s. −1 2 dt Km + s k−1 + k2 + k1 s +s k1
The obtained form suggests a much simpler reaction and a negligible role of the intermediate complex. Thus, the product formation is only limited by the kinetic constant k1 . The derived equation is also identical to the exact substrate ODE from (2) if we assume e = e0 , c ∼ 0, which would be natural given the kinetic parameter assumptions we have made.
6
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
3.2. The exponential decay model We may argue that under the conditions discussed above the second terms in the non-dimensionalized system (7)–(8) are almost 0. System (7)–(8) can be approximated as follows: du/dτ = −u,
dv/dτ =
u ; ε
u(0) = 1, v(0) = 0.
(9)
The solution u of this dynamical system is identical to the solution sd of the substrate dynamics induced by the simple reaction scheme: k
S + E −→ P + E,
(10)
wherein S is the substrate, E is the enzyme and P is the product. Reaction scheme (10) says, differently to (1), that during the transition of the substrate S into product P the enzyme E does not bound the substrate into a complex. Applying the mass action law, assuming e = [E] constant, the kinetic scheme (10) leads to the following “exponential decay” differential equation for the substrate concentration sd = [S]: dsd = −kesd . dt
(11)
Solving the Cauchy problem related to Eq. (11), having set e = e0 , we obtain k = k1 . Indeed, the solution of Eq. (11) is sd (t) = s0 e−e0 kt . If we then solve the Cauchy problem (9), we would obtain s(t) = s0 e−e0 k1 t after converting it back to the initial variables. Our analysis shows that model (11) deserves to be included in the comparison of the two models (2), (4). We shall further show that under the conditions K ≫ 1, k−1 ≪ k2 the three models (2), (4) and (11) do indeed provide almost identical results and that the substrate solution of the MM-model may provide a good approximation even though e0 ≪ s0 does not necessarily hold. 3.3. Computational examples: three cases Case 1. K ≈ 1 (s0 ≈ K m ) (also valid for the range K ≤ 1 (s0 ≥ K m )). Computational example 1.1. In this example the MM–HMM-solutions for the substrate are close when ε is “small” as seen in Fig. 2. The values of the initial conditions used in the models are as follows: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0. The values of the parameters in (4) are consistent with the initial conditions and rate parameters in (2) and are calculated using (5) as Vmax = k2 e0 , K m = (k−1 + k2 )/k1 . The parameter k in the “exponential decay” differential equation (11) is set to be equal to k1 as derived in the previous section. In all the examples that follow, only s0 and e0 will vary, all the other model parameters will have the same values as in this example or they will be calculated using the corresponding s0 , e0 . This numerical example demonstrates that when the condition ε ≪ 1 holds then the MM-model (4) can be a very good approximation of the original HMM-system (2). Note that in this example the Michaelis constant K m used for the computation of the approximate MM-solution is derived from the HMM-rate parameters k1 , k−1 , k2 . This means that both models describe the process dynamics using equivalent rate constants and, since we consider the HMM-system to be true, this implies that the approximate model is also valid under the assumption ε ≪ 1. Our next numerical examples aim to examine what happens whenever the assumption ε ≪ 1 does not hold. Computational example 1.2. Our second numerical example shows how the solutions start to deviate when s0 and e0 are close to each other (Fig. 3). We can observe that with ε ∼ 1 the transient phase, related to the timescale tc , takes much longer and the assumption that tc ≪ ts is violated, thus invalidating the use of the MM-model. The goodness of the approximation of the complex c can be a better indicator of the applicability of the MM-model than the approximation of the substrate s, as c undergoes a boundary layer effect in the beginning of the reaction. The values of the initial conditions used are as follows: s0 = 100, e0 = 50, c0 = 0, p0 = 0. As before, we use (5) to calculate Vmax = k2 e0 , K m = (k−1 + k2 )/k1 .
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
7
Fig. 2. The solutions of the substrate (for models (2), (4), (11)) and the complex (for models (2), (4)) in the case ε = 10−3 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 +k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0. Under these conditions the MM-model clearly serves as a good approximation of the exact HMM-system.
Computational example 1.3. In this numerical example ε > 1. The HMM-solution for the substrate is even closer to the exponential decay solution (Fig. 4). As observed from the approximations of the complex in the same figure, with ε > 1, the MM-model is inapplicable due to the much longer timescale of the complex buildup. The values of the initial conditions used for the given solution are as follows: s0 = 100, e0 = 400, c0 = 0, p0 = 0. Here, and everywhere in the sequel Vmax = k2 e0 , K m = (k−1 + k2 )/k1 using (5). A figure of the specific growth rate function µ(s) = Vmax s/(K m + s) is included for reference, see Fig. 5. Case 2. K > 1. Computational example 2.1. It is clear that for values of K larger than 1 (K ≈ 2 in this example) but not K ≫ 1, the relations between the substrate solutions of the three models (2), (4) and (11) are preserved (Fig. 6). The values of the initial conditions used in the models are as follows: s0 = 50, e0 = 0.1, c0 = 0, p0 = 0. Computational example 2.2. In the second numerical example of this case the substrate solutions deviate when s0 and e0 are close to each other (Fig. 7), just as in example 1.2 and we can draw the same conclusion from the behavior of the complex, namely the assumption that tc ≪ ts has been violated, invalidating the use of the MM-model. All these follow from our choice of s0 and e0 (ε = 1). The values of the initial conditions used for the given solution are as follows: s0 = 50, e0 = 50, c0 = 0, p0 = 0. Computational example 2.3. In this numerical example ε > 1 (Fig. 8). The solutions are identical to the solutions of example 1.3. The approximation of the complex c in Fig. 8 can lead us to the conclusion that the MM-model is
8
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Fig. 3. The solutions of the substrate (for models (2), (4), (11)) and the complex (for models (2), (4)) in the case ε = 0.5 ∼ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 500; initial conditions: s0 = 100, e0 = 50, c0 = 0, p0 = 0.
inapplicable due to the much longer timescale of the complex buildup. This again follows from our choice of initial conditions s0 = 50, e0 = 150, meaning ε = 3 > 1. The values of the initial conditions used for the given solution are as follows: s0 = 50, e0 = 150, c0 = 0, p0 = 0. A figure of the specific growth rate function µ(s) = Vmax s/(K m + s) is included for reference, see Fig. 9. Case 3. K ≫ 1. Computational example 3.1. In this numerical example ε ≪ 1 and the solutions of the three models (2), (4) and (11) are identical (Fig. 10). The values of the initial conditions used for the given solution are as follows: s0 = 0.1, e0 = 0.001, c0 = 0, p0 = 0. Computational example 3.2. Although we have set the initial conditions so that ε ≫ 1, the substrate solutions of the three models (2), (4) and (11) are still very close to each other (Fig. 11). The values of the initial conditions used in the models are as follows: s0 = 0.1, e0 = 2000, c0 = 0, p0 = 0. A figure of the specific growth rate function µ(s) = Vmax s/(K m + s) is included for reference, cf. Fig. 12. As described above, in case 3: K ≫ 1 (also keeping in mind the assumption k−1 ≪ k2 ), the MM-model can serve as a very good approximation to the Henri–Michaelis–Menten law in terms of substrate dynamics, regardless of the range of initial conditions. As seen from the solutions of the complex though (Fig. 11), the MM-model completely fails to provide a good approximation of the complex dynamics when ε ≫ 1. This is expected since the approximate model only provides the outer solutions of the system in terms of singular perturbation analysis and in this case the
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
9
Fig. 4. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 4 > 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 4000; initial conditions: s0 = 100, e0 = 400, c0 = 0, p0 = 0.
Fig. 5. The specific growth rate function µ(s) = Vmax s/(K m + s) from (4) calculated in the range s from 0 to s0 = 100. Case—ε = 10−3 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0.
10
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Fig. 6. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 0.002 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 50, e0 = 0.1, c0 = 0, p0 = 0.
inner solution has to be accounted for due to the longer timescale of the complex buildup process. According to the MM-model, the solution of the enzyme–substrate complex can be derived as c(t) = e0 s(t)/(K m + s(t)), which under the proposed conditions (K ≫ 1, k−1 ≪ k2 and ε ≫ 1) cannot serve as a good approximation of the exact complex behavior. If we consider the non-dimensionalized system (9) under the same conditions, we can conclude that dv/dτ ∼ 0 hence v = const. The specific growth rate function visualized on Fig. 12 is almost linear in this case (again, regardless of whether ε ≪ 1 or ε ≫ 1) in contrast to cases 1 and 2. It is thus natural to expect the approximate model to predict substrate solutions similar to exponential decay in this case. Finally, the criterion (6) is obviously violated for computational example 3.2. (in contrast to all other examples where it is evaluated to less than 1): k2 e0 = 19.9. k1 (s0 + K m )2 Nevertheless, since K ≫ 1, k−1 ≪ k2 , the proposed explanations in the beginning of the section are indeed manifested and the substrate solutions of all three considered models are identical. As we also noted, this does not necessarily mean that the solutions of the other molecular species are correct, e.g. the MM-model solutions for the complex are incorrect in case ε ≫ 1.
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
11
Complex
Fig. 7. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 500; initial conditions: s0 = 50, e0 = 50, c0 = 0, p0 = 0.
Conclusion. There exist theoretical cases when the decay model produces better approximation to the exact HMMscheme than the MM-model. Additionally, in certain cases in which the MM-model can provide good substrate approximations, it may completely fail to give correct estimates of the complex or enzyme behavior. 3.4. Is the Michaelis constant a reliable kinetic parameter? Examples 1.2., 1.3., 2.2. and 2.3. demonstrate that the approximate model’s solutions are far from those of the exact HMM-system despite of the consistency of the parameters used in (2) and (4) even though we derive the Vmax and K m constants from the solutions of the exact model. This, of course, follows from the fact that under the initial conditions we have set the approximations embedded in the MM-model do not hold. However, if we free the MMmodel parameters from the conditions binding them to the exact scheme (5), we may still be able to find a solution closer to the exact one (compare Figs. 13(a) and 14 to Figures resp. 4 and 8). On the other hand, even if such a solution exists, it would correspond to completely different dynamics as we could easily determine from the solutions of the complex c in Fig. 13(b). The duration of the transient phase, as observed in Fig. 13(b), is related to the parameter ε (ε > 1 in the considered case) and it indicates that the MM-model may not be suitable under the given conditions since the approximations embedded in it are in contradiction to the behavior of the system. As a consequence any parameter estimation would suggest values of the kinetic parameters much different than the real ones. Although other sources [26] have also revealed the flaws of the inappropriate application of such approximate models under conditions which do not hold, some authors still prefer to fit their experimental data with them. As in Fig. 13(a) the solutions may appear to be close to those of the exact scheme but the dynamics that the approximate
12
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
complex
Fig. 8. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 3 > 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1500; initial conditions: s0 = 50, e0 = 150, c0 = 0, p0 = 0.
Fig. 9. The specific growth rate function µ(s) = Vmax s/(K m + s) from (4) calculated in the range s from 0 to s0 = 50. Case — ε = 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 500; initial conditions: s0 = 50, e0 = 50, c0 = 0, p0 = 0.
model suggests could be completely different than the real one (i.e. the MM-model parameters used in Figs. 13(a), 14 are respectively (K m = 200, Vmax = 9000) for Case 1 and (K m = 85, Vmax = 1500) for Case 2; in both cases K m can be evaluated to K m = 100.1 while Vmax = 4000 in Case 1 and Vmax = 1500 in Case 2).
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
13
Complex 1.X10-6
3.X10-7
5.X10-7
4.X10-7
2.X10-7
10000
20000
30000
40000
50000
Time
Fig. 10. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 0.01 ≪ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 0.01; initial conditions: s0 = 0.1, e0 = 0.001, c0 = 0, p0 = 0.
As explained in the previous section it is always beneficial to verify the criterion for the validity of the approximate model before applying it. Some general analysis of the behavior of the ODE system, based on the initial conditions, could also be valuable. We can similarly find relatively good fits of the decay model (11) for Examples 1.1. and 2.1. if we unbind the rate constant k. Figs. 15 and 16 show that in some cases the exponential decay solution may still approximate the system behavior to a certain degree, although the rate constant k of the decay model may not be strictly connected to the parameters of the exact HMM-model. Of course, if ε is small enough so that the substrate uptake is almost linear, no solution of (10) would be close to the exact solutions produced by (2). Conclusions. The presented computational examples show that the concept of “reliability” of the approximate models (the MM-model and the decay model) is not well defined. This concept obtains different meanings depending on whether one keeps the parameters of the approximate models bound to the HMM-system parameters, or not. 4. Computation of the HMM-system rate parameters using time course data 4.1. General settings The dynamics of the following biochemical process is studied using the HMM-system. The evaluation of the models is based on the goodness of fit achieved with them, keeping in mind the measurement errors contained in the data. The correctness of the dynamics suggested by the model is assessed by deriving the Michaelis constant for
14
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Time Complex
Time
Fig. 11. The solutions of the substrate (for models (4), (2), (11)) and the complex (for models (4), (2)) in the case ε = 2 ∗ 104 ≫ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 2 ∗ 104 ; initial conditions: s0 = 0.1, e0 = 2000, c0 = 0, p0 = 0.
Fig. 12. The specific growth rate function µ(s) = Vmax s/(K m +s) from (4) calculated in the range s from 0 to s0 = 0.1. Case — ε = 2∗104 ≫ 1. Kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 2 ∗ 104 ; initial conditions: s0 = 0.1, e0 = 2000, c0 = 0, p0 = 0.
each experiment and comparing it to known values from the literature corresponding to the same biochemical process. The conclusion is that approximate models should be avoided as much as possible. In the next section we propose
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
15
Substrate concentration
Time
(a) A substrate solution of the MM-model closer to the exact one may exist if we free the MM-model parameters from the conditions binding them to the exact scheme (5). Complex
120 100 80 60 40 20
0.02
0.04
0.06
0.08
0.10
0.12
0.14
Time
(b) The poor approximation of the complex indicates that the MM-model may not be applicable in the presented case. Fig. 13. The solutions are as in Example 1.3 (ε = 4 > 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 4000; initial conditions: s0 = 100, e0 = 400, c0 = 0, p0 = 0), with the only difference in (4) being: K m = 200, Vmax = 9000.
contemporary tools allowing the computation of the HMM-system parameters directly from time course measurement data. The process we follow is similar to the parameter fitting procedures and reverse engineering described in other papers, such as [26,27,23]. 4.2. Case study: fitting the HMM-system against experimental data of acetylcholine hydrolysis We fit the HMM model to experimental data of acetylcholine hydrolysis by acetylcholinesterase that has been obtained using the derived rate constants from the HXA method of the experiments in [34]. The rate constants were used to compute the solutions of (2) and evaluate them for a number of time points ti . Finally, a certain amount of noise was added to the data according to the standard deviations given in [34]. We take into account the measurement errors in the data and they are displayed as concentration intervals for each point in the plots. Our computational problem can be formulated as follows. Given time course experimental data (together with measurement errors) for the substrate, enzyme and product concentrations find values for the parameters k−1 , k1 , k2 and initial values E 0 , s0 such that the solution of the HMM system (2) fit well against the experimental data and
16
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
Time
Fig. 14. The solutions are as in Case 2, Example 3 (ε = 3 > 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1500; initial conditions: s0 = 50, e0 = 150, c0 = 0, p0 = 0), with the only difference in (4) being: K m = 85. Substrate concentration
Time
Fig. 15. The solutions are as in Case 1, Example 1 (ε = 10−3 ≪ 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m =
(k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 100, e0 = 0.1, c0 = 0, p0 = 0), with the only difference in (11) being: k = 0.07, rather than k = 0.1.
possibly fit into the measurement intervals. Our procedure for solving this problem passes in two stages: first we find an initial rough “guess” for the parameter values, then we iteratively improve the parameter set, resp. the solutions, using some optimization routines of the numerical computing environment MATLAB. More precisely, we have used MATLAB’s lsqnonlin procedure (optimization algorithm defaults to Trust-Region-Reflective Algorithm), where we aim to minimize the difference between the experimental data and the computed enzyme kinetics model solutions (using ode23 or ode23s solvers) for a given set of rate constants and initial conditions. At the first step, the solution of the ODEs HMM-system (2) with initial “guess” parameters (Fig. 17) is shown in Fig. 18. The values of the parameters are as follows: k1 = 25000M−1 s −1 , k−1 = 0s −1 , k2 = 10s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M. The results of the iterative optimization procedure we applied using the initial parameters can be seen in Fig. 19. The values of the parameters used for the given solution are as follows: k1 = 16847M−1 s −1 , k−1 = 7s −1 , k2 = 12s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M. Using the parameter set we have found to solve the HMM-system we get solutions that fit nicely the experimental data we started with. Using K m = (k−1 +k2 )/k1 we can derive the Michaelis constant from the optimal parameters in the HMM-system, K m = 0.00112M. This value is close to known values of the Michaelis constant from other sources [9,11], hence we can conclude that the model estimates correctly the behavior of the underlying biochemical process.
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
17
Substrate concentration 50
40
30
20
10
0
0
100
200
300 Time
400
500
600
Fig. 16. The solutions are as in Case 2, Example 1 (ε = 0.002 ≪ 1; kinetic parameter values: k1 = 0.1, k−1 = 0.01, k2 = 10, K m = (k−1 + k2 )/k1 = 100.1, Vmax = k2 e0 = 1; initial conditions: s0 = 50, e0 = 0.1, c0 = 0, p0 = 0), with the only difference in (11) being: k = 0.08, rather than k = 0.1.
Fig. 17. An example of the usage of dynamical numerical solution of the ODE system and their graphical representation using the Mani pulate and N DSolve functions of CAS Mathematica. An initial rough fit of the experimental data is obtained in this example.
The MM-model (4) on the other hand provides very different results for the same experimental data (Fig. 20). Although the obtained parameters (Vmax = 7.12 ∗ 10−7 Ms −1 ; K m = 2.193 ∗ 10−3 M) are relatively close to the real ones, the solutions do not pass through the intervals of the data. As in Example 2 and Example 3 of the first two Cases, the approximate model’s solutions deviate from the real ones even though its parameters may be close to the exact model’s rate constants. This can easily be explained by the fact that the complex c has more complicated dynamics (see Fig. 19 for the dynamics of c and e) than what the MM-model can estimate, i.e. data points from the transient phase cannot be accounted for. The obtained parameters are thus unreliable as they do not correspond to a good approximation of the complex c and the enzyme e. A rather good fit using the same approximate model (4) can be achieved by omitting the first observation data point of the enzyme (Fig. 21). The values of the model parameters (Vmax = 6.6 ∗ 10−7 Ms −1 ; K m = 1.26 ∗ 10−3 M) are even closer to those found using the HMM model. In general, we cannot expect data points related to the “faster” reaction (taking place during the tc timescale) to be approximated well with the MM-model because it is entirely focused on the “slower” reaction and the solutions we get cannot possibly follow the dynamics during both timescales.
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation ( Enzyme Kinetics Model
Concentration (substrate,product [mM])
x10-3
)
–
x10-7 substrate - s product - p experimental data - s experimental data - p enzyme e complex c experimental data - e errors
3
2 0.6
0.4 1
concentration - enzyme, complex [mM]
18
0.2
0
0
2000
4000
6000
8000
10000
12000
0 14000
Fig. 18. The numerical solution of the HMM-system with an initial parameter “guess” set of the unknown parameters k1 = 25000M−1 s −1 , k−1 = 0s −1 , k2 = 10s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M. Enzyme Kinetics Model
x10-7 substrate - s product - p experimental data - s experimental data - p enzyme e complex c experimental data - e errors
2.5
2 0.6 1.5
0.5 0.4
1 0.3
concentration - enzyme, complex [mM]
Concentration (substrate,product [mM])
x10-3
0.2
0.5
0.1 0
0
2000
4000
6000
8000
10000
12000
14000
Fig. 19. The numerical solution of the HMM-system (2) after the fitting of the model. The parameters and initial values used for the given solution are k1 = 16847M−1 s −1 , k−1 = 7s −1 , k2 = 12s −1 ; s0 = 2.5 ∗ 10−3 M, e0 = 5.4 ∗ 10−8 M.
5. Concluding remarks In this work we demonstrate that the value of the Michaelis constant may strongly depend on the way one makes use of the experimental data and the method of calculation. On the other side the rate constants in the enzyme kinetic Henri’s reaction scheme are well defined and can be efficiently computed from time course experimental data. Moreover, the modeling of metabolic pathways makes use of multiple reactions and can be made more rigorous when based on Henri’s reactions instead of Michael–Menten approximate reactions, being especially imprecise when applied to living cells. Such arguments confirm the need of refining the general methods for the determination of the
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
Enzyme Kinetics Model
x10-3
19
– x10-7
experimental data - s experimental data - p enzyme e complex c errors experimental data - e 2.5
2
0.7 0.6
1.5
0.5 0.4
1 0.3 0.2
0.5
0
0
2000
4000
6000
8000
10000
12000
concentration - enzyme, complex [mM]
Concentration (substrate,product [mM])
substrate - s product - p
14000
Fig. 20. The numerical solution of the MM-model after fitting it to the experimental data. Model parameters obtained from the optimization procedure: Vmax = 7.12 ∗ 10−7 Ms −1 ; K m = 2.193 ∗ 10−3 M. Enzyme Kinetics Model
x10-3
x10-7
experimental data - p enzyme e complex c errors experimental data - e 2.5
2 0.6 1.5
0.5 0.4
1 0.3
concentration - enzyme, complex [mM]
Concentration (substrate,product [mM])
substrate - s product - p experimental data - s
0.2
0.5
0.1 0 0
2000
4000
6000
8000
10000
12000
14000
Fig. 21. The numerical solution of the MM-model after fitting it to the experimental data. The first experimental data point for the enzyme has been omitted in the optimization. Model parameters obtained from the optimization procedure are now Vmax = 6.6 ∗ 10−7 Ms −1 ; K m = 1.26 ∗ 10−3 M.
Henri’s reaction scheme rate constants. Such general methods using contemporary software tools are suggested and discussed in the present work. In addition we wish to point out the need of validation methods in combination with the availability of experimental data involving measurement errors [17]; a more detailed introduction of the use of ⁀ validated interval methods is presented in [8]. In this work we make use of CAS Mathematica as computational tool, however other suitable software tools can be used, such as MATLAB and COPASI [15,20]. Ensuring the validity of the Michaelis–Menten model has been identified as a significant step in its application which has led to the comprehensive studies of validation criteria [28].
20
S. Dimitrov, S. Markov / Mathematics and Computers in Simulation (
)
–
When modeling kinetic reactions we suggest to start with a description of the process by a biochemical reaction equations (schemes) and then pass to mathematical differential equations via the mass action principle. Such an approach is useful in keeping close to the biochemical mechanism of the modeling process; we recommend the use of this approach not only in fields like enzyme kinetic and metabolic networks but also in fields like cell growth and population dynamics [19]. Acknowledgments The authors are very grateful to the anonymous reviewers whose valuable reports contributed much to the quality of the paper. References [1] K.R. Albe, M.H. Butler, B.E. Wright, Cellular concentrations of enzymes and their substrates, J. Theoret. Biol. 143 (1990) 163–195. [2] M.N. Berberan-Santos, A general treatment of Henri–Michaelis–Menten enzyme kinetics: exact series solution and approximate analytical solutions, MATCH Commun. Math. Comput. Chem. 63 (2010) 283–318. [3] A.M. Bersani, E. Bersani, G. Dell’Acqua, M.G. Pedersen, New trends and perspectives in nonlinear intracellular dynamics: one century from Michaelis–Menten paper, Contin. Mech. Thermodyn. 27 (2014) 659–684. [4] A.M. Bersani, G. Dell’Acqua, Is there anything left to say on enzyme kinetic constants and QSSA? J. Math. Chem. 50 (2012) 335–344. [5] G.E. Briggs, J.B.S. Haldane, A note on the kinetic of enzyme action, Biochem. J. 19 (1925) 338–339. [6] W.W. Chen, M. Niepel, P.K. Sorger, Classic and contemporary approaches to modeling biochemical reactions, Genes. Dev. 24 (2010) 1861–1875. [7] U. Deichmann, S. Schuster, J.-P. Mazat, A. Cornish-Bowden, Commemorating the 1913 Michaelis–Menten paper Die Kinetik der Invertinwirkung: three perspectives, FEBS J. 281 (2014) 435–463. J.-P. Mazat, Part 3: pp.452–463. [8] S. Dimitrov, G. Velikova, V. Beschkov, S. Markov, On the numerical computation of enzyme kinetic parameters, Biomath Commun. 1 (2) (2014). [9] A. Forsberg, Kinetics for the inhibition of acetylcholinesterase from the electric eel by some organophosphates and carbamates, Eur. J. Biochem. 140 (1984) 153–156. [10] R. Grima, N.G. Walter, S. Schnell, Single molecule enzymology a` la Michaelis–Menten, FEBS J. 281 (2014) 518–530. [11] A. Hai, et al., Acetylcholinesterase-ISFET based system for the detection of acetylcholine and acetylcholinesterase inhibitors, Biosens. Bioelectron. 22 (2006) 605–612. [12] S.M. Hanson, S. Schnell, The reactant stationary approximation in enzyme kinetics, J. Phys. Chem. A 112 (2008) 8654–8658. [13] V. Henri, Recherches sur la loi de l’action de la sucrase, C. R. Hebd. Acad. Sci. 133 (1901) 891–899. [14] V. Henri, Lois g´en´erales de l’action des diastases, Librairie Scientifique A. Hermann, Paris, 1903. [15] S. Hoops, et al., COPASI—a complex pathway simulator, Bioinformatics 22 (24) (2006) 3067–3074. [16] K.A. Johnson, R.S. Goody, The original Michaelis constant: translation of the 1913 Michaelis–Menten paper, Biochemistry 50 (39) (2011) 8264–8269. [17] W. Kraemer, J. Wolff, v. Gudenberg (Eds.), Scientific Computing, Validated Numerics, Interval Methods, in: Proc. of SCAN-2000/Interval2000, Kluwer/Plenum, 2001. [18] K.J. Laidler, Theory of the transient phase in kinetics, with special reference to enzyme systems, Can. J. Chem. 33 (1955) 1614–1624. [19] S. Markov, Cell growth models using reaction schemes: batch cultivation, Biomath. 2 (2) (2013) 1312301–9. [20] P. Mendes, S. Hoops, S. Sahle, R. Gauges, J. Dada, U. Kummer, Computational modeling of biochemical networks using COPASI, Methods Mol. Biol. 500 (2009) 17–59. [21] L. Michaelis, M.L. Menten, Die Kinetik der Invertinwirkung, Biochem. Z. 49 (1913) 333–369. [22] J.D. Murray, Mathematical Biology: I. An Introduction, third ed., Springer, 2002. [23] N. Nikolova, K. Tenekedjiev, K. Kolev, Uses and misuses of progress curve analysis in enzyme kinetics, Cent. Eur. J. Biol. 3 (4) (2008) 345–350. [24] M.H.M. Olsson, P.E.M. Siegbahn, A. Warshel, Simulations of the large kinetic isotope effect and the temperature dependence of the hydrogen atom transfer in lipoxygenase, J. Am. Chem. Soc. 126 (9) (2004) 2820–2828. [25] M.G. Pedersen, A.M. Bersani, E. Bersani, The quasi steady-state approximation for fully competitive enzyme reactions, Bull. Math. Biol. 69 (2007) 433–457. [26] M.G. Pedersen, A.M. Bersani, E. Bersani, G. Cortese, The total quasi-steady-state approximation for complex enzyme reactions — a word of caution, Math. Comput. Simul. 79 (2008) 1010–1019. [27] T.D. Pollard, E.M. De La Cruz, Take advantage of time in your experiments: a guide to simple, informative kinetics assays, Mol. Biol. Cell. 24 (8) (2013) 1103–1110. [28] S. Schnell, Validity of the Michaelis–Menten equation — steady-state, or reactant stationary assumption: that is the question, FEBS J. 281 (2014) 464–472. [29] S. Schnell, P.K. Maini, Enzyme kinetics at high enzyme concentration, Bull. Math. Biol. 62 (2000) 483–499. [30] S. Schnell, P.K. Maini, A century of enzyme kinetics: Reliability of the K m and vmax estimates, Comments Theor. Biol. 8 (2003) 169–187. [31] M.E. Stroppolo, M. Falconi, A.M. Caccuri, A. Desideri, Superefficient enzymes, Cell. Mol. Life Sci. 58 (2001) 1451–1460. [32] L. Stryer, J.M. Berg, J.L. Tymoczko, Biochemistry, fifth ed., W H Freeman, New York, 2002. [33] I.B. Wilson, M.A. Harrison, Turnover number of acetylcholinesterase, J. Biol. Chem. 236 (8) (1961) 2292–2295. [34] P. Zdrazilova, et al., Kinetics of total enzymatic hydrolysis of acetylcholine and acetylthiocholine, Z. Nat.forsch. C. 61 (3–4) (2006) 289–294.