Stochastic Chemical Kinetics Research Articles

A differentiable Gillespie algorithm for simulating chemical kinetics, parameter estimation, and designing synthetic biological circuits.

The Gillespie algorithm is commonly used to simulate and analyze complex chemical reaction networks. Here, we leverage recent breakthroughs in deep learning to develop a fully differentiable variant of the Gillespie algorithm. The differentiable Gillespie algorithm (DGA) approximates dis-continuous operations in the exact Gillespie algorithm using smooth functions, allowing for the calculation of gradients using backpropagation. The DGA can be used to quickly and accurately learn kinetic parameters using gradient descent and design biochemical networks with desired properties. As an illustration, we apply the DGA to study stochastic models of gene promoters. We show that the DGA can be used to: (i) successfully learn kinetic parameters from experimental measurements of mRNA expression levels from two distinct E. coli promoters and (ii) design nonequilibrium promoter architectures with desired input-output relationships. These examples illustrate the utility of the DGA for analyzing stochastic chemical kinetics, including a wide variety of problems of interest to synthetic and systems biology.

(Invited) Local pH in the (Photo/Electro) Catalytic Environment

(Photo)electrodes designed for improved efficiency and selectivity for CO2 reduction are often nanostructured, for example highly roughened metal films, nanowire arrays, and catalyst particles within a gas diffusion electrode. In all cases, local pockets of water with dissolved CO2 and ions are necessary for the reduction reaction. Although these regions are thought of as having a pH, they are often too small to have a well-defined pH or be at thermodynamic equilibrium until a minimum electrolyte pool size is reached. Computational studies using stochastic chemical kinetics provide a window into this environment, bridging the nanometer scale to the bulk-like scale. The calculations show that the compositions of these regions are characterized by large concentration fluctuations including periods of time when some ions are not present at all. The calculations enable useful trends in electrolyte composition to be identified as a function of electrolyte pool size, providing insights to the complex environments surrounding the catalytic center.

Stochastic chemical kinetics of cell fate decision systems: From single cells to populations and back.

Stochastic chemical kinetics is a widely used formalism for studying stochasticity of chemical reactions inside single cells. Experimental studies of reaction networks are generally performed with cells that are part of a growing population, yet the population context is rarely taken into account when models are developed. Models that neglect the population context lose their validity whenever the studied system influences traits of cells that can be selected in the population, a property that naturally arises in the complex interplay between single-cell and population dynamics of cell fate decision systems. Here, we represent such systems as absorbing continuous-time Markov chains. We show that conditioning on non-absorption allows one to derive a modified master equation that tracks the time evolution of the expected population composition within a growing population. This allows us to derive consistent population dynamics models from a specification of the single-cell process. We use this approach to classify cell fate decision systems into two types that lead to different characteristic phases in emerging population dynamics. Subsequently, we deploy the gained insights to experimentally study a recurrent problem in biology: how to link plasmid copy number fluctuations and plasmid loss events inside single cells to growth of cell populations in dynamically changing environments.

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.

Optimal control of bioproduction in the presence of population heterogeneity.

Cell-to-cell variability, born of stochastic chemical kinetics, persists even in large isogenic populations. In the study of single-cell dynamics this is typically accounted for. However, on the population level this source of heterogeneity is often sidelined to avoid the inevitable complexity it introduces. The homogeneous models used instead are more tractable but risk disagreeing with their heterogeneous counterparts and may thus lead to severely suboptimal control of bioproduction. In this work, we introduce a comprehensive mathematical framework for solving bioproduction optimal control problems in the presence of heterogeneity. We study population-level models in which such heterogeneity is retained, and propose order-reduction approximation techniques. The reduced-order models take forms typical of homogeneous bioproduction models, making them a useful benchmark by which to study the importance of heterogeneity. Moreover, the derivation from the heterogeneous setting sheds light on parameter selection in ways a direct homogeneous outlook cannot, and reveals the source of approximation error. With view to optimally controlling bioproduction in microbial communities, we ask the question: when does optimising the reduced-order models produce strategies that work well in the presence of population heterogeneity? We show that, in some cases, homogeneous approximations provide remarkably accurate surrogate models. Nevertheless, we also demonstrate that this is not uniformly true: overlooking the heterogeneity can lead to significantly suboptimal control strategies. In these cases, the heterogeneous tools and perspective are crucial to optimise bioproduction.

