CERENA: ChEmical REaction Network Analyzer--A Toolbox for the Simulation and Analysis of Stochastic Chemical Kinetics

Catalyst: Fast and flexible modeling of reaction networks

We introduce Catalyst.jl, a flexible and feature-filled Julia library for modeling and high-performance simulation of chemical reaction networks (CRNs). Catalyst supports simulating stochastic chemical kinetics (jump process), chemical Langevin equation (stochastic differential equation), and reaction rate equation (ordinary differential equation) representations for CRNs. Through comprehensive benchmarks, we demonstrate that Catalyst simulation runtimes are often one to two orders of magnitude faster than other popular tools. More broadly, Catalyst acts as both a domain-specific language and an intermediate representation for symbolically encoding CRN models as Julia-native objects. This enables a pipeline of symbolically specifying, analyzing, and modifying CRNs; converting Catalyst models to symbolic representations of concrete mathematical models; and generating compiled code for numerical solvers. Leveraging ModelingToolkit.jl and Symbolics.jl, Catalyst models can be analyzed, simplified, and compiled into optimized representations for use in numerical solvers. Finally, we demonstrate Catalyst's broad extensibility and composability by highlighting how it can compose with a variety of Julia libraries, and how existing open-source biological modeling projects have extended its intermediate representation.

Shadow enhancers mediate trade-offs between transcriptional noise and fidelity

Enhancers are stretches of regulatory DNA that bind transcription factors (TFs) and regulate the expression of a target gene. Shadow enhancers are two or more enhancers that regulate the same target gene in space and time and are associated with most animal developmental genes. These multi-enhancer systems can drive more consistent transcription than single enhancer systems. Nevertheless, it remains unclear why shadow enhancer TF binding sites are distributed across multiple enhancers rather than within a single large enhancer. Here, we use a computational approach to study systems with varying numbers of TF binding sites and enhancers. We employ chemical reaction networks with stochastic dynamics to determine the trends in transcriptional noise and fidelity, two key performance objectives of enhancers. This reveals that while additive shadow enhancers do not differ in noise and fidelity from their single enhancer counterparts, sub- and superadditive shadow enhancers have noise and fidelity trade-offs not available to single enhancers. We also use our computational approach to compare the duplication and splitting of a single enhancer as mechanisms for the generation of shadow enhancers and find that the duplication of enhancers can decrease noise and increase fidelity, although at the metabolic cost of increased RNA production. A saturation mechanism for enhancer interactions similarly improves on both of these metrics. Taken together, this work highlights that shadow enhancer systems may exist for several reasons: genetic drift or the tuning of key functions of enhancers, including transcription fidelity, noise and output.

BioSANS: A software package for symbolic and numeric biological simulation

Modeling biochemical systems can provide insights into behaviors that are difficult to observe or understand. It requires software, programming, and understanding of the system to build a model and study it. Softwares exist for systems biology modeling, but most support only certain types of modeling tasks. Desirable features including ease in preparing input, symbolic or analytical computation, parameter estimation, graphical user interface, and systems biology markup language (SBML) support are not seen concurrently in one software package. In this study, we developed a python-based software that supports these features, with both deterministic and stochastic propagations. The software can be used by graphical user interface, command line, or as a python import. We also developed a semi-programmable and intuitively easy topology input method for the biochemical reactions. We tested the software with semantic and stochastic SBML test cases. Tests on symbolic solution and parameter estimation were also included. The software we developed is reliable, well performing, convenient to use, and compliant with most of the SBML tests. So far it is the only systems biology software that supports symbolic, deterministic, and stochastic modeling in one package that also features parameter estimation and SBML support. This work offers a comprehensive set of tools and allows for better availability and accessibility for studying kinetics and dynamics in biochemical systems.

Automated Exploration of Prebiotic Chemical Reaction Space: Progress and Perspectives

