📒 Saa 2015 | A General Framework for Thermodynamically Consistent Parameterization and Efficient Sampling of Enzymatic Reactions

A General Framework for Thermodynamically Consistent Parameterization and Efficient Sampling of Enzymatic Reactions1



  • Michaelis and Menten enzyme kinetics: quasi-steady-state assumption for intermediates
    • can be automated using King-Altman’s (KA) method
  • Related to the apparent equilibrium constant of the overall reaction by the Haldane relationships
  • 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. Computational sampling enables us to determine emergent enzyme properties from high dimensional parameter spaces

Allosteric regulation

  • Current parameterization and sampling approaches also fail to accurately capture allosteric regulation. Metabolic control is mostly achieved through allosteric and transcriptional regulatory interactions
  • symmetry model of Monod, Wyman and Changeux (MWC)
  • sequential model of Koshland, Némethy and Filmer (KNF)
  • assumed equilibrium between two conformational states: relaxed (R) and tense (T) state
  • the MWC model is regarded the model of choice for describing allosteric and cooperative interactions


General Reaction Assembly and Sampling Platform (GRASP)

  • parameterizing and sampling the kinetics of any oligomeric enzyme by using minimal reference and biochemical mechanistic data
  • the parameterization retains all intrinsic thermodynamic constraints between kinetic parameters, obeying the principle of microscopic reversibility


General framework for modelling metabolic reactions

  • employs the MWC model for an oligomeric enzyme (one active site), developed by Popova and Sel’kov
  • v is the product of two independent functions (catalytic: R; regulatory: R<-> T):
  • Q is a function that determines the current ratio between the R and T states

$$ \begin{aligned} v &= \Phi_{catalytic} \cdot \Psi_{regulatory} \cr \Phi_{catalytic} &= n \cdot v_{R} \cr \Psi_{regulatory} &= \frac{1 + (v_T / v_R) Q}{1 + Q} \end{aligned} $$

  • decomposing the reaction velocity into two independent functions
  • enables inclusion of fundamental thermodynamic relations between kinetic parameters
  • compatible with the elementary reaction formalism

Parameterization and sampling of the catalytic mechanism

  • Every enzymatic reaction can be broken into simple reversible steps called elementary reactions, by law of mass action: $$ \begin{array}{l}{v_{\text { elem }}=k \cdot x \cdot e \text { for binding steps }} \cr {v_{\text { elem }}=k \cdot e \text { for dissociation steps }}\end{array} $$
  • The absolute values of the metabolite and total enzyme concentrations are not known (although physiological ranges can be estimated), normalization of all the variables around a reference point (steady-state flux) is a convenient strategy. The normalized metabolite concentrations are unitary at the reference point (x˜ref=1)
    • For binding steps $$ v_{\text { elem }}=\left(k \cdot e_{\text { totel }}^{\text { ref }} \cdot x^{\text { ref }}\right) \cdot\left(\frac{x}{x^{\text { ref }}}\right) \cdot\left(\frac{e}{e_{\text { total }}^{\text { ref }}}\right)=\tilde{k} \cdot \tilde{x} \cdot \tilde{e} $$
    • For dissociation steps $$ v_{\mathrm{elem}}=\left(k \cdot e_{\mathrm{total}}^{\mathrm{ref}}\right) \cdot\left(\frac{e}{e_{\mathrm{total}}^{\mathrm{ref}}}\right)=\tilde{k} \cdot \tilde{e} $$
  • Rate laws can be derived from the enzyme mechanism and the microscopic rate constants using the King-Altman’s method.
  • 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

$$ v_{\text { elem }}^{\text { ref }}=P\left(\tilde{e}^{\text { ref }}\right) \cdot \tilde{k} $$

  • v = elementary fluxes, P = diagonal matrix function, k = vector of rate constants
  • The enzyme intermediate abundances sums to one and can be readily sampled using appropriate probability distributions. Equivalent as sampling from a multivariate Dirichlet distribution with hyper-parameter vector 1, using Gamma distributions. $$ \sum_{i=1}^{p} \tilde{e}_{i}=1 $$

