Despite their obvious advantages, classical parameteri-zations require large amounts of data to fit their parameters. Particularly, enzymes displaying complex reaction and regulatory (allosteric) mechanisms require a great number of parameters and are therefore often represented by approximate formulae, thereby facilitating the fitting but ignoring many real kinetic behaviours. Here, we show that full exploration of the plausible kinetic space for any enzyme can be achieved using sampling strategies provided a thermodynamically feasible parameterization is used. To this end, we developed a General Reaction Assembly and Sampling Platform (GRASP) capable of consistently parameterizing and sampling accurate kinetic models using minimal reference data. The former integrates the generalized MWC model and the elementary reaction formalism. By formulating the appropriate thermodynamic constraints, our framework enables parameterization of any oligomeric enzyme kinetics without sacrificing complexity or using simplifying assumptions. This thermodynamically safe parameterization relies on the definition of a reference state upon which feasible parameter sets can be efficiently sampled. Uniform sampling of the kinetics space enabled dissecting enzyme catalysis and revealing the impact of thermodynamics on reaction kinetics. Our analysis distinguished three reaction elasticity regions for common biochemical reactions: a steep linear region (O> AG, >-2 kJ/ mol), a transition region (-2> AG, >-20 kJ/mol) and a constant elasticity region (AG, <-20 kJ/mol). We also applied this framework to model more complex kinetic behaviours such as the monomeric cooperativity of the mammalian glucokinase and the ultrasensitive response of the phosphoenolpyruvate carboxylase of Escherichia coli. In both cases, our approach described appropriately not only the kinetic behaviour of these enzymes, but it also provided insights about the particular features underpinning the observed kinetics. Overall, this framework will enable systematic parameterization and sampling of enzymatic reactions.
Different frameworks have been proposed to parameterize enzymatic reactions using approximate expressions while maintaining thermodynamic consistency. Approximate expressions have been particularly sought and used, as kinetic expressions typically require large amounts of data to fit their parameters. The latter however ignores real kinetic behaviours and incurs in loss of generality. To overcome these limitations, here we present a novel framework GRASP for exploring the kinetic behaviour of enzymatic reactions under uncertainty based on parameter sampling. By formulating the appropriate thermodynamic constraints and using minimal biochemical reference data, our framework is capable of parameterizing the kinetics of any oligomeric enzyme without sacrificing complexity. To this task, we integrated the generalized MWC model and the elementary reaction formalism, providing a thermodynamically-safe and easy-to-sample parameterization. Application of our framework provided valuable insights into how reactions are regulated under non-equilibrium conditions. We also showed how our approach can be used to describe and understand complex kinetic behaviours of enzymes involved in key regulatory steps of cell metabolism. Overall, this framework enables systematic exploration of the feasible kinetic behaviour of enzymes.
The particular catalytic pattern of the enzyme dictates its mathematical representation, which can be obtained by solving equations for the enzyme intermediate concentrations . A quasi-steady-state assumption for these intermediates is commonly employed to this end , yielding a final expression function of microscopic rate constants and reactant concentrations. The whole process can be automated using King-Altman’s method . Some of the rate constants are related to the apparent equilibrium constant of the overall reaction by the Haldane relationships, directly linking kinetics with thermodynamics . Moreover the rate constants can be converted into macroscopic kinetic constants , which can be measured and estimated from subsequent enzymatic assays.
With more complex reaction mechanisms, the number of parameters becomes so large that a combination of model reduction and/ or computational sampling is required for study. Models can be reduced by introducing simplifying assumptions, e.g., steps being at equilibrium, lumped elementary steps, etc. Multiple expressions have been proposed to reasonably approximate kinetic behaviour using available data while maintaining essential thermodynamic consistency [7—12]; however they always incur some loss of generality. Computational sampling enables us to determine emergent enzyme properties from high dimensional parameter spaces. Several sampling strategies have been proposed, for example uniform parameter sampling within the S-formalism [13,14], elasticity sampling [15,16] and relative enzyme saturation sampling [17,18]. All use simplified kinetic expressions (loss of generality) and most ignore intrinsic thermodynamic constraints between kinetic parameters, hence they will sample infeasible parameter sets.
Even though mass inhibition/ activation can account for some degree of regulation, metabolic control is mostly achieved through allosteric and transcriptional regulatory interactions . Modelling allosteric behaviour requires the inclusion of conformational information, which enables description of both allosteric and cooperative interactions. The two most famous of such models are the symmetry model of Monod, Wyman and ChangeuX (MWC)  and the sequential model of Koshland, Némethy and Filmer (KNF) . Both models are based on the assumed equilibrium between two conformational states of the enzyme; a relaxed (R) and tense (T) state. The main difference between these theories rests in how conformational transitions proceed upon binding of ligands. While the MWC model demands maintenance of conformational symmetry, the KNF model does not and instead requires strict induced fit. Although these models can be considered as special cases of more general models, these generalizations have so far not proven to be useful. In fact, the MWC model is regarded the model of choice for describing allosteric and cooperative interactions . Moreover, it has been demonstrated that the MWC model can be cast in a convenient mathematical form , which can be combined with the elementary reaction formalism.
Parameterization combines the generalized MWC model for modelling the kinetics of oligomeric enzymes with the elementary reaction formalism for deriving thermodynamically consistent catalytic eXpressions. By employing a convenient normalization at the elementary reaction level , we are able to describe reaction kinetics provided a reference flux and the thermodynamic affinity of the reaction at the reference point. Using an accurate parameterization necessarily requires a large number of parameters. An advantage is that the parameterization retains all intrinsic thermodynamic constraints between kinetic parameters. We designed an efficient sampler producing thermodynamically consistent parameters obeying the principle of microscopic reversibility. Using an efficient Monte Carlo sampling technique that eXploits the shape of the sampling space, we ensure high parameter quality and low parameter rejection rates. The framework is demonstrated through exploration of the full kinetic space of reactions, assessing the impact of thermodynamic affinity, reaction molecularity and mechanisms, as well as modelling compleX kinetic behaviours.
The framework employs the MWC model, which provides the basis for modelling most cooperative and allosteric behaviours of multimeric enzymes . The classical formulation including detailed assumptions and limitations are described in 81 Text. Our framework is based on the recasting of the model developed by Popova and Sel’kov [23,25—27] , in which the velocity of reaction of any oligomeric enzyme is expressed as the product of two independent functions,
Equation 1 provides a general and simple interpretation of the kinetics of an oligomeric enzyme. Firstly, the shape of the catalytic function is determined only by the mechanism of elementary interactions between substrates and products with one active site of the enzyme (catalytic mechanism). Secondly, the regulatory function is invariant with respect to the action mechanism of the catalytic sites. Thus, if one has information about the conformational mechanism of the enzyme, Le. number of subunits, transitions, allosteric sites and effectors, and possesses an eXpression representing the catalytic mechanism for the conversion of substrates to products, the catalytic and regulatory functions can be written as follows , where vR and VT represent the velocity of reaction for the R and T conformations of the oligomer-ic enzyme, respectively (both states follow the same reaction mechanism as protomers are identical), and Q is a function that determines the current ratio between the R and T states (see later).
An important feature of this parameterization is that it enables inclusion of fundamental thermodynamic relations between kinetic parameters, as it is compatible with the elementary reaction formalism. Some of these relations are lost when using other parameterizations. A complete overview of the proposed framework is depicted in Fig. 1A. In the following, we present a systematic method for parameterizing and sampling thermodynamically consistent kinetics. Before considering compleX cooperative or allosteric mechanisms, we will consider a simple non-allosteric enzyme, i.e. n = 1 and L = 0 (no tense state) or L = 00 (no relaxed state), the resulting flux at any reference state is purely due to the catalytic function.
Every enzymatic reaction can be broken into simple reversible steps called elementary reactions. Using the laW of mass action, the rate of each elementary reaction is written as ,
Typically, the absolute values of the metabolite and total enzyme concentrations are not known (although physiological ranges can be estimated). To overcome this limitation, normalization of all the variables around a reference point (steady-state flux) is a convenient strategy. Following the scaling procedure employed by Tran et al. , normalization of these variables yields, The normalized metabolite concentrations are unitary at the reference point (5cref = 1), which is used extensively in developing an efficient sampling strategy.
In order to sample kinetics, we can in principle sample the rate constants directly. However, this leads to an inefficient sampler, where thermodynamic constraints can only be validated after sampling resulting in a high rejection rate. Instead, we design the sampling procedure to incorporate constraints directly without the need of assuming any particular distribution for the rate constants. Sampling enzyme intermediate abundances. It can be shown that Equation 6 can be written in a simple generic form for the computation of the rate parameters at the reference state (81 Text), Where vref is the vector of elementary fluxes, P(éref) is a diagonal matrix function of the en-elem zyme intermediate abundances vector éref and I; is a vector of rate constants. Notably, the enzyme intermediate abundances sums to one and can be readily sampled using appropriate probability distributions. Specifically, we seek to uniformly sample enzyme compleX abundances subject to This problem can be viewed as uniformly sampling from the surface of a p-dimensional simplex, which is equivalent as sampling from a multivariate Dirichlet distribution With hyper-parameter vector 1 . This distribution can be easily constructed and readily sampled using Gamma distributions . Sampling reversibilities. The rate constants I; are subject to thermodynamic constraints, Which can be eXploited in sampling by introducing a reversibility parameter (R1) for each elementary reaction step describing the ratio of the reverse and forward elementary fluxes .
Briefly, this principle states that for any reaction at equilibrium the frequency of transitions in both directions is the same for each individual reaction step. The same applies for non-sequen-tial reaction (random-order) mechanisms at equilibrium provided that the alternative paths are equivalent as far as kinetic order is concerned . In the reference state considered here, the enzymes are not at equilibrium and the probability ratio of the forward and reversed trajectories is not unity (minimum entropy condition), but rather is given by the eXponential of the total change in entropy along the forward process. We have extended the principle of microscopic reversibility to the non-equilibrium steady-steady flux (81 Text) producing the following criterion,
In Equation 10 the sum is made over all sets of reversibilities capable of carrying out the reaction. Such sets are called fundamental cycles and it has been demonstrated that a necessary and sufficient condition for the system satisfying the microscopic reversibility principle, is for each of these cycles to satisfy the criterion (Fig. 1B-C) . The fundamental cycles are found by traversing the edges of a graph connecting substrates to products and translating the paths into linear relations constraining different reversibility sets. A reversibility matriX (Qrev) is constructed that contains all the thermodynamic constraints in the reversibility vector ln(R). In order to sample the reversibilities, they are first normalised so that they can be sampled from the unitary simplex, where R2 represents the scaled reversibility. Using this transformation, Equation 11 can be cast into a more convenient form.
In order to sample ran-dom-order mechanisms, the Dirichlet distribution can be used in conjunction with Linear Programming (LP) techniques as shown in Equation 14. For each scaled reversibility, the lower bound is randomly sampled using the Dirichlet distribution and later Equation 13 is solved by maximizing the sum of the reversibilities. By fiXing the lower bound of each reversibility, we randomly generate points at different corners of the feasible hyperspace spanning all the solutions described by Equation 13. Using a sampled set of reversibilities satisfying Equation 11 and imposing the steady-state condition on each elementary reaction, 1'. e. v the elementary flux vector can be computed as the product of two elements, Where F is a diagonal matrix function of the reversibilities and relem is a branching vector taking into account the pattern stoichiometry (see below). The diagonal elements of F depend on the direction of the elementary reactions according to, Finally, the constant rate vector I} can be calculated by combining Equations 7 and 15 using the general equation.
For a complete derivation of this equation see the 81 Text. Calculation of the rate parameters as presented here ensures that they satisfy the fundamental thermodynamic principles under non-equilibrium conditions.
In order to determine the branching vector, relem, one needs to take into account the stoichiometry of the reaction pattern. This can be achieved by formulating and solving the mass balances for the forward and reverse elementary reactions (Equation 18).
1B). The solution space is infinite as cycles are per definition indeterminate. However, the sampling space can be constrained by formulating two additional restrictions: (1) elementary net fluxes are nonnegative and (2) the maximum elementary net flux in the pattern is equal to the reference flux for the R and T protomers (Equation 19).
Direct sampling v2me from this subspace can be difficult. An alternative approach is to eXpress all the paths solving Equation 18 as a linear combination of the null space basis of Spattern (Npattem), sample uniformly the weights (w) that yield a nonnegative solution and then normalize (Fig. 1B). This procedure is summarized in Equations 20 and computationally produces valid branching vectors more efficiently than direct sampling. In summary, for a monomeric enzyme we use the elementary reaction formalism and sample the microscopic rate constants indirectly by sampling the three elements in Equation 17: the reference state values of the enzyme intermediates, the reversibilities and the branching vector. This sampling strategy ensures that any parameter set sampled satisfies both pattern and kinetic constraints (Fig. lB-C).
Even if the reaction flux is known in a particular state, the particular values of both functions are unknown. To resolve this, we need to elucidate the contributions of the relaxed and tense conformations. The regulatory function is always less than or equal to 1, as the catalytic activity of the T state is less than the R state , and we introduce the activity ratio (aref) of a tense protomer at the reference state and sample this uniformly. Given are we can rearrange Equations 2 and 3 to calculate the contributions of each state The catalytic mechanism is assumed identical for the R and T state and we sample this mechanism tWice using the two different reference points to generate feasible parameterisa-tions for the R and T state. The final step is to generate a sample of parameters for the Q function that satisfy Equation 3 in the reference point. The Q function can be eXpressed as , Where L0 is the allosteric constant between the R and T states in the absence of ligands, xFJ-J-represent effector concentrations binding to specific allosteric sites, KRM and KTJ-J- denote the effectors dissociation constants for each state, éRO and em are the free enzyme fractions in both conformational states as function of the respective rate parameters (kR,kT) and reactant concentrations (x), 111 represents the number of allosteric sites, and finally, r and tare the number of positive and negative effectors binding to the allosteric sites in the R and T states, respectively. Here we have assumed that the allosteric activators and inhibitors only bind to the R and T, respectively , although this constraint can be relaxed. In the reference point, Q does not depend on the reactant and effector concentrations as they are defined as unitary at the reference point. We can determine all parameters in the Q function based on sampled enzyme abundances (refer to Equation 8).
These transitions point to a higher relative abundance of the free enzyme in the T state in the absence of substrate, Whereas in the presence of substrate the R state is more favoured . Ultimately, the latter yields L>1 and a ratio of affinity constants for both states Kfff/K;ff > 1 or equivalently KT/KR > 1 in the case of dissociation constants (Fig. 1D). Colosimo et a1.  have derived a simple expression to determine the allosteric constant in the absence of ligands assuming symmetric binding for the two states (Equation 24). To estimate the dissociation constants in Equation 24, we can make use of the equilibrium assumption between the R and T states and the unitary ligand concentration normalized at the reference state.
The same principle can be used to calculate the effectors dissociation constants. Calculation of each constant is achieved using the following general formula for the allosteric site 111, and effector r binding to the R state and effector tbinding to the T state. ref ref In particular, the absolute concentrations of the allosteric effectors (x x ) need not to be m,r7 m,t known, as they are scaling factors for the absolute effector concentrations in Equation 23 Which are unitary at the reference state. Parameter set accuracy check
Since this parameterization is built upon a reference point, this can be validated by confirming that the parameter set produces the reference flux at the reference point, i.e. When the normalized metabolite concentrations are unitary. We control the numerical accuracy of parameter sets by accepting only sets achieving the reference flux Within a tolerance, e. g. 8 2 10'8, In general, rejected instances are insignificant and normally represent less than 0.01% of the sampled models.
Microscopic rate constants are commonly difficult to measure; therefore standard macroscopic constants are preferred to parameterize rate expressions. Transformation of rate constants into macroscopic rate constants can be performed following Cleland’s rules , Which are consistent With the Haldane relationships. In this manner, estimated macroscopic constants can readily be compared With available eXperimental data.
In order to determine the impact of reactant perturbations on the reaction rate, we estimated the reaction elasticities upon an infinitesimal variation in the concentration of substrates and products . The partial derivatives were calculated using a central difference approximation, Where 5cref represents the perturbed normalized reactant concentration, 1: denotes the rate constants vector and Ah is the size of the perturbation. Given that the perturbations for reactant concentrations are performed in the vicinity of the reference state, i. e. Scref = 1, a uniform step size of 10'2 equivalent to a 1% change was employed for all calculations.
Definition of the reversibility matrix was achieved using appropriate MATLAB functions from the Bioinformatics toolbox. Automated derivation of rate laws was achieved using King-Altman’s method for finding valid reaction patterns . Qi et al.  have recently reported an efficient algorithm (KAPattem) employing topological theory of linear graphs for accomplishing this goal. By employing this algorithm, we derived the enzyme intermediates abundance functions which we then assembled to build the final velocity rate. All computations were run on a Dell OptiPlex 990 Desktop (Intel Core i5-2400, 4 GB ram, Microsoft Windows 7, x86-based architecture).
1). In the following, we Will describe several applications of this platform to assess the impact of thermodynamic affinity, reaction molecularity and mechanisms, as well as modelling compleX kinetic behaviours.
Depending on the mechanism of reaction, these relations relate the values of the macroscopic kinetic constants, i.e. dissociation and catalytic constants, with the apparent equilibrium constant of the reaction. It would be interesting to derive similar relations under non-equilibrium steady-state conditions, provided that biological systems operate in this regime. It has been demonstrated that the analysis of the relation between the thermodynamic affinity (represented by AGr/RT) and the reaction velocity can still be performed as if it were at equilibrium, by displacing the reference point from the equilibrium to another appropriate reference state . As such, relations analogous to the Haldane relationships can be derived at this reference point (81 Text). For example, for a Uni-Uni reaction converting a substrate A into a product P the following relation can be derived, Where I; denote the rate constants at the reference state, I~<A and I~<P represent the normalized and I; the scaled catalytic constants for the forward and reverse reactions. It can be demonstrated that dissociation constants for A and P at the reference state, respectively, and I; are the following generalized equation holds true for any reaction mechanism, where s and p denote the number of substrates and products, respectively. We have designated the first term on the right-hand side of Equation 30 the catalytic (turnover) term, while the second is called the binding (saturation) term. Notably, the catalytic term is independent of the reactant concentrations at the reference state, as opposed to the saturation term. In fact, it can be shown that I~<A = A/Aref, which can be regarded as the degree of saturation of the enzyme for reactant A. In this way, Equation 30 enables the energetic analysis of any reaction by decomposing it in two contributions: (1) catalysis (how efficient is the enzyme converting substrates into products at the reference state) and (2) binding (how saturated is the enzyme at the reference state). More importantly, this relation imposes a natural tradeoff in enzyme catalysis. Greater contributions of the catalytic term suggest a higher enzymatic efficiency in the conversion of substrates to products, but a suboptimal saturation of the enzyme, 1'. e. K A/Aref high. On the contrary, high saturation contributions suggest lower conversion efficiency of substrates into products, i.e. kcak / 12w? + is relatively higher compared to the binding term.
SI Table shows the definition of the dissociation and catalytic constants in terms of the rate constants using Cleland’s definitions. We performed the analysis considering conditions close to equilibrium (-1 kI/mol) to far from equilibrium (-80 kI/mol) (Fig. 2). Notably, the average conversion and saturation contributions ratio remains constant for different AGr/RT for all the kinetics sampled. As eXpected, the higher the molecularity of the reaction, the higher the contribution of the binding term. The Uni-Uni mechanism eXhibits a low average binding contribution (32%), which increases for the BiBi mechanism (60%) and further for the Ter-Ter mechanism (71%). These results suggest that catalysis in multi-substrate enzymes is a process strongly driven by the degree of saturation of the enzyme and to a lesser extent by the actual conversion of substrates into products.
For the aforementioned cases, the sums of the blue and red areas represent the 95% confidence region of all the thermodynamically feasible parameter sets. Interestingly, for all the tested mechanism, the feasible area increases with higher driving force (AGr/RT more negative). Such increased area point to a greater diversity of feasible parameter sets under more thermodynamically favourable conditions. On the contrary, more homogeneous parameter sets can be found closer to equilibrium, i.e. different parameter sets have similar energetic contributions. The latter suggests that sampled parameter sets might have similar kinetic behaviours closer to equilibrium. The next analysis will provide additional support for this assertion.
We next analysed the impact of thermodynamics on enzyme kinetics. To that end, we uniformly sampled the kinetic space of reactions with same molecularity. Since the majority of enzymes found in nature catalyse bimolecular reactions  , we focussed our attention on the most representative bimolecular mechanisms. There are two general mechanisms for bimolecular reactions, the ping-pong mechanism in which a product is released before both substrates have reacted with the enzyme, and the ternary complex (sequential) mechanism in which the enzyme combines with both substrates before products are formed . Sequential mechanisms can be further divided into random and ordered mechanisms. In random-order mechanisms, reactants can bind in either order, while in the ordered type sequential process one reactant always binds to a certain site before a second reactant binds to the other site . The representation of all three mechanisms is shown in Fig. 3A. The impact of thermodynamics was evaluated in each case by calculating the reaction sensitivities for substrates and products, i.e. substrate and product elasticities, at different Gibbs free energy differences ranging from -1 (close to equilibrium) to -80 kI/mol (practically irreversible). The 95% confidence regions for the calculated reaction elasticities for the above mentioned reaction mechanisms are shown in Fig. 3B-C.
Within the range of AG, analysed, two reaction elasticity regions can be distinguished: a variable elasticity region (O>AGr>-20 kI/mol) and a constant elasticity region (AGr<-20 kI/mol). The former region can be further subdivided in two: a linear regime with a steep slope (O>AGr >-2 kI/mol) and a transition regime (-2>AGr>-20 kI/mol). Notably, previous works have demonstrated an almost perfect linear correspondence between the thermodynamic affinity and the reaction flux close to equilibrium, i.e. O>AGr > -1.5 kI/mol [38,39]. Such relationship has been shown to hold for AGr>-7 kI/mol with less than 15% error . Our simulation results are thus in line with what would be eXpected to be the kinetic behaviour close to equilibrium.
In the case of the reaction elasticities close to equilibrium, our sampling results show the substrate and product elasticities are respectively monotonically decreasing and increasing, reaching their respective maximum and minimum at equilibrium (Fig. 3B). This behaviour is consistent With previous analyses on the behaviour of the reaction elasticities in this region  and supports our intuition: reactions operating close to equilibrium are more susceptible to slight changes in the reactant concentrations, eXhibiting greater changes upon these perturbations. On the other hand, our analysis far from equilibrium revealed both substrate and product elasticities reach a plateau at highly negative thermodynamic affinities for all the reaction mechanisms analysed. In the case of the substrate elasticity, this plateau stabilizes close to zero for highly negative AGr. Notably, the average substrate elasticity for thermodynamically favoured reactions is approximately 0.22 and not 0 as it would be eXpected. The latter is a direct result of our sampling strategy, Which seeks to uniformly sample all the possible parameter sets that are consistent With the Haldane relationships at a chosen reference point. On the contrary, in the case of the product elasticity, the average elasticity consistently approaches zero for all the mechanisms analysed under favoured conditions. Furthermore, our sampling results suggest that product inhibition is almost negligible for thermodynamically favoured reactions, i.e. AGr< -30 kI/mol (Fig. 3C). 81 Text provides an illustrative explanation of the asymptotic behaviour of the reaction elasticities near and far from equilibrium. Additionally, it is also important to highlight that the allowable elasticity regions becomes tighter as we move closer to equilibrium (Fig. 3B-C). The latter supports our earlier findings of sampling more homogenous parameter sets close to equilibrium. Indeed, the closer to equilibrium the tighter the feasible parameter space becomes, and thus, a more similar response upon reactant perturbations is obtained.
However, the average response of ordered and random mechanisms is more similar than the response of the ping-pong mechanism. This is most readily appreciated when comparing the substrate elasticities (Fig. 3B). As previously mentioned, the ping-pong mechanism executes a fundamentally different mechanism in which the release of the first product takes place before binding of the second substrate. The latter explains the slightly lower average substrate elasticity. Interestingly, negative substrate elasticities can be encountered in the random order mechanism for high thermodynamic affinities (approx. AGr<- 18 kI/mol), i.e. the addition of substrate decreases the velocity of reaction (elasticities below the red dashed line in Fig. 3B). This behaviour has been demonstrated previously for random-order mechanisms . Theoretical analysis of these mechanisms has shown that they can give rise to apparent substrate inhibition or substrate activation . For example, depending on the kinetic parameter values and the substrate concentrations, the reaction rate for the enzymes following this mechanism can display an apparent cooperative behaviour (sigmoidal reaction rate) upon addition of one substrate maintaining the other constant, while in the opposite situation they can exhibit substrate inhibition (the reaction rate pass through a maximum) . Although this behaviour is not very common, there is evidence of such in the literature . More importantly, our framework enabled unbiased sampling of all feasible kinetic behaviours, thereby revealing the impact of the thermodynamic affinity on the reaction rate.
In mammals, this step is performed by four different isozymes (EC 188.8.131.52) located in different tissues. Hexokinases I-III are mainly located in the brain, muscle and erythrocytes , While hexoki-nase IV (glucokinase) is primarily located in the liver and pancreatic B-cells . In the latter tissues, glucose phosphorylation is the rate-limiting step of glucose metabolism, and more importantly, it is ultimately responsible for the release of insulin into the bloodstream Which maintains glucose homeostasis in the body.
It displays a sigmoidal response upon increasing glucose concentration, i.e. positive cooperativity, but has a hyperbolic behaviour upon increasing MgATPZ', i.e. Michaelis-Menten kinetics . This positive cooperativity for glucose is intriguing as the enzyme is monomeric, Which contradicts standard cooperativity models and thus requires a conceptually different eXplanation. One of the simplest models capable of explaining this behaviour is the mnemonic model . Briefly, this model proposes that the conformation of an enzyme following product release could be different from the initial enzyme state, i.e. the enzyme has memory. In the case of glucokinase, this model suggests a slow conformational transition from a low-affinity state (E*) to a high-affinity state (E), Which ultimately carries out both catalysis and product release (Fig. 4). This transition can be enhanced With increasing glucose concentrations, yielding the observed positive cooperative behaviour on this substrate (kinetic cooperativity). Other possible explanations like a reaction mechanism With random addition of substrates has been
Schematic representation of the mnemonic model for the mammalian glucokinase. This model proposes a slow conformational transition from a low-affinity state (E*) to a high-affinity state (E) that can be enhanced with increasing glucose concentration yielding the observed cooperative behaviour. discarded, as isotope exchange studies have demonstrated an ordered mechanism with glucose binding first .
In this case, we are interested in estimating the frequency of positive cooperativity for glucose of the glucokinase. To this end, we uniformly sampled the kinetic space for this enzyme and counted the frequency of parameter sets displaying positive cooperativity for glucose (estimated Hill coefficient, nH>1). In order to mimic enzymatic assays conditions, i.e. high substrate concentrations and almost no products, a reference Gibbs free energy difference of - 100 kI/mol was used. Very similar results were found using down to -50 kI/mol of AG.. The latter was eXpected as reaction elasticities were found to reach a plateau for AGr<-30 kI/mol (Fig. 3B-C). Nevertheless, in order to ensure initial velocity conditions, we chose the former difference. The reference steady-state flux under this condition was set to the experimental value of 0.064 mM/min found in the literature .
5A). However, a small portion of the models (~7%) displayed an apparent negative cooperative behaviour. The latter is possible due to the existence of competing parallel pathways in the reaction mechanism, i.e. E —> E-glc —> E-glc-atp —> E-g6p-adp —> E-adp —> E and E —> E* —> E-glc —> E-glc-atp —> E-g6p-adp —> E-adp —> E (Fig. 4). As previously mentioned, reaction mechanisms with branched steps can eXhibit apparent substrate inhibition depending on the kinetic parameter values and the substrate concentrations. In this case, this will depend on the value of the branching factor for the E —> E-glc and E —> E* steps. In particular, if a large elementary net flux is sampled for the E —> E* step at the reference state, then the successive addition of glucose will inhibit the velocity of the reaction as opposed to enhancing it.
5B). This kinetic behaviour has been extensively studied and there is abundant supporting evidence . Remarkably, the sampled kinetics contained the experimentally observed Hill coefficient for this enzyme (nHJeal = 1.70 :01 ), within one standard deviation (nH,sampled = 1.77:0.5), which confirms the suitability of the mnemonic model for modelling this kinetic behaviour. Fig. 5C shows a comparison of the fitted kinetic model by Storer and Cornish-Bowden  and the sampled kinetics using our framework. Better agreement between the models is encountered at high glucose concentrations (>20 mM), i.e. where the cooperative behaviour is less pronounced, while at lower glucose concentrations (<4 mM) slightly higher discrepancies can be observed. As eXpected, a perfect match between both models is found at the reference state, as by construction the proposed framework builds the parameterization at this point (vref = 0.064 mM/min). Interestingly, the variability of the sampled kinetics is not proportional to the glucose concentration as opposed to fitted model (Fig. 5C). For some glucose concentration regions, the kinetic behaviour of the sampled models displays a relative greater variability for low glucose concentrations (<20 mM) and a relative lower variability for higher concentrations. The latter shows that thermodynamically consistent parameter sampling can provide means for effectively accounting for the feasible kinetic space without the need for excessive data. Moreover, the addition of extra data will further reduce the feasible kinetic space for this particular mechanism. Taken together, these results show that by using this framework, eXploration of the kinetic space for compleX reactions can be achieved provided that the reaction mechanism is known and by setting the right thermodynamic constraints.
In glucose-limited Escherichia coli cultures, PEPC is involved in the only active anaplerotic reaction replenishing pools of intermediary metabolites of the tricarboxylic acid (TCA) cycle  (Fig. 6A). This enzyme catalyses the conversion of phos-phoenolpyruvate (PEP) and bicarbonate into oxaloacetate and orthophosphate in the presence of Mg2+. This reaction is highly exergonic (AG? 2 -43.2 kI/mol), making it practically irreversible under physiological conditions. The catalytic mechanism of this enzyme has been elucidated With an ordered binding of phosphoenolpyruvate first and releasing of oxaloacetate as the final step (Fig. 6B) . The regulatory mechanism behind the operation of this enzyme is much more complex. The enzyme is allosterically activated by acetyl-CoA, long-chain fatty acids and their acyl coenzyme A derivatives, fructose 1,6-bisphosphate (FBP) and the nucleo-tides guanosine-5’-triphosphate and cytidine-5’-diphosphate, whereas it is inhibited by L-aspartate and L-malate [53—56] (Fig. 6A). In particular, the mechanism of activation of this enzyme upon the combined action of FBP and acetyl-CoA has been studied in detail. These two activators bind to different allosteric sites, synergistically activating the tetrameric active form . Alone FBP has been shown to exert little activation of the enzyme upon binding of PEP; however acetyl-CoA alone can greatly enhance the activation induced by FBP . A plausible synergistic model capable of explaining this behaviour has been proposed by Smith et al.  and it is shown in Fig. 6C. Notably, this model is a hybrid between that of Monod et al.  and that of Koshland et al.  in that the mechanism of transitions requires the presence of the activators and/ or ligand (PEP) for the enzyme to be active (induced fit), however it also assumes isomerization of all subunits to describe cooperative interactions (concerted model). The particular kinetic features of the PEPC have been shown lately to play a key role in the regulation of anapleurosis in E. coli. Xu et al.  have demonstrated that E. coli is able to turn off PEP consumption quickly upon glucose removal thanks to the ultrasensitive response of the PEPC upon FBP depletion. Following glucose removal, depletion of FBP from 15 mM to 0.45 mM has been shown to almost entirely turn off PEP consumption . This rapid response enabled accumulation of PEP for the rapid import of glucose when it becomes available again. Using the present framework, we sampled the complex regulatory behaviour using the mechanistic information depicted in Fig. 6B-C. For the regulatory mechanism, however, we also considered the tense form to be capable of performing the reaction, although with a lower activity compared to the relaxed form. Indeed, the model of Smith et al.  can be regarded as a special case of our parameterization with aref
In this case, we are interested in describing the ultrasensitive response of PEPC in the presence/absence of acetyl-CoA under changing FBP concentrations. To this end, experimental data from Xu et al.  was used as reference data to build and sample feasible kinetics. The thermodynamic reference state was chosen to ensure initial velocity as done for the mammalian glucokinase. To assess the performance of our framework, we compared our results against an empirical model developed by Lee et al.  and calibrated using data collected under similar conditions . The empirical model was adjusted to the same reference state used during sampling to ensure a fair comparison. Our sampling strategy accurately described the kinetic behaviour of the PEPC for different FBP concentrations in the presence of acetyl-CoA at physiological concentrations (Fig. 7A). Moreover, it eXhibited a slightly better performance than the model of Lee et al.  under the same condition. However, a worse performance of our approach is observed in the absence of acetyl-CoA (Fig. 7A). This was eXpected as our framework builds kinetic models around one reference state, which in this case was set to 0.63 mM acetyl-CoA and 2 mM PEP. When ace-tyl-CoA is absent, a larger diversity of plausible kinetics is predicted by our sampling approach, although displaying a more sigmoidal behaviour. In order to improve the fitting to this condition, we can perform a rejection step during the sampling so that every accepted parameter set agrees with the experimental data under this condition. This strategy is typically used in Bayesian Inference by Approximate Bayesian Computation (ABC) methods . In particular, the above method corresponds to the simplest ABC method, the rejection sampler . There are other more efficient samplers implemented within the ABC setting to compute the parameters’ distributions [62,63]; however we opted for the simplest as a-proof-of-principle to demonstrate our strategy. As it can be observed in Fig. 7B, the inclusion of additional experimental data further constraints the plausible kinetic space. Interestingly, even with the addition of extra information during the sampling, the most accurate description for this kinetics displays a sigmoidal behaviour as opposed to the observed hyperbolic kinetics. Indeed, PEPC activation is a fairly complex phenomenon and involves the synergistic interplay of two classes of effectors (type I, e.g. FBP, and type II, e.g. acetyl-CoA) . In particular, PEPC activation can be achieved by the sole action of FBP or combined with acetyl-CoA. This behaviour has been previously attributed to play a key role in the rapid adaptation of E. coli from normal-growing culture conditions to carbon-starvation or acetate switch conditions .
7C, see SI Text for generation of Hill curves). In order to compare our results, we also show experimentally measured Hill coefficients reported by Izui et al.  obtained under similar conditions. Firstly, our results are consistent with the theory as they show bell-shaped Hill curves reaching an asymptotic value of unity at either PEP —> 0 or PEP —> oo . Furthermore, they suggest acetyl-CoA is the most powerful activator. As the curves move to the left, the apparent affinity constants increase. The largest displacements are observed in the presence of acetyl-CoA, which supports its role as an invariant activator for the synergistic co-activation of PEPC . Another consequence of the latter is the increase of the enzyme forms in the relaxed state. As the apparent affinity increases, more enzyme forms transition from the tense to the relaxed state. In terms of the prediction of the cooperativity under the different conditions studied, the sampled kinetics is in close agreement with the experimental data. Notably, a theoretical higher cooperativity for PEPC upon PEP binding is observed in the presence of FBP rather than acetyl-CoA (max(nH,AcetY1_C0A(_),FBp(+)) = 2.05, max (11H, Acety1_Co A(+)FBP(_)) = 1.4) (Fig. 7C). These results agree with previous reports indicating low cooperativity for acetyl-CoA alone (1.07< nH < 1.2 for 0.02 < acetyl-CoA < 1.0 mM and 0.1 < PEP < 50 mM [55,57]), whereas significant cooperativity for the FBP interaction [53,55,65]. Moreover, our results predict a maximal Hill coefficient of approx. 1.15 for PEPC in the absence of both FBP and acetyl-CoA which is close to the reported value of 1.2 reported by Izui et al. .
Under carbon-starvation conditions (low FBP) and at acetyl-CoA physiological concentrations, PEPC was only activated in ~10% (Fig. 7D). The presence of both effectors at physiological concentrations as seen in glu-cose-limited cultures enhanced PEPC activity to almost 99%. Altogether, our results predict a drop of ~90% in the activity of PEPC upon shifting from normal culture conditions to carbon starvation. Experimental values determined for this enzyme simulating these conditions yielded an approx. 94% decrease in its activity  , close to the one predicted by our model. More importantly, the employed parameterization provides the means to understand how this enzyme is being regulated. As such, this framework is not only capable of sampling and modelling complex kinetic behaviours, but it also provides useful insights into the regulatory mechanisms underpinning those kinetics.
Parameterization relies on the integration of the generalized MWC model for modelling the kinetics of oligomeric enzymes with the elementary reaction formalism for deriving the catalytic rates and thermodynamic constraints between rate parameters. As a result, this framework enables exploration of the feasible kinetic space of models following a particular reaction mechanism, displaying a given reference flux under a specific thermodynamic condition. Exploration of the kinetic space in this way provided the necessary tools for evaluating the impact of thermodynamics on enzyme kinetics, as well as the consequences of particular regulatory features on the kinetic behaviour.
To that end, a simple relation between conversion terms (catalytic constants) and saturation terms (dissociation constants) was derived. As expected, higher molecularities increase the relative importance of the saturation term on average compared to the conversion term. The latter suggest that higher impacts on the catalysis of multi-substrate enzymes would be expected by modifying its ligand binding and releasing properties rather than the conversion rate of substrates into products. In fact, most common approaches for multi-substrate enzyme engineering involve modifying substrate specificity and selectivity , although the success of these strategies will ultimately depend on the particular reaction mechanism and enzyme. Notably, our approach could be employed to sample the rate constants of a particular reaction and determine the flux control coefficient distributions of each step in the mechanism . The latter would provide a broader picture of the pattern regulation, enabling the identification of rate-limiting steps whose modification would most likely improve desired kinetic properties.
It has been shown that the reaction elasticities strongly depend on this variable. Our analysis distinguished three elasticity regions as a function of the thermodynamic affinity independent of the bimolecular mechanism analysed: a linear sensitivity region with a steep slope (O> AGr>-2 kI/mol), a transition region (-2> AGr>-20 kI/mol) and a constant sensitivity region (AGr< -20 kI/mol). Notably, substrate and product elasticities reached their highest absolute values close to equilibrium, which agrees with the tendency of substrate elasticities to approach infinity at equilibrium . This framework also enabled analysis far from equilibrium. In all reaction mechanism tested, both substrate and product elasticities reached a plateau, reflecting the saturation state of the enzyme. In terms of the reaction mechanisms, all of them displayed elasticities across the thermodynamic reference states, although a particular behaviour could be distinguished for the random mechanism. The latter suggests that knowledge of the thermodynamic state can be more valuable than exact determination of the reaction mechanism at least for nonrandom reaction mechanisms.
To illustrate this, we sampled the kinetic behaviours of the mammalian glucokinase and the phosphoenolpyruvate carboxylase of E. coli. In the case of the mammalian glucokinase, we were able to model its positive cooperativity for glucose by employing the mechanistic mnemonic model. Random sampling of the kinetic behaviour for this enzyme showed that the experimental Hill coefficient for this enzyme lies surprisingly fairly close to the average Hill coefficient sampled, highlighting the importance of the architecture of its reaction mechanism.
In particular, application of our framework to other monomeric proteins exhibiting allosteric behaviours, e.g. thrombin, myoglobin [68,69] , might be a valuable tool to better understand the mechanisms underpinning their kinetics. In the case of PEPC, eXploration of the kinetic space provided valuable information related to its regulation. Not only were we able to describe the compleX co-activation behaviour of this enzyme by FBP and ace-tyl-CoA, but also to determine the overall impact of this co-activation on its cooperativity for PEP. Furthermore, our strategy offered insights into the ultrasensitive regulation of PEPC upon binding of FBP and acetyl-CoA. As in the case of mammalian glucokinase, our results emphasises the importance of using mechanistic/phenomenological models for describing and interpreting kinetic behaviours.
The sampling of enzymatic reactions within metabolic pathways seems particularly promising. The difficulty of parameterizing metabolic pathways is widely recognised with the main difficulty being the large number of parameters and the relative few data available. Rather than grossly simplifying kinetics to enable fitting, this work suggests that it is possible to build feasible, accurate kinetic parameterizations with limited data by integrating phenomenological models with efficient sampling techniques. In particular, we believe that eXperimentalists could greatly benefit from this framework in those cases where the fitting is difficult or requires large amounts of data. It will be therefore our next step to deploy GRASP as a software package to assist in such complicated tasks.
81 Table. Definition of kinetic constants in terms of rate constants for the different mecha-81 Text. Mathematical derivations of the formalism, kinetic and thermodynamic expressions.
Performed the experiments: PS. Analyzed the data: PS. Contributed reagents/materials/ analysis tools: PS. Wrote the paper: PS LKN. Rubinstein RY (1982) Generating Random Vectors Uniformly Distributed inside and on the Surface of Different Regions. Eur J Oper Res 10: 205—209.
The Journal of biological chemistry 255: 1635—1642. PMID: 6986377 Del Moral P, Doucet A, Jasra A (2011) An adaptive sequential Monte Carlo method for approximate Bayesian computation. Statistics and Computing 22: 1009—1020. Wyman J, Gill SJ (1990) Binding and Linkage Functional Chemistry of Biological Macromolecules. Mill Valley, California: University Science Books.
See all papers in April 2015 that mention rate constants.
See all papers in PLOS Comp. Biol. that mention rate constants.
Back to top.
See all papers in April 2015 that mention parameter sets.
See all papers in PLOS Comp. Biol. that mention parameter sets.
Back to top.
See all papers in April 2015 that mention substrates.
See all papers in PLOS Comp. Biol. that mention substrates.
Back to top.
See all papers in April 2015 that mention kinetic parameters.
See all papers in PLOS Comp. Biol. that mention kinetic parameters.
Back to top.
See all papers in April 2015 that mention experimental data.
See all papers in PLOS Comp. Biol. that mention experimental data.
Back to top.
See all papers in April 2015 that mention free energy.
See all papers in PLOS Comp. Biol. that mention free energy.
Back to top.
See all papers in April 2015 that mention Hill coefficient.
See all papers in PLOS Comp. Biol. that mention Hill coefficient.
Back to top.
See all papers in April 2015 that mention E. coli.
See all papers in PLOS Comp. Biol. that mention E. coli.
Back to top.
See all papers in April 2015 that mention energy difference.
See all papers in PLOS Comp. Biol. that mention energy difference.
Back to top.
See all papers in April 2015 that mention steady-state.
See all papers in PLOS Comp. Biol. that mention steady-state.
Back to top.
See all papers in April 2015 that mention ligands.
See all papers in PLOS Comp. Biol. that mention ligands.
Back to top.