Prebiotic chemistry often involves the study of complex systems of chemical reactions that form large networks with a large number of diverse species. Such complex systems may have given rise to emergent phenomena that ultimately led to the origin of life on Earth. The environmental conditions and processes involved in this emergence may not be fully recapitulable, making it difficult for experimentalists to study prebiotic systems in laboratory simulations. Computational chemistry offers efficient ways to study such chemical systems and identify the ones most likely to display complex properties associated with life. Here, we review tools and techniques for modelling prebiotic chemical reaction networks and outline possible ways to identify self-replicating features that are central to many origin-of-life models.

Improved estimations of stochastic chemical kinetics by finite-state expansion

Stochastic reaction networks are a fundamental model to describe interactions between species where random fluctuations are relevant. The master equation provides the evolution of the probability distribution across the discrete state space consisting of vectors of population counts for each species. However, since its exact solution is often elusive, several analytical approximations have been proposed. The deterministic rate equation (DRE) gives a macroscopic approximation as a compact system of differential equations that estimate the average populations for each species, but it may be inaccurate in the case of nonlinear interaction dynamics. Here we propose finite-state expansion (FSE), an analytical method mediating between the microscopic and the macroscopic interpretations of a stochastic reaction network by coupling the master equation dynamics of a chosen subset of the discrete state space with the mean population dynamics of the DRE. An algorithm translates a network into an expanded one where each discrete state is represented as a further distinct species. This translation exactly preserves the stochastic dynamics, but the DRE of the expanded network can be interpreted as a correction to the original one. The effectiveness of FSE is demonstrated in models that challenge state-of-the-art techniques due to intrinsic noise, multi-scale populations and multi-stability.

MomentClosure.jl: automated moment closure approximations in Julia

SUMMARY

MomentClosure.jl is a Julia package providing automated derivation of the time-evolution equations of the moments of molecule numbers for virtually any chemical reaction network using a wide range of moment closure approximations. It extends the capabilities of modelling stochastic biochemical systems in Julia and can be particularly useful when exact analytic solutions of the chemical master equation are unavailable and when Monte Carlo simulations are computationally expensive.

AVAILABILITY AND IMPLEMENTATION

MomentClosure.jl is freely accessible under the MIT licence. Source code and documentation are available at https://github.com/augustinas1/MomentClosure.jl.

Compartor: A toolbox for the automatic generation of moment equations for dynamic compartment populations

Summary

Many biochemical processes in living organisms take place inside compartments that can interact with each other and remodel over time. In a recent work, we have shown how the stochastic dynamics of a compartmentalized biochemical system can be effectively studied using moment equations. With this technique, the time evolution of a compartment population is summarized using a finite number of ordinary differential equations, which can be analyzed very efficiently. However, the derivation of moment equations by hand can become time-consuming for systems comprising multiple reactants and interactions. Here we present Compartor, a toolbox that automatically generates the moment equations associated with a user-defined compartmentalized system. Through the moment equation method, Compartor renders the analysis of stochastic population models accessible to a broader scientific community.

Availability and implementation

Compartor is provided as a Python package and is available at https://pypi.org/project/compartor/. Source code and usage tutorials for Compartor are available at https://github.com/zechnerlab/Compartor.

Comparison of Queueing Data-Structures for Kinetic Monte Carlo Simulations of Heterogeneous Catalysts

On-lattice kinetic Monte Carlo (KMC) is a computational method used to simulate (among others) physicochemical processes on catalytic surfaces. The KMC algorithm propagates the system through discrete configurations by selecting (with the use of random numbers) the next elementary process to be simulated, e.g., adsorption, desorption, diffusion, or reaction. An implementation of such a selection procedure is the first-reaction method in which all realizable elementary processes are identified and assigned a random occurrence time based on their rate constant. The next event to be executed will then be the one with the minimum interarrival time. Thus, a fast and efficient algorithm for selecting the most imminent process and performing all of the necessary updates on the list of realizable processes post execution is of great importance. In the current work, we implement five data-structures to handle the elementary process queue during a KMC run: an unsorted list, a binary heap, a pairing heap, a one-way skip list, and finally, a novel two-way skip list with a mapping array specialized for KMC simulations. We also investigate the effect of compiler optimizations on the performance of these data-structures on three benchmark models, capturing CO oxidation, a simplified water gas shift mechanism, and a temperature-programmed desorption run. Excluding the least efficient and impractical for large-problems unsorted list, we observe a 3× speedup of the binary or pairing heaps (most efficient) compared to the one-way skip list (least efficient). Compiler optimizations deliver a speedup of up to 1.8×. These benchmarks provide valuable insight into the importance of, often-overlooked, implementation-related aspects of KMC simulations, such as the queueing data-structures. Our results could be particularly useful in guiding the choice of data-structures and algorithms that would minimize the computational cost of large-scale simulations.