Sampling reversibilities

  • The rate constants k˜ are subject to thermodynamic constraints with reversibility parameter (R_i) for each elementary reaction step $$ R_{i}=\left(\frac{v_{\text { elem }, i}^{\text { reverse }}}{v_{\text { elem }, i}^{\text { forward }}}\right)^{\operatorname{sgn}\left(v^{\text { ref }}\right)} $$

  • for any reaction at equilibrium the frequency of transitions in both directions is the same for each individual reaction step

  • extended the principle of microscopic reversibility to the non-equilibrium steady-steady flux (S1 Text) producing the following criterion: $$ \sum_{i \in \text { fundmental cycle }} ln(R_i)=\operatorname{sgn}\left(v^{\mathrm{ref}}\right) \cdot \frac{\Delta G_{r}}{R T} $$

  • A reversibility matrix (Ωrev) is constructed that contains all the thermodynamic constraints in the reversibility vector ln(R) $$ \Omega_{\mathrm{rev}} \cdot \ln (\boldsymbol{R})=\operatorname{sgn}\left(v^{\mathrm{ref}}\right) \cdot \frac{\Delta G_{r}}{R T} $$

  • normalised form: $$ \ln \left(\widehat{R}_{i}\right)=\frac{\ln \left(R_{i}\right)}{\operatorname{sgn}\left(v^{\mathrm{ref}}\right) \cdot \Delta G_{r} / R T} $$ $$ \begin{array}{l}{\Omega_{\mathrm{rev}} \cdot \ln (\widehat{\boldsymbol{R}})=\mathbf{1}} \cr {0 \leq \ln \left(\widehat{R}_{i}\right) \leq 1}\end{array} $$

  • For compulsory-order mechanisms, Ωrev = 1

  • In order to sample random-order mechanisms, the Dirichlet distribution can be used in conjunction with Linear Programming (LP) techniques. Lower bound (lb) is randomly sampled using the Dirichlet distribution.

  • the elementary flux vector can be computed as the product of two elements $$ v_{\mathrm{elem}}^{\mathrm{ref}}=\Gamma \cdot r_{\mathrm{elem}} $$

    • The diagonal elements of Γ depend on the direction of the elementary reactions according to $$ \Gamma_{i, i}^{\text { forward }} :=\frac{1}{1-R_{i}^{s g n\left(v^{\mathrm{ref}}\right)}} $$ $$ \Gamma_{i, i}^{\mathrm{reverse}} :=\frac{R_{i}^{s g n\left(v^{\mathrm{ref}}\right)}}{1-R_{i}^{s g n\left(v^{\mathrm{ref}}\right)}} $$
  • the constant rate vector k˜ can be calculated by combining Equations 7 and 15 using the general equation, satisfing the fundamental thermodynamic principles under non-equilibrium conditions. $$ \tilde{k}=P^{-1}\left(\tilde{e}^{\mathrm{ref}}\right) \cdot \Gamma(R) \cdot r_{\mathrm{elem}} $$

Sampling elementary flux vectors

  • solving the mass balances for the forward and reverse elementary reactions $$ S_{\text { pattern }} \cdot\left(v_{\text { elem }}^{\text { forward }}-v_{\text { elem }}^{\text { reverse }}\right)=S_{\text { pattern }} \cdot v_{\text { elem }}^{\text { net }}=0 $$
  • two additional restrictions: (1) elementary net fluxes are non-negative and (2) the maximum elementary net flux in the pattern is equal to the reference flux for the R and T protomers $$ \begin{array}{l}{v_{\text { elem }}^{\text { net }} \geq 0} \cr {\max \left(v_{\text { elem }}^{\text { net }}\right)=v^{\text { ref }}}\end{array} $$
  • expressing all the paths solving Equation 18 as a linear combination of the null space basis of S pattern (N pattern) $$ \begin{aligned} \mathbf{v}_{\text { elem }}^{\text { net }} &=\mathbf{N}_{\text { pattern }} \cdot \mathbf{w} \cr \mathbf{r}_{\text { elem }} &=\boldsymbol{v}^{\text { ref }} \cdot \frac{\mathbf{v}_{\text { elem }}^{\text { net }}}{\max \left(\mathbf{v}_{\text { elem }}^{\text { net }}\right)} \end{aligned} $$

