Optimal enzyme utilization suggests that concentrations and thermodynamics determine binding mechanisms and enzyme saturations

Living organisms continuously adapt to genetic and environmental perturbations, and understanding these responses requires characterization of enzymatic reaction mechanisms and kinetic properties across metabolic networks. Kinetic models relate reaction rates to metabolite concentrations and kinetic parameters, yet comprehensive experimental kinetic datasets are scarce even for central metabolism, necessitating estimation or sampling methods. Because enzyme parameters are outcomes of natural selection, appropriate evolutionary fitness functions can illuminate design principles governing catalytic rates. Prior work suggests evolutionary pressure toward efficient enzyme utilization, with the ratio of specific flux to enzyme concentration (vnet/Etot) as a key determinant. Modeling efforts have explored how evolution shapes enzyme and metabolite concentrations and kinetic parameters under various objectives and optimization paradigms. However, existing approaches either simplify kinetics or are limited to ordered mechanisms and local, nonlinear optimizations, making extension to complex random or ping-pong mechanisms difficult. This study addresses these gaps by developing OpEn, a mixed-integer linear programming (MILP) framework that maximizes net steady-state flux at fixed total enzyme, to estimate optimal kinetic parameters and to reveal how thermodynamic displacements, enzyme saturations, and binding mechanisms depend on substrate/product concentrations and thermodynamics.

The authors review evolutionary optimization principles applied to enzyme kinetics and metabolism, from isolated enzymes to network-level objectives (e.g., minimizing fluxes or concentrations, maximizing thermodynamic efficiency). Meta-analyses indicate most enzymes are evolved to be “good enough,” not perfect, but the constraints and pressures shaping kinetic parameters remain incompletely understood. Population-based methods have explored catalytic optimality but often neglect detailed kinetics and thermodynamics or rely on simplified mechanisms and extensive hyperparameter tuning without guarantees of optimality. Heinrich, Klipp, and Wilhelm developed nonlinear optimization approaches for unbranched, ordered mechanisms and showed that reactant concentrations partition parameter space into regions with distinct optimal binding characteristics and that reactant concentrations and Michaelis constants co-vary over evolutionary time. These studies, however, are restricted to ordered mechanisms and require pre-deriving candidate solution structures followed by local optimization, limiting scalability to diverse mechanisms found across metabolic networks. The literature thus motivates a general, computationally efficient, globally optimal framework capable of handling arbitrary elementary mechanisms with thermodynamic rigor.

OpEn formulates optimal enzyme utilization as a MILP that maximizes net steady-state flux (vnet) for a given total enzyme concentration. Inputs: (i) the elementary reaction mechanism, (ii) intracellular substrate/product concentrations, and (iii) thermodynamic properties (standard Gibbs free energy; Keq or overall displacement Γ). Variables include elementary rate constants (forward/backward), elementary thermodynamic displacements (γi), and enzyme state abundances (ei). Constraints implement: (1) quasi–steady-state mass-action kinetics decomposing reversible steps into forward and backward fluxes, (2) conservation of total enzyme states (Σei=1), (3) thermodynamic consistency linking forward/backward flux ratios to displacements (Γ=∏Yi over steps or fundamental cycles) with Γ≤γi≤1 for forward operation, and (4) biophysical upper bounds on rate constants (bimolecular 10^8–10^10 M−1 s−1; monomolecular 10^4–10^6 s−1). Variables and parameters are normalized to yield dimensionless quantities: rate constants by their biophysical limits, metabolite concentrations by a characteristic concentration [C]ch, enzyme states by total enzyme concentration. To remove nonlinearities, bilinear products of rate constants and enzyme states are replaced by new variables (zij) with bounding constraints; remaining nonlinearities from γi are handled by approximating selected independent elementary displacements or their mechanistically meaningful combinations via piecewise-constant (0th order) functions with binary expansions, ensuring the overall thermodynamic constraint (product over cycles equals Γ) is satisfied. The resulting mixed bilinear terms (binary×continuous) are linearized using Petersen’s scheme, yielding a MILP that ensures global optimality and can enumerate alternative solutions. For random-ordered mechanisms, branch fluxes are aggregated to vnet and a splitting ratio α is defined as the fraction of net flux through one branch (e.g., upper branch where substrate A binds first). Variability analysis quantifies multiplicity at the optimal vnet, and sampling (e.g., ACHR or optGpSampler) of the linearized problem with fixed γi allows exploration of alternative modes and suboptimal landscapes (with constraint vnet≥c·vnet*). Macroscopic parameters (KM, kcat) can be back-calculated from microscopic solutions via Cleland’s notation or in silico initial-rate experiments. The framework is implemented in Python using optlang with commercial MILP solvers (CPLEX/Gurobi), and code/data are publicly available.

