Constant-complexity stochastic simulation algorithm with optimal binning

Advanced methods for gene network identification and noise decomposition from single-cell data

Central to analyzing noisy gene expression systems is solving the Chemical Master Equation (CME), which characterizes the probability evolution of the reacting species' copy numbers. Solving CMEs for high-dimensional systems suffers from the curse of dimensionality. Here, we propose a computational method for improved scalability through a divide-and-conquer strategy that optimally decomposes the whole system into a leader system and several conditionally independent follower subsystems. The CME is solved by combining Monte Carlo estimation for the leader system with stochastic filtering procedures for the follower subsystems. We demonstrate this method with high-dimensional numerical examples and apply it to identify a yeast transcription system at the single-cell resolution, leveraging mRNA time-course experimental data. The identification results enable an accurate examination of the heterogeneity in rate parameters among isogenic cells. To validate this result, we develop a noise decomposition technique exploiting time-course data but requiring no supplementary components, e.g., dual-reporters.

GillesPy2: A Biochemical Modeling Framework for Simulation Driven Biological Discovery

Stochastic modeling has become an essential tool for studying biochemical reaction networks. There is a growing need for user-friendly and feature-complete software for model design and simulation. To address this need, we present GillesPy2, an open-source framework for building and simulating mathematical and biochemical models. GillesPy2, a major upgrade from the original GillesPy package, is now a stand-alone Python 3 package. GillesPy2 offers an intuitive interface for robust and reproducible model creation, facilitating rapid and iterative development. In addition to expediting the model creation process, GillesPy2 offers efficient algorithms to simulate stochastic, deterministic, and hybrid stochastic-deterministic models.

Modeling Protein Aggregation Kinetics: The Method of Second Stochasticization

The nucleation of protein aggregates and their growth are important in determining the structure of the cell's membraneless organelles as well as the pathogenesis of many diseases. The large number of molecular types of such aggregates along with the intrinsically stochastic nature of aggregation challenges our theoretical and computational abilities. Kinetic Monte Carlo simulation using the Gillespie algorithm is a powerful tool for modeling stochastic kinetics, but it is computationally demanding when a large number of diverse species is involved. To explore the mechanisms and statistics of aggregation more efficiently, we introduce a new approach to model stochastic aggregation kinetics which introduces noise into already statistically averaged equations obtained using mathematical moment closure schemes. Stochastic moment equations summarize succinctly the dynamics of the large diversity of species with different molecularity involved in aggregation but still take into account the stochastic fluctuations that accompany not only primary and secondary nucleation but also aggregate elongation, dissociation, and fragmentation. This method of "second stochasticization" works well where the fluctuations are modest in magnitude as is often encountered in vivo where the number of protein copies in some computations can be in the hundreds to thousands. Simulations using second stochasticization reveal a scaling law that correlates the size of the fluctuations in aggregate size and number with the total number of monomers. This scaling law is confirmed using experimental data. We believe second stochasticization schemes will prove valuable for bridging the gap between in vivo cell biology and detailed modeling. (The code is released on https://github.com/MYTLab/stoch-agg.).

Stochastic simulation algorithms for computational systems biology: Exact, approximate, and hybrid methods

Nowadays, mathematical modeling is playing a key role in many different research fields. In the context of system biology, mathematical models and their associated computer simulations constitute essential tools of investigation. Among the others, they provide a way to systematically analyze systems perturbations, develop hypotheses to guide the design of new experimental tests, and ultimately assess the suitability of specific molecules as novel therapeutic targets. To these purposes, stochastic simulation algorithms (SSAs) have been introduced for numerically simulating the time evolution of a well-stirred chemically reacting system by taking proper account of the randomness inherent in such a system. In this work, we review the main SSAs that have been introduced in the context of exact, approximate, and hybrid stochastic simulation. Specifically, we will introduce the direct method (DM), the first reaction method (FRM), the next reaction method (NRM) and the rejection-based SSA (RSSA) in the area of exact stochastic simulation. We will then present the τ-leaping method and the chemical Langevin method in the area of approximate stochastic simulation and an implementation of the hybrid RSSA (HRSSA) in the context of hybrid stochastic-deterministic simulation. Finally, we will consider the model of the sphingolipid metabolism to provide an example of application of SSA to computational system biology by exemplifying how different simulation strategies may unveil different insights into the investigated biological phenomenon. This article is categorized under: Models of Systems Properties and Processes > Mechanistic Models Analytical and Computational Methods > Computational Methods.

Hybrid approaches for multiple-species stochastic reaction-diffusion models

Reaction-diffusion models are used to describe systems in fields as diverse as physics, chemistry, ecology and biology. The fundamental quantities in such models are individual entities such as atoms and molecules, bacteria, cells or animals, which move and/or react in a stochastic manner. If the number of entities is large, accounting for each individual is inefficient, and often partial differential equation (PDE) models are used in which the stochastic behaviour of individuals is replaced by a description of the averaged, or mean behaviour of the system. In some situations the number of individuals is large in certain regions and small in others. In such cases, a stochastic model may be inefficient in one region, and a PDE model inaccurate in another. To overcome this problem, we develop a scheme which couples a stochastic reaction-diffusion system in one part of the domain with its mean field analogue, i.e. a discretised PDE model, in the other part of the domain. The interface in between the two domains occupies exactly one lattice site and is chosen such that the mean field description is still accurate there. In this way errors due to the flux between the domains are small. Our scheme can account for multiple dynamic interfaces separating multiple stochastic and deterministic domains, and the coupling between the domains conserves the total number of particles. The method preserves stochastic features such as extinction not observable in the mean field description, and is significantly faster to simulate on a computer than the pure stochastic model.