Sampling functional contributions: catalytic and regulatory effects

  • the activity ratio (a ref) of a tense protomer at the reference state and sample this uniformly. $$ a^{\mathrm{ref}}=\frac{v_{\mathrm{T}}^{\mathrm{ref}}}{v_{\mathrm{R}}^{\mathrm{ref}}} \sim \operatorname{Uniform}(0,1) $$

  • calculate the contributions of each state $$ \begin{array}{l}{v_{\mathrm{R}}^{\mathrm{ref}}=\frac{\Phi_{\mathrm{catalytic}}^{\mathrm{ref}}}{n}=\frac{v^{\mathrm{ref}}}{n \Psi_{\mathrm{refulatory}}^{\mathrm{ref}}}} \cr {v_{\mathrm{T}}^{\mathrm{ref}}=a^{\mathrm{ref}} \cdot v_{\mathrm{R}}^{\mathrm{ref}}}\end{array} $$

  • The Q function

    • L_0 is the allosteric constant between the R and T states in the absence of ligands
    • x_F,i,j represent effector concentrations binding to specific allosteric sites
    • K_R,i,j and K T,i,j denote the effectors dissociation constants for each state
    • e˜R0 and e˜T0 are the free enzyme fractions in both conformational states as function of the respective rate parameters
    • m represents the number of allosteric sites
    • r and t are the number of positive and negative effectors binding to the allosteric sites in the R and T states
  • 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

  • L_0 is esimated as follows: $$ L_{0}=\left(\frac{\tilde{e}_{\mathrm{T} 0}^{\mathrm{ref}} / \tilde{e}_{\mathrm{T} 1}^{\mathrm{ref}}}{\tilde{e}_{\mathrm{R} 0}^{\mathrm{ref}} / \tilde{e}_{\mathrm{R} 1}^{\mathrm{ref}}}\right)^{n / 2} $$ (e˜ref0) and (e˜ref1) denote the enzyme fractions free and bound to the ligand for the R and T states at the reference state

  • Calculation of dissociation constants

Parameter set accuracy check

  • confirming that the parameter set produces the reference flux at the reference point, within a tolerance ϵ $$ \left| \boldsymbol{\Phi}_{\text { catalytic }}^{\text { ref }} \cdot \Psi_{\text { regulatory }}^{\text { ref }}-v^{\text { ref }}\right|<\varepsilon $$

  • If the reference values for the total enzyme and metabolite concentrations are known, they can be used to transform the set of scaled rate constant into absolute constants. Standard macroscopic constants are preferred to parameterize rate expressions by Cleland’s rules, which are consistent with the Haldane relationships.

Elasticity analysis of the velocity of reaction

  • impact of reactant perturbations on the reaction rate, using a central difference approximation

Computational implementation

  • MATLAB 2013a, Bioinformatics toolbox + KAPattern


Dissecting enzyme catalysis: assessing the impact of reaction molecularity

  • The connection between reaction thermodynamics and kinetics can be normally found in the form of the Haldane relationships
  • the following generalized equation holds true for any reaction mechanism, with catalytic (turnover) term and binding (saturation) term
    • the catalytic term is independent of the reactant concentrations at the reference state
    • 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)

Sampling three ordered mechanisms: Uni-Uni, Bi-Bi and Ter-Ter

  • The average conversion and saturation contributions ratio remains constant for different ΔG r/RT for all the kinetics sampled.
  • The higher the molecularity of the reaction, the higher the contribution of the binding term.
  • 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.
  • the feasible area increases with higher driving force: greater diversity of feasible parameter sets under more thermodynamically favourable conditions

