Understanding the dynamics of biological systems in evolving environments is a challenge due to their scale and complexity. Here, we present a computational framework for the timescale decomposition of biochemical reaction networks to distill essential patterns from their intricate dynamics. This approach identifies timescale hierarchies, concentration pools, and coherent structures from time-series data, providing a system-level description of reaction networks at physiologically important timescales. We apply this technique to kinetic models of hypothetical and biological pathways, validating it by reproducing analytically characterized or previously known concentration pools of these pathways. Moreover, by analyzing the timescale hierarchy of the glycolytic pathway, we elucidate the connections between the stoichiometric and dissipative structures of reaction networks and the temporal organization of coherent structures. Specifically, we show that glycolysis is a cofactor-driven pathway, the slowest dynamics of which are described by a balance between high-energy phosphate bond and redox trafficking. Overall, this approach provides more biologically interpretable characterizations of network dynamics than large-scale kinetic models, thus facilitating model reduction and personalized medicine applications. IMPORTANCE Complex interactions within interconnected biochemical reaction networks enable cellular responses to a wide range of unpredictable environmental perturbations. Understanding how biological functions arise from these intricate interactions has been a long-standing problem in biology. Here, we introduce a computational approach to dissect complex biological systems dynamics in evolving environments. This approach characterizes the timescale hierarchies of complex reaction networks, offering a system-level understanding at physiologically relevant timescales. Analyzing various hypothetical and biological pathways, we show how stoichiometric properties shape the way energy is dissipated throughout reaction networks. Notably, we establish that glycolysis operates as a cofactor-driven pathway, where the slowest dynamics are governed by a balance between high-energy phosphate bonds and redox trafficking. This approach enhances our understanding of network dynamics and facilitates the development of reduced-order kinetic models with biologically interpretable components.

Mathematical models of metabolic networks utilize simulation to study system-level mechanisms and functions. Various approaches have been used to model the steady state behavior of metabolic networks using genome-scale reconstructions, but formulating dynamic models from such reconstructions continues to be a key challenge. Here, we present the Mass Action Stoichiometric Simulation Python (MASSpy) package, an open-source computational framework for dynamic modeling of metabolism. MASSpy utilizes mass action kinetics and detailed chemical mechanisms to build dynamic models of complex biological processes. MASSpy adds dynamic modeling tools to the COnstraint-Based Reconstruction and Analysis Python (COBRApy) package to provide an unified framework for constraint-based and kinetic modeling of metabolic networks. MASSpy supports high-performance dynamic simulation through its implementation of libRoadRunner: the Systems Biology Markup Language (SBML) simulation engine. Three examples are provided to demonstrate how to use MASSpy: (1) a validation of the MASSpy modeling tool through dynamic simulation of detailed mechanisms of enzyme regulation; (2) a feature demonstration using a workflow for generating ensemble of kinetic models using Monte Carlo sampling to approximate missing numerical values of parameters and to quantify biological uncertainty, and (3) a case study in which MASSpy is utilized to overcome issues that arise when integrating experimental data with the computation of functional states of detailed biological mechanisms. MASSpy represents a powerful tool to address challenges that arise in dynamic modeling of metabolic networks, both at small and large scales.

Allosteric regulation has traditionally been described by mathematically-complex allosteric rate laws in the form of ratios of polynomials derived from the application of simplifying kinetic assumptions. Alternatively, an approach that explicitly describes all known ligand-binding events requires no simplifying assumptions while allowing for the computation of enzymatic states. Here, we employ such a modeling approach to examine the "catalytic potential" of an enzyme-an enzyme's capacity to catalyze a biochemical reaction. The catalytic potential is the fundamental result of multiple ligand-binding events that represents a "tug of war" among the various regulators and substrates within the network. This formalism allows for the assessment of interacting allosteric enzymes and development of a network-level understanding of regulation. We first define the catalytic potential and use it to characterize the response of three key kinases (hexokinase, phosphofructokinase, and pyruvate kinase) in human red blood cell glycolysis to perturbations in ATP utilization. Next, we examine the sensitivity of the catalytic potential by using existing personalized models, finding that the catalytic potential allows for the identification of subtle but important differences in how individuals respond to such perturbations. Finally, we explore how the catalytic potential can help to elucidate how enzymes work in tandem to maintain a homeostatic state. Taken together, this work provides an interpretation and visualization of the dynamic interactions and network-level effects of interacting allosteric enzymes.

Knowing the catalytic turnover numbers of enzymes is essential for understanding the growth rate, proteome composition, and physiology of organisms, but experimental data on enzyme turnover numbers is sparse and noisy. Here, we demonstrate that machine learning can successfully predict catalytic turnover numbers in Escherichia coli based on integrated data on enzyme biochemistry, protein structure, and network context. We identify a diverse set of features that are consistently predictive for both in vivo and in vitro enzyme turnover rates, revealing novel protein structural correlates of catalytic turnover. We use our predictions to parameterize two mechanistic genome-scale modelling frameworks for proteome-limited metabolism, leading to significantly higher accuracy in the prediction of quantitative proteome data than previous approaches. The presented machine learning models thus provide a valuable tool for understanding metabolism and the proteome at the genome scale, and elucidate structural, biochemical, and network properties that underlie enzyme kinetics.

Computational models have great potential to accelerate bioscience, bioengineering, and medicine. However, it remains challenging to reproduce and reuse simulations, in part, because the numerous formats and methods for simulating various subsystems and scales remain siloed by different software tools. For example, each tool must be executed through a distinct interface. To help investigators find and use simulation tools, we developed BioSimulators (https://biosimulators.org), a central registry of the capabilities of simulation tools and consistent Python, command-line and containerized interfaces to each version of each tool. The foundation of BioSimulators is standards, such as CellML, SBML, SED-ML and the COMBINE archive format, and validation tools for simulation projects and simulation tools that ensure these standards are used consistently. To help modelers find tools for particular projects, we have also used the registry to develop recommendation services. We anticipate that BioSimulators will help modelers exchange, reproduce, and combine simulations.