Model Reconstruction for Moment-Based Stochastic Chemical Kinetics

Switch like responses appear as common strategies in the regulation of cellular systems. Here we present a method to characterize bistable regimes in biochemical reaction networks that can be of use to both direct and reverse engineering of biological switches. In the design of a synthetic biological switch, it is important to study the capability for bistability of the underlying biochemical network structure. Chemical Reaction Network Theory (CRNT) may help at this level to decide whether a given network has the capacity for multiple positive equilibria, based on their structural properties. However, in order to build a working switch, we also need to ensure that the bistability property is robust, by studying the conditions leading to the existence of two different steady states. In the reverse engineering of biological switches, knowledge collected about the bistable regimes of the underlying potential model structures can contribute at the model identification stage to a drastic reduction of the feasible region in the parameter space of search. In this work, we make use and extend previous results of the CRNT, aiming not only to discriminate whether a biochemical reaction network can exhibit multiple steady states, but also to determine the regions within the whole space of parameters capable of producing multistationarity. To that purpose we present and justify a condition on the parameters of biochemical networks for the appearance of multistationarity, and propose an efficient and reliable computational method to check its satisfaction through the parameter space.

Automated Design of Enzyme-Driven DNA Circuits

Quantifying effects of stochastic focusing and defocusing on noise propagation in biochemical reaction networks

MotivationStochastic molecular processes are a leading cause of cell-to-cell variability. Their dynamics are often described by continuous-time discrete-state Markov chains and simulated using stochastic simulation algorithms. As these stochastic simulations are computationally demanding, ordinary differential equation models for the dynamics of the statistical moments have been developed. The number of state variables of these approximating models, however, grows at least quadratically with the number of biochemical species. This limits their application to small- and medium-sized processes.ResultsIn this article, we present a scalable moment-closure approximation (sMA) for the simulation of statistical moments of large-scale stochastic processes. The sMA exploits the structure of the biochemical reaction network to reduce the covariance matrix. We prove that sMA yields approximating models whose number of state variables depends predominantly on local properties, i.e. the average node degree of the reaction network, instead of the overall network size. The resulting complexity reduction is assessed by studying a range of medium- and large-scale biochemical reaction networks. To evaluate the approximation accuracy and the improvement in computational efficiency, we study models for JAK2/STAT5 signalling and NFκB signalling. Our method is applicable to generic biochemical reaction networks and we provide an implementation, including an SBML interface, which renders the sMA easily accessible.Availability and implementationThe sMA is implemented in the open-source MATLAB toolbox CERENA and is available from https://github.com/CERENADevelopers/CERENA.Supplementary informationSupplementary data are available at Bioinformatics online.

BackgroundStochastic modeling and simulation provide powerful predictive methods for the intrinsic understanding of fundamental mechanisms in complex biochemical networks. Typically, such mathematical models involve networks of coupled jump stochastic processes with a large number of parameters that need to be suitably calibrated against experimental data. In this direction, the parameter sensitivity analysis of reaction networks is an essential mathematical and computational tool, yielding information regarding the robustness and the identifiability of model parameters. However, existing sensitivity analysis approaches such as variants of the finite difference method can have an overwhelming computational cost in models with a high-dimensional parameter space.ResultsWe develop a sensitivity analysis methodology suitable for complex stochastic reaction networks with a large number of parameters. The proposed approach is based on Information Theory methods and relies on the quantification of information loss due to parameter perturbations between time-series distributions. For this reason, we need to work on path-space, i.e., the set consisting of all stochastic trajectories, hence the proposed approach is referred to as “pathwise”. The pathwise sensitivity analysis method is realized by employing the rigorously-derived Relative Entropy Rate, which is directly computable from the propensity functions. A key aspect of the method is that an associated pathwise Fisher Information Matrix (FIM) is defined, which in turn constitutes a gradient-free approach to quantifying parameter sensitivities. The structure of the FIM turns out to be block-diagonal, revealing hidden parameter dependencies and sensitivities in reaction networks.ConclusionsAs a gradient-free method, the proposed sensitivity analysis provides a significant advantage when dealing with complex stochastic systems with a large number of parameters. In addition, the knowledge of the structure of the FIM can allow to efficiently address questions on parameter identifiability, estimation and robustness. The proposed method is tested and validated on three biochemical systems, namely: (a) a protein production/degradation model where explicit solutions are available, permitting a careful assessment of the method, (b) the p53 reaction network where quasi-steady stochastic oscillations of the concentrations are observed, and for which continuum approximations (e.g. mean field, stochastic Langevin, etc.) break down due to persistent oscillations between high and low populations, and (c) an Epidermal Growth Factor Receptor model which is an example of a high-dimensional stochastic reaction network with more than 200 reactions and a corresponding number of parameters.