Revealing the impact of thermodynamics on enzyme kinetics

  • uniformly sampled the kinetic space of reactions with same molecularity
  • Ordered vs ping-pong vs random-order
  • calculating the reaction sensitivities for substrates and products @ at different Gibbs free energy differences
  • The substrate and product elasticities strongly depend on the chosen thermodynamic reference state: variable vs constant
    • variable region: linear regime with a steep slope and transition regime
  • Agrees with previous work: almost perfect linear correspondence between the thermodynamic affinity and the reaction flux close to equilibrium: reactions operating close to equilibrium are more susceptible to slight changes in the reactant concentrations
  • (New in this paper): Far from equilibrium revealed both substrate and product elasticities reach a plateau (constant elasticity), A: approximately 0.22 and not 0; P: approaches 0
    • product inhibition is almost negligible for thermodynamically favoured reactions
  • allowable elasticity regions becomes tighter as we move closer to equilibrium
  • the average response of ordered and random mechanisms is more similar than the response of the ping-pong mechanism
  • negative substrate elasticities can be encountered in the random order mechanism for high thermodynamic affinities (apparent substrate inhibition)

Sampling a rare kinetic event: Monomeric cooperativity of mammalian glucokinase

  • Hexokinases I-III are mainly located in the brain, muscle and erythrocytes [43], while hexokinase IV (glucokinase) is primarily located in the liver and pancreatic β-cellsm where it’s the rate-limiting step of glycolysis
  • glucokinase displays a sigmoidal response upon increasing glucose concentration, but hyperbolic behaviour (MM kinetics) with MgATP
  • the enzyme is monomeric, which contradicts standard cooperativity models and thus requires a conceptually different explanation
  • uniformly sampled the kinetic space for this enzyme and counted the frequency of parameter sets displaying positive cooperativity for glucose (estimated Hill coefficient, n H>1). a reference Gibbs free energy difference of -100 kJ/mol was used.
  • 93%: positive cooperativity for glucose. 7% negative => existence of competing parallel pathways in the reaction mechanism
  • The great majority of the models displaying cooperativity for glucose are consistent with a slow transition from the low to the high-affinity enzyme state (98.1%)
  • the sampled kinetics contained the experimentally observed Hill coefficient for this enzyme (1.70)
  • thermodynamically consistent parameter sampling can provide means for effectively accounting for the feasible kinetic space without the need for excessive data

Allosteric regulation: Modelling co-activation of Escherichia coli PEP carboxylase (PEPC)

  • PEPC: CO2-fixing enzymes for plants and bacteria: PEP + HCO3- = OAA + Pi
    • Not PEPCK (PEP carboxykinase), which works in reverse and is involved in gluconeogenesis: OAA + GTP = PEP + GDP + CO2
  • This reaction is highly exergonic (ΔGor = -43.2kJ/mol[52]), making it practically irreversible under physiological conditions
  • Allosterically regulated:
  • sampled the complex regulatory behaviour using the mechanistic information. Also considered the tense form to be capable of performing the reaction, although with a lower activity compared to the relaxed form. Used to describe the ultrasensitive response of PEPC in the presence/absence of acetyl-CoA under changing FBP concentrations
  • accurately described the kinetic behaviour of the PEPC for different FBP concentrations in the presence of acetyl-CoA at physiological concentrations
  • A: a worse performance of our approach is observed in the absence of acetyl-CoA (the framework builds kinetic models around one reference state with acetyl-CoA present)
  • B: perform a rejection step during the sampling so that every accepted parameter set agrees with the experimental data under this condition, typically used in Bayesian Inference by Approximate Bayesian Computation (ABC) methods.
    • most accurate description for this kinetics displays a sigmoidal behaviour as opposed to the observed hyperbolic kinetics.
  • PEPC activation can be achieved by the sole action of FBP or combined with acetyl-CoA
  • C: Hill curves
  • The 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 → ∞
  • acetyl-CoA is the most powerful activator
  • the sampled kinetics is in close agreement with the experimental data
  • 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%


  • Parameterization relies on the integration of the generalized MWC model.
  • Elementary reaction formalism
  • Catalytic rates and thermodynamic constraints between rate parameters
  • Exploration of the feasible kinetic space of models
  • Revealing the impact of the thermodynamic reference state on enzyme kinetics. reaction elasticities
  • Exploration of complex kinetic behaviours can also be achieved


  1. Saa P, Nielsen LK. A general framework for thermodynamically consistent parameterization and efficient sampling of enzymatic reactions. PLoS Comput Biol. 2015;11(4):e1004195. Published 2015 Apr 14. doi:10.1371/journal.pcbi.1004195 PMC4397067 ↩︎