Estimating kinetic constants in the Michaelis–Menten model from one enzymatic assay using Approximate Bayesian Computation

The Michaelis–Menten equation is one of the most extensively used models in biochemistry for studying enzyme kinetics. However, this model requires at least a couple (e.g., eight or more) of measurements at different substrate concentrations to determine kinetic parameters. Here, we report the discovery of a novel tool for calculating kinetic constants in the Michaelis–Menten equation from only a single enzymatic assay. As a consequence, our method leads to reduced costs and time, primarily by lowering the amount of enzymes, since their isolation, storage and usage can be challenging when conducting research.

In biochemistry, the Michaelis-Menten model [1][2][3] is one of the best-known and useful approaches to enzyme kinetics [4][5][6][7]. It takes the form of an equation describing the rate of enzymatic reaction by relating a rate of formation of the product to a substrate concentration. The system involves two reactions where a substrate binds reversibly to the enzyme to form an enzyme-substrate complex, which then reacts irreversibly to generate a product and to regenerate the original enzyme. In the beginning of the early 20th century, Victor Henri discovered that enzyme reactions are initiated by a bond between the enzyme and the substrate [8]. Later, his work was taken up by Leonor Michaelis and Maud Menten [3], who measured the initial velocity as a function of sucrose concentration and derived an equation that approximates the modern version of the Michaelis-Menten equation, which is widely used in enzyme-substrate interaction study [9]. Moreover, this approach is applied to other biochemical problems, including antigen-antibody binding, protein-protein interaction [2] as well as pharmacokinetics [10], and to areas outside biochemical interaction like clearance of blood alcohol [11], alveolar clearance of dusts [12], and bacteriophage infection [13]. The model provides valuable knowledge for researchers in the form of kinetic parameters that explain the properties of a particular enzyme. The standard method for determining the kinetic parameters involves running a series of enzyme assays at multiple substrate concentrations, and measuring the initial rate of the reaction [2]. By plotting reaction rate against substrate concentration, and applying a nonlinear regression to the Michaelis-Menten equation, the kinetic parameters could be obtained. Until recently, the graphical methods, such as the one described by Lineweaver and Burk, which utilized the double reciprocal plot, were used for this purpose [2]. Nowadays, various software programs provide fast tools to obtain kinetic parameters from nonlinear regression of the Michaelis-Menten equation.
In order to minimize experimental work that involves several enzymatic assays, we propose to use the Approximate Bayesian Computation (ABC) techniques [14][15][16][17][18] to determine kinetic constant in the Michaelis-Menten model. Commonly, biochemical systems could be described using differential equations or a set of nondifferentiable set of rules with parameters having physical interpretation [6,19]. Then, for given initial conditions and specific values of parameters it is possible to generate data, that is, provide a solution of the set of differential equations. However, an analytical form of the distribution of solutions is highly complex and, thus, it is impossible to formulate the likelihood function. ABC techniques offer an almost automated solution in situations where evaluation of the posterior is computationally prohibitive, or whenever suitable likelihoods are not available. The main idea behind ABC is to utilize a mathematical model for a biochemical system and apply an approximate inference (e.g., Monte Carlo methods) to calculate the posterior over parameters. A classical example of applying ABC is a model of spreading tuberculosis [19] that contains a set of probabilistic rules of events, such as, illness transmission, mutation, death, and recovery. The model is relatively simple yet it is analytically intractable. Since the probabilistic inference is infeasible due to lack of an analytical form of the model, a rejection sampling procedure with a prespecified acceptance threshold was used to infer about parameters. Namely, values of parameters are first sampled from a prior distribution and then data are generated by the model. If the generated data are close to the real data in terms of a chosen distance measure (i.e., the distance is smaller than the threshold e), the parameters' values are accepted and further used to approximate the posterior. Later, the ABC techniques were widely used in other applications, such as, biology [20,21], evolution and ecology [15,22], the evolution of genomes [23], the dynamics of gene regulation [18], the demographic spread of species [24][25][26][27], or mRNA self-regulation [28]. It is worth mentioning that [28] used Parallel Tempering and Metropolis-Hastings sampling techniques to determine parameters of the Michaelis-Menten model. However, they used a synthetic likelihood function [29,30] assuming 1% Gaussian error while our approach is likelihood-free.
In this paper, we present a novel application of the ABC computational tool for calculating kinetic constants in the Michaelis-Menten equation using one enzymatic assay at a single substrate concentration. This extremely useful framework gives wider possibilities in determination of kinetic parameters by immense reduction of the cost, primarily by lowering the amount of enzyme, substrate and other reagents as well as time.