Slack reactants: A state-space truncation framework to estimate quantitative behavior of the chemical master equation

State space truncation methods are widely used to approximate solutions of the chemical master equation. While most methods of this kind focus on truncating the state space directly, in this work, we propose modifying the underlying chemical reaction network by introducing slack reactants that indirectly truncate the state space. More specifically, slack reactants introduce an expanded chemical reaction network and impose a truncation scheme based on desired mass conservation laws. This network structure also allows us to prove inheritance of special properties of the original model, such as irreducibility and complex balancing. We use the network structure imposed by slack reactants to prove the convergence of the stationary distribution and first arrival times. We then provide examples comparing our method with the stationary finite state projection and finite buffer methods. Our slack reactant system appears to be more robust than some competing methods with respect to calculating first arrival times.

Robust calibration of hierarchical population models for heterogeneous cell populations

Cellular heterogeneity is known to have important effects on signal processing and cellular decision making. To understand these processes, multiple classes of mathematical models have been introduced. The hierarchical population model builds a novel class which allows for the mechanistic description of heterogeneity and explicitly takes into account subpopulation structures. However, this model requires a parametric distribution assumption for the cell population and, so far, only the normal distribution has been employed. Here, we incorporate alternative distribution assumptions into the model, assess their robustness against outliers and evaluate their influence on the performance of model calibration in a simulation study and a real-world application example. We found that alternative distributions provide reliable parameter estimates even in the presence of outliers, and can in fact increase the convergence of model calibration.

A Reaction-Based Model of the State Space of Chemical Reaction Systems Enables Efficient Simulations

The choice of the state space representation of a system can turn into a prominent advantage or burden in any endeavour to mathematically model dynamical systems since it entails the analytical tractability of the related modelling formalism and the efficiency of the numerical computation. The Reaction-Based Model (RBM) of the state space, which is presented in this article, is a novel formalization of the kinetics of a system of interacting molecules. According to our representation, the state S_μ of a system of M reactions and N molecular species, is identified with the occurrence of the reaction R_μ ( μ = 1, ..., M). The transition between any two states S_μ and S_ν is modelled as a first-order reaction S_μ → S_ν and described by mass action-like equation for the partial time derivative of the variables P(S_μ, t) and P(S_ν, t), which denote the probabilities that the system lies in the two states, respectively. The rate equations for the state probabilities are coupled with those for the abundance of molecular species. Altogether, the rate equations along with the specification of the initial conditions define the Cauchy problem whose solution describes the time-evolution of the system. The RBM has been applied to a typical severely stiff biological case study. The numerical solutions of the system's dynamics turned out to be computationally more efficient and in agreement with the results of the stochastic and hybrid stochastic/deterministic simulation algorithms.

Increasing Neural Stem Cell Division Asymmetry and Quiescence Are Predicted to Contribute to the Age-Related Decline in Neurogenesis

Adult murine neural stem cells (NSCs) generate neurons in drastically declining numbers with age. How cellular dynamics sustain neurogenesis and how alterations with age may result in this decline are unresolved issues. We therefore clonally traced NSC lineages using confetti reporters in young and middle-aged adult mice. To understand the underlying mechanisms, we derived mathematical models that explain observed clonal cell type abundances. The best models consistently show self-renewal of transit-amplifying progenitors and rapid neuroblast cell cycle exit. In middle-aged mice, we identified an increased probability of asymmetric stem cell divisions at the expense of symmetric differentiation, accompanied by an extended persistence of quiescence between activation phases. Our model explains existing longitudinal population data and identifies particular cellular properties underlying adult NSC homeostasis and the aging of this stem cell compartment.