Succinct, declarative, and domain-specific modeling languages have many advantages when creating simulation models. However, it is often challenging to efficiently execute models defined in such languages. We use code generation for model-specific simulators. Code generation has been successfully applied for high-performance algorithms in many application domains. By generating tailored simulators for specific simulation models defined in a domain specific language, we get the best of both worlds: a succinct, declarative and formal presentation of the model and an efficient execution. We illustrate this based on a simple domain-specific language for biochemical reaction networks as well as on the network representation of the established BioNetGen language. We implement two approaches adopting the same simulation algorithms: one generic simulator that parses models at runtime and one generator that produces a simulator specialized to a given model based on partial evaluation and code generation. Akin to profile-guided optimization we also use dynamic execution of the model to further optimize the simulators. The performance of the approaches is carefully benchmarked using representative models of small to mid-sized biochemical reaction networks. The generic simulator achieves a performance similar to state of the art simulators in the domain, whereas the specialized simulator outperforms established simulation algorithms with a speedup of more than an order of magnitude. This repository contains the code generation software as described in the 2020 PADS paper Partial evaluation via code generation for static stochastic reaction network models.

Formal methods have enabled breakthroughs in many fields, such as in hardware verification, machine learning and biological systems. The key object of interest in systems biology, synthetic biology, and molecular programming is chemical reaction networks (CRNs) which formalizes coupled chemical reactions in a well-mixed solution. CRNs are pivotal for our understanding of biological regulatory and metabolic networks, as well as for programming engineered molecular behavior. Although it is clear that small CRNs are capable of complex dynamics and computational behavior, it remains difficult to explore the space of CRNs in search for desired functionality. We use Alloy, a tool for expressing structural constraints and behavior in software systems, to enumerate CRNs with declaratively specified properties. We show how this framework can enumerate CRNs with a variety of structural constraints including biologically motivated catalytic networks and metabolic networks, and seesaw networks motivated by DNA nanotechnology. We also use the framework to explore analog function computation in rate-independent CRNs. By computing the desired output value with stoichiometry rather than with reaction rates (in the sense that X → Y+Y computes multiplication by 2), such CRNs are completely robust to the choice of reaction rates or rate law. We find the smallest CRNs computing the max, minmax, abs and ReLU (rectified linear unit) functions in a natural subclass of rate-independent CRNs where rate-independence follows from structural network properties.

Reaction networks in thermodynamic equilibrium under isothermal and isobaric conditions minimize the Gibbs free energy, but chemical reactions in living organisms operate typically far from equilibrium. Currently, there is no general optimization principle for nonequilibrium systems which can be used in the analysis of biochemical networks. Motivated by the avalailabity of whole genome reconstructions of metabolic reactions, the thermodynamics of biochemical stoichiometric networks has made significant progress in the last decade. These include the consistent formulation of conservation conditions resembling Kirchhoff’s law for electrical networks. In addition, Beard and Qian suggested that the flow force relationship Δμ = RT log(J + /J − ) between the forward and backward fluxes J + and J − and the chemical potential difference of a chemical reaction can be extended from mass action kinetics to more general reactions schemes. In this tutorial review we summarize the recent literature on reaction network thermodynamics and discuss its implications to the analysis of large biochemical systems. In addition, we discuss some recent work on flow-force relationships and global variational principles characterizing nonequilibrium steady states of reaction networks.

To elucidate the dynamic features of a biologically relevant large-scale reaction network, we constructed a computational model of minimal protein synthesis consisting of 241 components and 968 reactions that synthesize the Met-Gly-Gly (MGG) peptide based on an Escherichia coli-based reconstituted in vitro protein synthesis system. We performed a simulation using parameters collected primarily from the literature and found that the rate of MGG peptide synthesis becomes nearly constant in minutes, thus achieving a steady state similar to experimental observations. In addition, concentration changes to 70% of the components, including intermediates, reached a plateau in a few minutes. However, the concentration change of each component exhibits several temporal plateaus, or a quasi-stationary state (QSS), before reaching the final plateau. To understand these complex dynamics, we focused on whether the components reached a QSS, mapped the arrangement of components in a QSS in the entire reaction network structure, and investigated time-dependent changes. We found that components in a QSS form clusters that grow over time but not in a linear fashion, and that this process involves the collapse and regrowth of clusters before the formation of a final large single cluster. These observations might commonly occur in other large-scale biological reaction networks. This developed analysis might be useful for understanding large-scale biological reactions by visualizing complex dynamics, thereby extracting the characteristics of the reaction network, including phase transitions.

Biochemical reaction networks are expected to encode an efficient representation of the function of cells in a variable environment. It is thus important to see how these networks do learn and implement such representations. The first step in this direction is to characterize the function and learning capabilities of basic artificial reaction networks. In this study, we consider multilayer networks of reversible reactions that connect two layers of signal and response species through an intermediate layer of hidden species. We introduce a stochastic learning algorithm that updates the reaction rates based on the correlation values between reaction products and responses. Our findings indicate that the function of networks with random reaction rates, as well as their learning capacity for random signal-response activities, are critically determined by the number of reactants and reaction products. Moreover, the stored patterns exhibit different levels of robustness and qualities as the reaction rates deviate from their optimal values in a stochastic model of defect evolution. These findings can help suggest network modules that are better suited to specific functions, such as amplifiers or dampeners, or to the learning of biologically relevant signal-response activities.