Implementation
The proposed method was implemented in PYTHON using the NumPy package (www.numpy.org). In all experiments, we used our own implementation of the Runge-Kutta (RK4) solver. The code is available at: https://github.com/ e-weglarz-tomczak/mmabc.
In the implementation of our method, we used two heuristics to obtain samples from the posterior distribution over parameters: • For computational convenience, the sampling procedure is run 100 times to sample 10 000 parameters. • If the procedure is stopped before reaching 100 iterations (runs), then we increase the threshold of the rejection sampling e and repeat the procedure.
We noticed that such approach worked very well in practice.

Parameter estimation from synthetic data
In order to verify the usefulness of the proposed approach, we first evaluate our method on synthetic data. We simulate 6000 different synthetic measurements using the following setting:

Parameter estimation from real data
The utility of the proposed method is further demonstrated on kinetic parameters measurement data for three real enzymes, namely human aminopeptidase (hAPN), Sus scrofa APN (ssAPN), and human endoplasmic reticulum aminopeptidase 2 (hERAP2) [31,32]. For each enzyme eight measurements were collected at different levels of substrate.
Concentration In order to obtain the values of the kinetic constants, the standard method was used [2]. Further details how the experiments were carried out are outlined in Refs [31,32].

Metrics
The following metrics are used to quantify the performance of the presented approach: • Deviation specifies how the mean estimated value, h est , deviates from the real one, h real : • Accuracy indicates how often the real value h real is within one standard deviation interval from the mean estimated across all simulations.

Results and Discussion
The Michaelis-Menten model The Michaelis-Menten model takes the form of an equation describing the rate of enzymatic reactions, by relating reaction rate to the concentration of a substrate. The enzymatic reaction in this model involves reversible reaction where an enzyme, E, binds to a substrate, S, to form a complex, ES, and irreversible releasing a product, P, and the free enzyme. This may be represented schematically as follows: By applying the law of mass action, which states that the rate of any chemical reaction is directly proportional to the product of the masses/concentrations of the reacting substances [33], system of four ordinary differential equations that define the rate of change of reactants over time are obtained, namely: In their original analysis, Michaelis and Menten assumed that the substrate is in instantaneous chemical equilibrium with the complex enzyme-substrate [3]. An alternative analysis was proposed 10 years later by Briggs and Haldane, known as quasi-steady-state of the system, that assumed the concentration of the intermediate complex is approximately steady during progression of the reaction [34]. Both approximations assume irreversible dissociation of ES complex to the product and free enzyme that yields the following equation of the velocity of the reaction: where P is the concentration of the product, S denotes the concentration of the substrate, V max is the maximal velocity, E is the initial concentration of the enzyme, k cat denotes the constant of the conversion to the product (the catalytic rate constant), and K M is called the Michaelis constant. Additionally, the initial value of the substrate, S 0 is known. The model provides crucial information about the nature of the enzyme in the form of the kinetic parameters. The Michaelis constant (K M ) is the concentration of substrate at which the reaction rate is half of the maximum rate. This is an inverse measure of the substrate's affinity for the enzyme. The smaller the value of K M , the higher the affinity; hence, the rate will approach maximal velocity with lower S than reactions with a lower K M . The catalytic rate (k cat ) shows how fast reaction is. The k cat /K M constant provides the knowledge how efficiently an enzyme converts a substrate into a product. The curve describing the relationship between the velocity and the substrate concentration in the Michaelis-Menten kinetics can be formulated as a nonlinear regression model. The commonly used model is expressed using the exponential function: where b is a real-valued parameter, or by using a polynomial of the k-th order: where a 0 , a 1 , . . . , a k are constants.