A Compact and Power-Efficient Noise Generator for Stochastic Simulations.

This paper describes an adaptive noise generator circuit suitable for on-chip simulations of stochastic chemical kinetics. The circuit uses amplified BJT white noise and adaptive low-pass filtering to emulate the power spectrum and autocorrelation of random telegraph signals (RTS) with Poisson-distributed level transitions. A current-mode implementation in the AMS 0.35 μm BiCMOS process shows excellent agreement with theoretical results from the Gillespie stochastic simulation algorithm over a 60 dB range in mean current levels (modeling molecule count numbers). The circuit has an estimated layout area of 0.032 mm2 and typically consumes 400 μA, which are 73% and 50% less, respectively, than prior implementations. Moreover, it does not require any off-chip capacitors. Experimental results from a discrete board-level implementation of the circuit are in good agreement with theoretical predictions.

Estimation of monkeypox spread in a nonendemic country considering contact tracing and self-reporting: A stochastic modeling study.

In May 2022, monkeypox started to spread in nonendemic countries. To investigate contact tracing and self-reporting of the primary case in the local community, a stochastic model is developed. An algorithm based on Gillespie's stochastic chemical kinetics is used to quantify the number of infections, contacts, and duration from the arrival of the primary case to the detection of the index case (or until there are no more local infections). Different scenarios were set considering the delay in contact tracing and behavior of infectors. We found that the self-reporting behavior of a primary case is the most significant factor affecting outbreak size and duration. Scenarios with a self-reporting primary case have an 86% reduction in infections (average: 5-7, in a population of 10 000) and contacts (average: 27-72) compared with scenarios with a non-self-reporting primary case (average number of infections and contacts: 27-72 and 197-537, respectively). Doubling the number of close contacts per day is less impactful compared with the self-reporting behavior of the primary case as it could only increase the number of infections by 45%. Our study emphasizes the importance of the prompt detection of the primary case.

Stochastic chemical kinetics of cell fate decision systems: from single cells to populations and back

Stochastic chemical kinetics of cell fate decision systems: from single cells to populations and back

Stochastic chemical kinetics of cell fate decision systems: from single cells to populations and back

Stochastic chemical kinetics of cell fate decision systems: from single cells to populations and back

(Invited) Kinetics of Dye-Sensitized Water Oxidation Catalysis in Natural Sunlight

The couplings of ultrafast electronic excitation, charge injection and charge transfer processes to the catalytic cycles involved in water oxidation using molecular dye-catalyst photoelectrosynthesis systems are examined using detailed stochastic chemical kinetics simulations. By modeling the full timescale from subpicoseconds to many minutes and assuming illumination by the sun (AM 1.5), we have gained insights to how 4-photon catalytic cycle works when stimulated by relatively rare interactions with light. In this talk the model will be described, and its predictions for how catalyst and dye state populations evolve over time and how the flow of the catalytic cycle works will be discussed.

Revisiting moment-closure methods with heterogeneous multiscale population models

Stochastic chemical kinetics at the single-cell level give rise to heterogeneous populations of cells even when all individuals are genetically identical. This heterogeneity can lead to nonuniform behaviour within populations, including different growth characteristics, cell-fate dynamics, and response to stimuli. Ultimately, these diverse behaviours lead to intricate population dynamics that are inherently multiscale: the population composition evolves based on population-level processes that interact with stochastically distributed single-cell states. Therefore, descriptions that account for this heterogeneity are essential to accurately model and control chemical processes. However, for real-world systems such models are computationally expensive to simulate, which can make optimisation problems, such as optimal control or parameter inference, prohibitively challenging. Here, we consider a class of multiscale population models that incorporate population-level mechanisms while remaining faithful to the underlying stochasticity at the single-cell level and the interplay between these two scales. To address the complexity, we study an order-reduction approximations based on the distribution moments. Since previous moment-closure work has focused on the single-cell kinetics, extending these techniques to populations models prompts us to revisit old observations as well as tackle new challenges. In this extended multiscale context, we encounter the previously established observation that the simplest closure techniques can lead to non-physical system trajectories. Despite their poor performance in some systems, we provide an example where these simple closures outperform more sophisticated closure methods in accurately, efficiently, and robustly solving the problem of optimal control of bioproduction in a microbial consortium model.