SimKinet: A free educational tool based on an electrical analogy to solve chemical kinetic equations

In this article we introduce the software SimKinet, a free tool specifically designed to solve systems of differential equations without any programming skill. The underlying method is the so-called Network Simulation Method, which designs and solves an electrical network equivalent to the mathematical problem. SimKinet is versatile, fast, presenting a real user-friendly interface, and can be employed for both educational and researching purposes. It is particularly useful in the first courses of different scientific degrees, mainly Chemistry and Physics, especially when facing non-analytic or complex-dynamics problems. Moreover, SimKinet would help students to understand fundamental concepts, being an opportunity to improve instruction in Chemistry, Mathematics, Physics and other Sciences courses, with no need of advanced knowledge in differential equations. The potency of SimKinet is demonstrated via two applications in chemical kinetics: the photochemical destruction of stratospheric ozone and the chaotic dynamics of the peroxidase-oxidase reaction.

Chemical Kinetics Roots and Methods to Obtain the Probability Distribution Function Evolution of Reactants and Products in Chemical Networks Governed by a Master Equation

In this paper first, we review the physical root bases of chemical reaction networks as a Markov process in multidimensional vector space. Then we study the chemical reactions from a microscopic point of view, to obtain the expression for the propensities for the different reactions that can happen in the network. These chemical propensities, at a given time, depend on the system state at that time, and do not depend on the state at an earlier time indicating that we are dealing with Markov processes. Then the Chemical Master Equation (CME) is deduced for an arbitrary chemical network from a probability balance and it is expressed in terms of the reaction propensities. This CME governs the dynamics of the chemical system. Due to the difficulty to solve this equation two methods are studied, the first one is the probability generating function method or z-transform, which permits to obtain the evolution of the factorial moment of the system with time in an easiest way or after some manipulation the evolution of the polynomial moments. The second method studied is the expansion of the CME in terms of an order parameter (system volume). In this case we study first the expansion of the CME using the propensities obtained previously and splitting the molecular concentration into a deterministic part and a random part. An expression in terms of multinomial coefficients is obtained for the evolution of the probability of the random part. Then we study how to reconstruct the probability distribution from the moments using the maximum entropy principle. Finally, the previous methods are applied to simple chemical networks and the consistency of these methods is studied.

SELANSI: a toolbox for simulation of stochastic gene regulatory networks

Motivation

Gene regulation is inherently stochastic. In many applications concerning Systems and Synthetic Biology such as the reverse engineering and the de novo design of genetic circuits, stochastic effects (yet potentially crucial) are often neglected due to the high computational cost of stochastic simulations. With advances in these fields there is an increasing need of tools providing accurate approximations of the stochastic dynamics of gene regulatory networks (GRNs) with reduced computational effort.

Results

This work presents SELANSI (SEmi-LAgrangian SImulation of GRNs), a software toolbox for the simulation of stochastic multidimensional gene regulatory networks. SELANSI exploits intrinsic structural properties of gene regulatory networks to accurately approximate the corresponding Chemical Master Equation with a partial integral differential equation that is solved by a semi-lagrangian method with high efficiency. Networks under consideration might involve multiple genes with self and cross regulations, in which genes can be regulated by different transcription factors. Moreover, the validity of the method is not restricted to a particular type of kinetics. The tool offers total flexibility regarding network topology, kinetics and parameterization, as well as simulation options.

Availability and implementation

SELANSI runs under the MATLAB environment, and is available under GPLv3 license at https://sites.google.com/view/selansi.

A scalable moment-closure approximation for large-scale biochemical reaction networks

Motivation

Stochastic 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.

Results

In 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 implementation

The sMA is implemented in the open-source MATLAB toolbox CERENA and is available from https://github.com/CERENADevelopers/CERENA .

Contact

jan.hasenauer@helmholtz-muenchen.de or atefeh.kazeroonian@tum.de.

Supplementary information

Supplementary data are available at Bioinformatics online.

Robust parameter estimation for dynamical systems from outlier-corrupted data

Motivation