The standard approach to Michaelis-Menten kinetics
The typical method for determining the constants involves performing a series of in vitro experiments at varying substrate concentrations. Next, from enzyme assays, the initial reaction rates are obtained. The initial reaction rate is taken to mean because the equilibrium or quasi-steady-state approximation remains valid. By plotting reaction rate against concentration of the substrate, and using a nonlinear regression of the Michaelis-Menten equation, the kinetic parameters can be calculated. Nowadays, available software programs for enzyme kinetics calculations allow to fit the exponential Eqn (8) and/or the polynomial (9) to measured values. Before computational era, graphical methods involving linearization of the equation were used to perform nonlinear regression, such as, the Eadie-Hofstee diagram, Hanes-Woolf plot, and Lineweaver-Burk plot as most commonly used [2]. The standard approach for finding k cat and K M assumes the following three steps [2]: • Calculate the initial rate of the reaction measured at varying substrate concentration from the linear portion of the reaction progress curve (product vs. time) (see step 1 in Figure 1). • Fit the Eqn (8) or the polynomial (9) to the plotted reaction rate against concentration of the substrate (see step 2 in Figure 1). • Calculate the kinetic parameters K M and k cat from the fitted curve (see step 3 in Figure 1).
The procedure is straightforward and rather accurate [2]; however, it is prone to noise in measurements and requires multiple (e.g., eight or more) measurements at different substrate concentration level in order to properly fit the curve and, eventually, determine the kinetic constants. Since carrying out several enzymatic assays is a time-consuming process, the whole procedure becomes rather slow. Additionally, performing all experiments requires substantial amount of enzyme that is typically costly and time-demanding to obtain. As a result, costs and time needed to determine the kinetic constants are usually high and become a considerable bottleneck in the whole research process.

ABC approach to the Michaelis-Menten kinetics
One of the main inference frameworks is the Bayesian approach where the prior knowledge about parameters, h, is further updated by observations v (the posterior) [35]. According to Bayes' rule the posterior distribution over the model parameters, p(h|v) is proportional to the product of the likelihood function of the observed data, p(v|h), and the prior distribution over the parameters, p(h), namely: Typically, the likelihood function and the prior are chosen from a known family of distributions. However, in many cases, the form of the distribution p(v|h) remains unknown, but we can sample from it (i.e., generate new data). In the ABC literature, such model is called a simulator [15].
In the considered case, we treat the Michaelis-Menten model in (7) as the simulator of enzyme kinetics with parameters h = {k cat , K M }. For given values of the kinetic parameters k cat and K M sampled from some assumed prior (e.g., the uniform prior), the ODE can be solved using a numerical ODE solver (e.g., RK4 solver). As a result, simulated data v sim can be further compared with the real measurements v real . If a distance (e.g., the Euclidean distance, jj jj 2 2 ) between the simulated data and the real observations is smaller than a prespecified threshold e (e.g., e = 0.001), then we can keep the values of parameters. Otherwise, we disregard them. Eventually, we could use all the accepted parameters to estimate the posterior distribution. This procedure is called the ABC rejection algorithm [21].
Since the initial conditions are known (E 0 and S 0 ), we could easily run the simulator multiple times in order to estimate the posterior distribution over the kinetic constants. Moreover, the measurement could be noisy and we can still estimate the parameters by enlarging the threshold e. Therefore, we propose to utilize the ABC rejection method to determine the kinetic constants of the Michaelis-Menten model. For given e, E 0 and S 0 , the procedure is the following (Figure 2): • Sample h = {k cat , K M } from the prior p(h) being the uniform distribution over a prespecified closed region (e.g., k cat 2 [0, 100] and K M 2 [0, 200]). • Simulate data v sim using the Michaelis-Menten model and sampled parameters h in Step 1 by running the RK4 solver. • For all simulated data in Step 2, if jjv real À v sim jj 2 2 \e, then accept h that was used to simulate v sim .
• Approximate the posterior distribution over the kinetic constants using all accepted hs.
The benefit of the proposed approach is apparent for two practical reasons. First, instead of collecting multiple measurements, only one enzymatic assay is sufficient. As a result, the whole process can be sped up a couple of times. Second, there is a significant saving of used enzyme and substrate for determining the kinetic constants. This results in lower costs and shorter time spent on performing experiments. Further, our approach is based on a probabilistic method; therefore, we are able to provide the mean value with an uncertainty estimate, for example, by using the standard deviation.