Topological data analysis of high dimensional probability landscapes of biochemical reaction networks using persistent homology

Biomolecular species such as proteins, DNAs, and RNA interact in cell to form reaction networks that regulate cellular functions. At small copy numbers, stochasticity arises from thermal fluctuations and plays important roles on cellular phenotypes. The underlying stochastic chemical kinetics (SCK) is governed by the discrete chemical master equation (dCME). Recent developments of the ACME algorithm enables the exact solution to dCME. With the constructions of exactly computed time-evolving and steady state probability landscapes, it is now possible to investigate the landscape properties of many stochastic networks, including where the probability basins are located, their significance, how they are connected, and where cycles of various dimensions appear.

Beyond the chemical master equation: Stochastic chemical kinetics coupled with auxiliary processes.

The chemical master equation and its continuum approximations are indispensable tools in the modeling of chemical reaction networks. These are routinely used to capture complex nonlinear phenomena such as multimodality as well as transient events such as first-passage times, that accurately characterise a plethora of biological and chemical processes. However, some mechanisms, such as heterogeneous cellular growth or phenotypic selection at the population level, cannot be represented by the master equation and thus have been tackled separately. In this work, we propose a unifying framework that augments the chemical master equation to capture such auxiliary dynamics, and we develop and analyse a numerical solver that accurately simulates the system dynamics. We showcase these contributions by casting a diverse array of examples from the literature within this framework and applying the solver to both match and extend previous studies. Analytical calculations performed for each example validate our numerical results and benchmark the solver implementation.

Uncertain chemical reaction equation

Reaction rate is a particularly important research object in chemical kinetics, and it is a measure of how fast a chemical reaction goes. In order to illustrate and clarify the evolution of concentration of a substance involved in the reaction, this paper derives an uncertain chemical reaction equation based on the theory of uncertain differential equation. By using the actual observations, one can estimate the parameters presented in the uncertain chemical reaction equation. As an application, the half-life of reaction is investigated. Finally, a paradox for stochastic chemical kinetics is given.

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.

Stochastic Limit-Cycle Oscillations of a Nonlinear System Under Random Perturbations

Dynamical systems with $\epsilon$ small random perturbations appear in both continuous mechanical motions and discrete stochastic chemical kinetics. The present work provides a detailed analysis of the central limit theorem (CLT), with a time-inhomogeneous Gaussian process, near a deterministic limit cycle in $\mathbb{R}^n$. Based on the theory of random perturbations of dynamical systems and the WKB approximation respectively, results are developed in parallel from both standpoints of stochastic trajectories and transition probability density and their relations are elucidated. We show rigorously the correspondence between the local Gaussian fluctuations and the curvature of the large deviation rate function near its infimum, connecting the CLT and the large deviation principle of diffusion processes. We study uniform asymptotic behavior of stochastic limit cycles through the interchange of limits of time $t\to\infty$ and $\epsilon\to 0$. Three further characterizations of stochastic limit cycle oscillators are obtained: (i) An approximation of the probability flux near the cycle; (ii) Two special features of the vector field for the cyclic motion; (iii) A local entropy balance equation along the cycle with clear physical meanings. Lastly and different from the standard treatment, the origin of the $\epsilon$ in the theory is justified by a novel scaling hypothesis via constructing a sequence of stochastic differential equations.

Sensitivity analysis of the reaction occurrence and recurrence times in steady-state biochemical networks

Continuous-time stationary Markov jump processes among discrete sites are encountered in disparate biochemical ambits. Sites and connecting dynamical events form a ‘network’ in which the sites are the available system’s states, and the events are site-to-site transitions, or even neutral processes in which the system does not change site but the event is however detectable. Examples include conformational transitions in single biomolecules, stochastic chemical kinetics in the space of the molecules copy numbers, and even macroscopic steady-state reactive mixtures if one adopts the viewpoint of a tagged molecule (or even of a molecular moiety) whose state may change when it is involved in a chemical reaction. Among the variety of dynamical descriptors, here we focus on the first occurrence times (starting from a given site) and on the recurrence times of an event of interest. We develop the sensitivity analysis for the lowest moments of the statistical distribution of such times with respect to the rate constants of the network. In particular, simple expressions and inequalities allow us to establish a direct relationship between selective variation of rate constants and effect on average times and variances. As illustrative cases we treat the substrate inhibition in enzymatic catalysis in which a tagged enzyme molecule jumps between three states, and the basic two-site model of stochastic gene expression in which the single gene switches between active and inactive forms.