• Michaelis-Menten mechanism (three-step reversible): OpEn reproduces known optimal rate constant patterns and shows that optimal enzyme saturation depends strongly on reactant concentrations, increasing with substrate and product concentrations relative to [C]ch. The optimal saturation increases rapidly when both reactants are below [C]ch and more slowly at higher concentrations. Notably, this saturation behavior appears largely independent of the overall thermodynamic displacement Γ. Analysis across saturation regimes shows systematic redistribution of thermodynamic driving forces: at low saturation, ~60% of ΔG′ drives substrate association, most of the remainder drives product dissociation, and little drives the interconversion; at intermediate saturation, driving force shifts toward the interconversion with reduced substrate association drive; at high saturation, forces are more evenly distributed and free enzyme becomes the least abundant species, though a minimal free enzyme pool remains necessary for optimality.
• Ordered Bi-Uni mechanism: At optimality, saturation increases with the concentration of the substrate that binds first to the enzyme, whereas the second substrate concentration has little effect on saturation. Product concentration still increases overall saturation (trend similar to Michaelis-Menten).
• Random-ordered Bi-Uni mechanism (branched binding): Optimal enzyme utilization exhibits concentration-dependent binding preferences. The splitting ratio α (fraction of net flux through the branch where A binds first) typically lies between 0.3 and 0.7 under physiological conditions (e.g., P=0.1 mM), indicating that a random-ordered mechanism is superior to any purely ordered mechanism. Key phenomenology: (i) antisymmetry: interchanging substrate concentrations (A,B ↔ B,A) swaps α such that αBA=1−αAB; (ii) uniqueness vs. flexibility: for unequal substrates (A≠B) the optimal α is unique; for equal substrates (A=B) α is flexible around 0.5, yielding alternative elementary constant sets with identical vnet, saturations, and thermodynamic force distributions for non-branch steps; (iii) the relative product concentration governs which substrate binds first: when the least abundant substrate is below P, the most abundant substrate binds first; when the least abundant substrate is comparable to or above P, the least abundant substrate binds first. At low P (e.g., P=0.1), the lowest substrate dictates binding preference; at high P (e.g., P=5), the most abundant substrate dictates preference. Symmetric operating points yield symmetric modes (interchange of EA and EB occupancies and branch fluxes), identical overall saturation, and conserved substrate- plus product-bound fractions across symmetric solutions.
• General: The framework ensures global optima for vnet, characterizes enzyme state distributions and thermodynamic force allocations at optimality, and identifies condition-specific unique or alternative operating modes dependent on reactant concentrations and thermodynamic constraints.

The study addresses how evolutionary pressures for efficient resource use translate into optimal enzyme utilization by maximizing vnet at fixed Etot under biophysical and thermodynamic constraints. Findings reveal that optimal saturation and binding mechanisms are primarily determined by reactant concentrations relative to a characteristic scale and by their relation to the thermodynamic driving forces, rather than by Γ alone. For Michaelis-Menten enzymes, the required saturation regime shifts with reactant availability, with thermodynamic driving forces reallocated across association, interconversion, and dissociation steps to sustain maximal throughput. For multi-substrate enzymes, concentration-dependent and product-relative effects determine binding order, with random-ordered mechanisms outperforming ordered ones under physiological conditions, clarifying why many enzymes may have evolved flexible binding sequences. The MILP formulation provides mechanistic insight into enzyme state occupancies, flux partitioning, and force distributions at optimality and can systematically enumerate alternative optimal kinetic designs when they exist. Practically, OpEn can fill kinetic gaps in metabolic models by estimating condition-specific optimal parameters from mechanism, thermodynamics, and metabolite levels, integrate with thermodynamics-based flux analysis to obtain feasible concentration/flux profiles, and explore suboptimal fitness landscapes to understand trade-offs shaping moderately efficient enzymes.

OpEn, a MILP-based framework, generalizes catalytic optimality analysis to arbitrary elementary mechanisms with explicit thermodynamic consistency, enabling estimation of optimal kinetic parameters, enzyme state distributions, and thermodynamic force allocations from minimal inputs. It reproduces classic results for ordered mechanisms and extends them to random-ordered multi-substrate systems, uncovering that optimal saturations and binding preferences are determined by reactant concentrations and their relation to product levels, with random mechanisms often optimal in bimolecular reactions under physiological conditions. The approach offers condition-specific theoretical upper bounds on catalytic efficiency, supports variability and sampling to map suboptimal fitness landscapes, and can guide enzyme engineering by revealing mode-of-operation signatures at optimality. Future work includes extending to crowded intracellular environments, integrating additional evolutionary objectives (e.g., protein cost, robustness), scaling to pathways and networks (linear complexity in reaction count), and leveraging combined omics-constrained models to estimate context-specific kinetic parameters in vivo.

The analyses focus on maximal net steady-state flux (optimal state), while real enzymes may operate suboptimally; although the framework can sample suboptimal solutions, these were not the main focus. The formulation assumes dilute-solution kinetics and quasi-steady state; crowding effects and rapid dynamics are not explicitly modeled here (though extensions are envisioned). Elementary displacements are approximated piecewise-constantly, introducing discretization error (controlled by resolution). Upper bounds for mono- and bimolecular rate constants are applied uniformly in most cases and do not distinguish specific mono-molecular step types, though this can be generalized. The approach requires specification of reaction mechanism, thermodynamics, and metabolite concentrations, which may be uncertain and context-dependent.