Abstract. One notable observation of enzymatic chemical reactions is that, for a given abundance of enzymes and substrates, temperature increases cause reaction rates to first increase consistent with the Arrhenius relationship, then plateau, and finally fall off quickly to zero at high temperatures. While many mathematical functions have been used to describe this pattern, we here propose a chemical kinetics theory which successfully replicates this observation and provides insights into the processes responsible for these dynamics. The chemical kinetics theory combines the law of mass action, von Smoluchowski's diffusion-limited chemical reaction theory, and Eyring's transition state theory. This new theory reveals that the thermally reversible enzyme denaturation ensured by the ceaseless thermal motion of molecules and ions in an enzyme solution explains the plateau and subsequent decrease in chemical reaction rates with increasing temperature. The temperature-dependent affinity parameter (K) that relates enzymes and substrates through their binding also affects the shape of the emergent temperature response. We demonstrate that with an increase in substrate availability, K shifts the optimal temperature, where reaction rates plateau, towards higher values. Further, we show that the chemical kinetics theory accurately represents 12 sets of published enzyme assay data and includes the popular mechanistic model by Ratkowsky et al. (2005) as a special case. Given its good performance and solid theoretical underpinning, we believe this new theory will facilitate the construction of more mechanistic-based environmental biogeochemical models.

Motivation: Biochemical reaction networks in the form of coupled ordinary differential equations (ODEs) provide a powerful modeling tool for understanding the dynamics of biochemical processes. During the early phase of modeling, scientists have to deal with a large pool of competing nonlinear models. At this point, discrimination experiments can be designed and conducted to obtain optimal data for selecting the most plausible model. Since biological ODE models have widely distributed parameters due to, e.g. biologic variability or experimental variations, model responses become distributed. Therefore, a robust optimal experimental design (OED) for model discrimination can be used to discriminate models based on their response probability distribution functions (PDFs).Results: In this work, we present an optimal control-based methodology for designing optimal stimulus experiments aimed at robust model discrimination. For estimating the time-varying model response PDF, which results from the nonlinear propagation of the parameter PDF under the ODE dynamics, we suggest using the sigma-point approach. Using the model overlap (expected likelihood) as a robust discrimination criterion to measure dissimilarities between expected model response PDFs, we benchmark the proposed nonlinear design approach against linearization with respect to prediction accuracy and design quality for two nonlinear biological reaction networks. As shown, the sigma-point outperforms the linearization approach in the case of widely distributed parameter sets and/or existing multiple steady states. Since the sigma-point approach scales linearly with the number of model parameter, it can be applied to large systems for robust experimental planning.Availability: An implementation of the method in MATLAB/AMPL is available at http://www.uni-magdeburg.de/ivt/svt/person/rf/roed.html.Contact: flassig@mpi-magdeburg.mpg.deSupplementary information: Supplementary data are are available at Bioinformatics online.

The evolutionary origin of universal statistics in biochemical reaction networks is studied, to explain the power-law distribution of reaction links and the power-law distributions of chemical abundances. Using cell models with catalytic reaction networks, we confirmed that the power-law distribution in abundances of chemicals emerges by the selection of cells with higher growth speeds, as suggested in our previous study. Through the further evolution, this inhomogeneity in chemical abundances is shown to be embedded in the distribution of links, leading to the power-law distribution. These findings provide novel insights into the nature of network evolution in living cells.KeywordsReaction PathReaction NetworkGrowth SpeedBiochemical NetworkChemical AbundanceThese keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

The rejection-based simulation technique has been applying to improve the computational efficiency of the stochastic simulation algorithm (SSA) in simulating large reaction networks, which are required for a thorough understanding of biological systems. We compare two recently proposed simulation methods, namely the composition-rejection algorithm (SSA-CR) and the rejection-based SSA (RSSA), aiming for this purpose. We discuss the right interpretation of the rejection-based technique used in these algorithms in order to make an informed choice when dealing with different aspects of biochemical networks. We provide the theoretical analysis as well as the detailed runtime comparison of these algorithms on concrete biological models. We highlight important factors that are omitted in previous analysis of these algorithms. The numerical comparison shows that for reaction networks where the search cost is expensive then SSA-CR is more efficient, and for reaction networks where the update cost is dominant, often the case in practice, then RSSA should be the choice.

Weak reversibility is a crucial structural property of chemical reaction networks (CRNs) with mass action kinetics, because it has major implications related to the existence, uniqueness and stability of equilibrium points and to the boundedness of solutions. In this paper, we present two new algorithms to find dynamically equivalent weakly reversible realizations of a given CRN. They are based on linear programming and thus have polynomial time-complexity. Hence, these algorithms can deal with large-scale biochemical reaction networks, too. Furthermore, one of the methods is able to deal with linearly conjugate networks, too.