Discussion
In our first experiment, we aimed at verifying whether it is possible to correctly determine the kinetic constants using a single observation. For this purpose, we generated synthetic data as described in Materials and methods. In order to quantify the performance of our approach, we used the deviation and accuracy. We expressed both metrics as a percentage. Ideally, the deviation should be equal 0% while the accuracy should be equal 100%.
The results of simulations (1000 simulations per single value of S 0 ) are presented in Figure 3 and Table A1 (see Appendix 1). First we noticed that the average deviation for k cat is low (~1-2%).Similarly, for K M , the deviation is also low (~5-6%). These results indicate that the proposed probabilistic approach allows to obtain good estimates of the real parameters. Moreover, in more than 90% of cases, the  real value is within one standard deviation from the estimated mean value. The remaining cases might be concerning; however, the deviation metric indicates that the uncertainty interval was simply very narrow and this could be alleviated by taking slightly larger e. We also experimented with taking three standard deviations and then the accuracy was 100%. The obtained results confirm that it is indeed possible to find very accurate estimate of the kinetic constants using a single enzymatic assay. Our second observation from this experiment is that there are preferred values of S 0 between and where our method attains the best scores of both metrics for k cat and K M . We believe that this effect could be explained as follows. If we look at the Michaelis-Menten plot (Figure 1), taking too large value of S would result in a value of the rate close to the maximal rate V max . In such case we can easily determine k cat , but we have less information about K M . On the other hand, taking too small value of S 0 would be not sufficiently informative in respect to k cat that is determined by the V max /2. The ideal value would be around the substrate concentration for which the rate is close to V max . This tells us most about K M and k cat . Obviously, finding such value of S might be challenging; however, the values around 100-200 seem to be sufficient in most cases. We will take advantage of this fact in the experiment on real enzymes.   We also inspected whether there is a dependency between a range of values of k cat and K M and the considered metric. We divided values of both parameters into four bins each and checked whether obtained values of metrics depend on the values of the parameters. Our conclusion is that there is no evident dependency between the values of parameters and the values of the metrics.
In the second experiment, we used data containing the progression of reaction catalyzed by enzyme from our previous study [31,32]. We considered a single substrate concentration assay for the ABC. The data involved experiments on isolated proteolytic enzymes for which the assumption of irreversibility of released product is fulfilled. We evaluated our approach on data from three real enzymes: hAPN, ssAPN, and hERAP2 [31,32].
The results for our method and the standard approach are presented in Figure 4 and Table A2 (see Appendix 1). In our method, we utilized a single enzymatic assay for S 0 = 200 since for this value we got the best results in the previous experiment. Interestingly, the results obtained using the standard approach and our approach do not differ too much. However, it is a known fact that the standard approach tends to under-or overestimate the true value [2] and, thus, it should be treated as a reference point rather than a real value. Nevertheless, the differences between both methods are rather small that allows us to conclude that our approach provides good estimates. Importantly, our method needs only one enzymatic assay while the standard procedure requires at least couple (e.g., eight) of them.

Conclusions
Application of statistical computational methods to biochemistry, biology, and chemistry allows to understand the underlying phenomena but they could also assist researchers in their daily routine. In this paper, we proposed a probabilistic method for determining the kinetic constants in the Michaelis-Menten model.
The advantage of applying the ABC rejection algorithm to this biochemistry relevant model is the possibility of using only one enzymatic assay. The single enzymatic assay provides information about progression of the reaction catalyzed by enzyme at one initial substrate concentration and we claim that it is sufficient to properly estimate the kinetic parameters. In comparison to the standard procedure that requires at least several measurements, our method notably reduces the whole research process. The obtained results on both synthetic and real data provide empirical evidence in favor of our claims. In our future works, we aim at focusing on more sophisticated ABC procedures [16,17,36] that could improve uncertainty estimation and sampling efficiency. Another interesting topic for future research is to develop new quality measures of the ABC methods. Finally, encouraged by the obtained results, we find investigating other kinetics models as a challenging future research direction.