Dynamics of cellular processes are often studied using mechanistic mathematical models. These models possess unknown parameters which are generally estimated from experimental data assuming normally distributed measurement noise. Outlier corruption of datasets often cannot be avoided. These outliers may distort the parameter estimates, resulting in incorrect model predictions. Robust parameter estimation methods are required which provide reliable parameter estimates in the presence of outliers.

Results

In this manuscript, we propose and evaluate methods for estimating the parameters of ordinary differential equation models from outlier-corrupted data. As alternatives to the normal distribution as noise distribution, we consider the Laplace, the Huber, the Cauchy and the Student's t distribution. We assess accuracy, robustness and computational efficiency of estimators using these different distribution assumptions. To this end, we consider artificial data of a conversion process, as well as published experimental data for Epo-induced JAK/STAT signaling. We study how well the methods can compensate and discover artificially introduced outliers. Our evaluation reveals that using alternative distributions improves the robustness of parameter estimates.

Availability and Implementation

The MATLAB implementation of the likelihood functions using the distribution assumptions is available at Bioinformatics online.

Contact

jan.hasenauer@helmholtz-muenchen.de.

Supplementary information

Supplementary material are available at Bioinformatics online.

Generalized method of moments for estimating parameters of stochastic reaction networks

BACKGROUND

Discrete-state stochastic models have become a well-established approach to describe biochemical reaction networks that are influenced by the inherent randomness of cellular events. In the last years several methods for accurately approximating the statistical moments of such models have become very popular since they allow an efficient analysis of complex networks.

RESULTS

We propose a generalized method of moments approach for inferring the parameters of reaction networks based on a sophisticated matching of the statistical moments of the corresponding stochastic model and the sample moments of population snapshot data. The proposed parameter estimation method exploits recently developed moment-based approximations and provides estimators with desirable statistical properties when a large number of samples is available. We demonstrate the usefulness and efficiency of the inference method on two case studies.

CONCLUSIONS

The generalized method of moments provides accurate and fast estimations of unknown parameters of reaction networks. The accuracy increases when also moments of order higher than two are considered. In addition, the variance of the estimator decreases, when more samples are given or when higher order moments are included.

Inference for Stochastic Chemical Kinetics Using Moment Equations and System Size Expansion

Quantitative mechanistic models are valuable tools for disentangling biochemical pathways and for achieving a comprehensive understanding of biological systems. However, to be quantitative the parameters of these models have to be estimated from experimental data. In the presence of significant stochastic fluctuations this is a challenging task as stochastic simulations are usually too time-consuming and a macroscopic description using reaction rate equations (RREs) is no longer accurate. In this manuscript, we therefore consider moment-closure approximation (MA) and the system size expansion (SSE), which approximate the statistical moments of stochastic processes and tend to be more precise than macroscopic descriptions. We introduce gradient-based parameter optimization methods and uncertainty analysis methods for MA and SSE. Efficiency and reliability of the methods are assessed using simulation examples as well as by an application to data for Epo-induced JAK/STAT signaling. The application revealed that even if merely population-average data are available, MA and SSE improve parameter identifiability in comparison to RRE. Furthermore, the simulation examples revealed that the resulting estimates are more reliable for an intermediate volume regime. In this regime the estimation error is reduced and we propose methods to determine the regime boundaries. These results illustrate that inference using MA and SSE is feasible and possesses a high sensitivity.

In this manuscript, we introduce efficient methods for parameter estimation for stochastic processes. The stochasticity of chemical reactions can influence the average behavior of the considered system. For some biological systems, a microscopic, stochastic description is computationally intractable but a macroscopic, deterministic description too inaccurate. This inaccuracy manifests itself in an error in parameter estimates, which impede the predictive power of the proposed model. Until now, no rigorous analysis on the magnitude of the estimation error exists. We show by means of two simulation examples that using mesoscopic descriptions based on the system size expansions and moment-closure approximations can reduce this estimation error compared to inference using a macroscopic description. This reduction is most pronounced in an intermediate volume regime where the influence of stochasticity on the average behavior is moderately strong. For the JAK/STAT pathway where experimental data is available, we show that one parameter that was not structurally identifiable when using a macroscopic description becomes structurally identifiable when using a mesoscopic description for parameter estimation.