Ultrafast Relaxations in Ruthenium Polypyridyl Chromophores Determined by Stochastic Kinetics Simulations.

Maximizing the efficiency of solar energy conversion using dye assemblies rests on understanding where the energy goes following absorption. Transient spectroscopies in solution are useful for this purpose, and the time-resolved data are usually analyzed with a sum of exponentials. This treatment assumes that dynamic events are well separated in time, and that the resulting exponential prefactors and phenomenological lifetimes are related directly to primary physical values. Such assumptions break down for coincident absorption, emission, and excited state relaxation that occur in transient absorption and photoluminescence of tris(2,2'-bipyridine)ruthenium(2+) derivatives, confounding the physical meaning of the reported lifetimes. In this work, we use inductive modeling and stochastic chemical kinetics to develop a detailed description of the primary ultrafast photophysics in transient spectroscopies of a series of Ru dyes, as an alternative to sums of exponential analysis. Commonly invoked three-level schemes involving absorption, intersystem crossing (ISC), and slow nonradiative relaxation and incoherent emission to the ground state cannot reproduce the experimentally measured spectra. The kinetics simulations reveal that ultrafast decay from the singlet excited state manifold to the ground state competes with ISC to the triplet excited state, whose efficiency was determined to be less than unity. The populations predicted by the simulations are used to estimate the magnitudes of transition dipoles for excited state excitations and evaluate the influence of specific ligands. The mechanistic framework and methodology presented here are entirely general, applicable to other dye classes, and can be extended to include charge injection by molecules bound to semiconductor surfaces.

Typical Stochastic Paths in the Transient Assembly of Fibrous Materials.

When chemically fueled, molecular self-assembly can sustain dynamic aggregates of polymeric fibers-hydrogels-with tunable properties. If the fuel supply is finite, the hydrogel is transient, as competing reactions switch molecular subunits between active and inactive states, drive fiber growth and collapse, and dissipate energy. Because the process is away from equilibrium, the structure and mechanical properties can reflect the history of preparation. As a result, the formation of these active materials is not readily susceptible to a statistical treatment in which the configuration and properties of the molecular building blocks specify the resulting material structure. Here, we illustrate a stochastic-thermodynamic and information-theoretic framework for this purpose and apply it to these self-annihilating materials. Among the possible paths, the framework variationally identifies those that are typical-loosely, the minimum number with the majority of the probability. We derive these paths from computer simulations of experimentally-informed stochastic chemical kinetics and a physical kinetics model for the growth of an active hydrogel. The model reproduces features observed by confocal microscopy, including the fiber length, lifetime, and abundance as well as the observation of fast fiber growth and stochastic fiber collapse. The typical mesoscopic paths we extract are less than 0.23% of those possible, but they accurately reproduce material properties such as mean fiber length.

Modeling Genetic Circuit Behavior in Transiently Transfected Mammalian Cells.

Binning cells by plasmid copy number is a common practice for analyzing transient transfection data. In many kinetic models of transfected cells, protein production rates are assumed to be proportional to plasmid copy number. The validity of this assumption in transiently transfected mammalian cells is not clear; models based on this assumption appear unable to reproduce experimental flow cytometry data robustly. We hypothesize that protein saturation at high plasmid copy number is a reason previous models break down and validate our hypothesis by comparing experimental data and a stochastic chemical kinetics model. The model demonstrates that there are multiple distinct physical mechanisms that can cause saturation. On the basis of these observations, we develop a novel minimal bin-dependent ODE model that assumes different parameters for protein production in cells with low versus high numbers of plasmids. Compared to a traditional Hill-function-based model, the bin-dependent model requires only one additional parameter, but fits flow cytometry input-output data for individual modules up to twice as accurately. By composing together models of individually fit modules, we use the bin-dependent model to predict the behavior of six cascades and three feed-forward circuits. The bin-dependent models are shown to provide more accurate predictions on average than corresponding (composed) Hill-function-based models and predictions of comparable accuracy to EQuIP, while still providing a minimal ODE-based model that should be easy to integrate as a subcomponent within larger differential equation circuit models. Our analysis also demonstrates that accounting for batch effects is important in developing accurate composed models.