Reaction-diffusion Master Equation Research Articles

Macroscopic limit for stochastic chemical reactions involving diffusion and spatial heterogeneity

To model bio-chemical reaction systems with diffusion one can either use stochastic, microscopic reaction–diffusion master equations or deterministic, macroscopic reaction–diffusion system. The connection between these two models is not only theoretically important but also plays an essential role in applications. This paper considers the macroscopic limits of the chemical reaction–diffusion master equation for first-order chemical reaction systems in highly heterogeneous environments. More precisely, the diffusion coefficients as well as the reaction rates are spatially inhomogeneous and the reaction rates may also be discontinuous. By carefully discretizing these heterogeneities within a reaction–diffusion master equation model, we show that in the limit we recover the macroscopic reaction–diffusion system with inhomogeneous diffusion and reaction rates.

Kinetics of spontaneous phase transitions from wüstite to magnetite in superparamagnetic core-shell nanocubes of iron oxides.

Iron oxide nanoparticles with a wüstite structure have been prepared by thermal decomposition. In air, they undergo a spontaneous transition into a thermodynamically more stable magnetite structure that grows from the surface. The thickness of this magnetite shell increases with time, thereby producing a series of core-shell nanoparticles. We investigated the kinetics of this phase transition in 23 nm nanocubes using time-resolved XRD, from which the fractions of individual phases were determined by the Rietveld refinement. This kinetics is described theoretically using three coupled reaction-diffusion master equations for the concentrations of oxygen, wüstite, and magnetite, in which both the diffusion of oxygen and its reaction with wüstite are thermally activated. The coefficients of these terms were adjusted so that the predictions of the model reproduce the XRD data at 298 K and 353 K, whereas the predictive capability of the model was assessed by comparing its predictions with measurements at 403 K.

Modeling and simulation of Lac-Operon using reaction-diffusion master equation on heterogeneous parallel platforms

The Lac-Operon is a pivotal molecular mechanism controlling gene expression in response to presence or absence of lactose. Mastering its dynamics aids cellular comprehension and gene circuit design. Employing the multi-particle diffusion reaction-diffusion master equation (MPD-RDME), we model gene expression kinetics. This couples with the Gillespie algorithm-driven Next sub-volume method for particle reactions and diffusion, enabling predictive time-course molecule concentration insights via RDME and Gillespie. Such insights inform experimentation and elucidate cellular operations. Particularly, inducible operon kinetics, like Lac-Operon, hold applications in recombinant protein production and gene therapy. Our work constructs a Lac-Operon gene expression model utilizing MPD-RDME and Next sub-volume on a Heterogeneous Parallel Platform (HPP). Implementation employs the combination of host and GPU, and host and homogeneous GPU pairs. Simulation speed, affected by external inducers, repressor count, and lattice spacing, is compared. Protein and mRNA count under various conditions were simulated on host, GPU, and GPU pairs. Results display 75 %–80 % speed enhancement on the HPP compared to host and GPU combination. This adaptable method extends to explore diverse gene circuits, unraveling molecular-level biological system behaviors.

4D dynamical whole-cell simulations of a growing minimal cell.

JCVI-syn3A is a genetically minimal bacterial cell with a single chromosome consisting of 493 genes that has retained few regulatory proteins or small RNAs. We recently developed a 3D fully dynamical whole-cell model that simulates the first 20 minutes of the cell cycle for this minimal cell before replication and growth. This model, based on experimental characterization of Syn3A and measurements from numerous biochemical experiments in related organisms, consisted of complete reaction networks for genetic information processing and metabolism totaling over 7200 chemical reactions. We now extend this 4D model to simulate the trajectory of a complete cell cycle while we dynamically update the shape of the membrane, chromosome configuration (including replication state), and polysome structures. To capture the broad range of length- and time-scales of the various dynamic processes, we periodically communicate cell state information among several coupled simulation techniques, including: reaction-diffusion and chemical master equation simulations for stochastic reactions, numerical integration of ODE's for metabolic reactions, and coarse-grained Brownian dynamics for the chromosome and membrane. Membrane shapes during division are inferred from experimental fluorescent images of dividing cells. The chromosome configurations are placed using physically realistic energy functions. From the model, we determine time-dependent concentrations, spatial distributions, and reaction fluxes that offer insight into the principles of life for this minimal cell, for example: prediction of doubling time matching the experimental value, realistic distributions of mRNA half-lives, and time-dependent balance of ATP costs and production for all reactions in the cell.

A Gaussian jump process formulation of the reaction-diffusion master equation enables faster exact stochastic simulations.

We propose a Gaussian jump process model on a regular Cartesian lattice for the diffusion part of the Reaction-Diffusion Master Equation (RDME). We derive the resulting Gaussian RDME (GRDME) formulation from analogy with a kernel-based discretization scheme for continuous diffusion processes and quantify the limits of its validity relative to the classic RDME. We then present an exact stochastic simulation algorithm for the GRDME, showing that the accuracies of GRDME and RDME are comparable, but exact simulations of the GRDME require only a fraction of the computational cost of exact RDME simulations. We analyze the origin of this speedup and its scaling with problem dimension. The benchmarks suggest that the GRDME is a particularly beneficial model for diffusion-dominated systems in three dimensional spaces, often occurring in systems biology and cell biology.

Invasiveness of a Growth-Migration System in a Two-dimensional Percolation cluster: A Stochastic Mathematical Approach.

We developed a framework based on the software Unstructured Reaction-Diffusion Master Equation (URDME) to address tumor cells' proliferation and migration in a heterogeneous space, herein a 2D percolation cluster. A mitogenic paracrine signaling pathway is utilized phenomenologically to reveal how cells cooperate with one another. We modeled the emerging Allee effect using low seeding density culture (LSDC) assays to fit the model parameters. A Finite time scaling (FTS) function has been formulated to quantitatively analyze invasiveness of a virtual Growth-Migration (GM) system in mimicking the cancer cell growth. Through such simulation, we analyzed the GM dynamics of virtual model in mimicking the growth of BT-474 cancer cell populations in vitro in a 2D percolation cluster and calculated the successful penetration rate (SPR). By analyzing the temporal trajectories of the SPR, we could determine the critical exponents of the critical SPR scaling relation. The SPR transition point ([Formula: see text]), which is a fundamentally different from a conventional percolation transition point, is found to be negatively correlated with the invasiveness of this cancer cell. The [Formula: see text] of the three variations of the virtual GM system distinctly designated by varying paracrine-regulated Allee (PAllee) model phenotypes is 0.3408, 0.3675, and 0.4454, respectively. FTS algorithm thereon may serve as an approach to quantify invasiveness of tumor cells. Through a phenomenological paracrine model, inter-cell cooperation and mutual mitogenic boosting are enabled to elicit the Allee effect in the GM systems. The rationale behind such computationally tunable virtual mechanism can be applied to other circumstances concerning emerging processes.

CellDynaMo-stochastic reaction-diffusion-dynamics model: Application to search-and-capture process of mitotic spindle assembly.

We introduce a Stochastic Reaction-Diffusion-Dynamics Model (SRDDM) for simulations of cellular mechanochemical processes with high spatial and temporal resolution. The SRDDM is mapped into the CellDynaMo package, which couples the spatially inhomogeneous reaction-diffusion master equation to account for biochemical reactions and molecular transport within the Langevin Dynamics (LD) framework to describe dynamic mechanical processes. This computational infrastructure allows the simulation of hours of molecular machine dynamics in reasonable wall-clock time. We apply SRDDM to test performance of the Search-and-Capture of mitotic spindle assembly by simulating, in three spatial dimensions, dynamic instability of elastic microtubules anchored in two centrosomes, movement and deformations of geometrically realistic centromeres with flexible kinetochores and chromosome arms. Furthermore, the SRDDM describes the mechanics and kinetics of Ndc80 linkers mediating transient attachments of microtubules to the chromosomal kinetochores. The rates of these attachments and detachments depend upon phosphorylation states of the Ndc80 linkers, which are regulated in the model by explicitly accounting for the reactions of Aurora A and B kinase enzymes undergoing restricted diffusion. We find that there is an optimal rate of microtubule-kinetochore detachments which maximizes the accuracy of the chromosome connections, that adding chromosome arms to kinetochores improve the accuracy by slowing down chromosome movements, that Aurora A and kinetochore deformations have a small positive effect on the attachment accuracy, and that thermal fluctuations of the microtubules increase the rates of kinetochore capture and also improve the accuracy of spindle assembly.

A probabilistic framework for particle-based reaction–diffusion dynamics using classical Fock space representations

The modeling and simulation of stochastic reaction–diffusion processes is a topic of steady interest that is approached with a wide range of methods. At the level of particle-resolved descriptions, where chemical reactions are coupled to the spatial diffusion of individual particles, there exist comprehensive numerical simulation schemes, while the corresponding mathematical formalization is relatively underdeveloped. The aim of this paper is to provide a framework to systematically formulate the probabilistic evolution equation, termed chemical diffusion master equation (CDME), that governs particle-based stochastic reaction–diffusion processes. To account for the non-conserved and unbounded particle number of this type of open systems, we employ a classical analogue of the quantum mechanical Fock space that contains the symmetrized probability densities of the many-particle configurations in space. Following field-theoretical ideas of second quantization, we introduce creation and annihilation operators that act on single-particle densities and provide natural representations of symmetrized probability densities as well as of reaction and diffusion operators. These operators allow us to consistently and systematically formulate the CDME for arbitrary reaction schemes. The resulting form of the CDME further serves as the foundation to derive more coarse-grained descriptions of reaction–diffusion dynamics. In this regard, we show that a discretization of the evolution equation by projection onto a Fock subspace generated by a finite set of single-particle densities leads to a generalized form of the well-known reaction–diffusion master equation, which supports non-local reactions between grid cells and which converges properly in the continuum limit.

A novel stochastic simulation approach enables exploration of mechanisms for regulating polarity site movement.

Cells polarize their movement or growth toward external directional cues in many different contexts. For example, budding yeast cells grow toward potential mating partners in response to pheromone gradients. Directed growth is controlled by polarity factors that assemble into clusters at the cell membrane. The clusters assemble, disassemble, and move between different regions of the membrane before eventually forming a stable polarity site directed toward the pheromone source. Pathways that regulate clustering have been identified but the molecular mechanisms that regulate cluster mobility are not well understood. To gain insight into the contribution of chemical noise to cluster behavior we simulated clustering using the reaction-diffusion master equation (RDME) framework to account for molecular-level fluctuations. RDME simulations are a computationally efficient approximation, but their results can diverge from the underlying microscopic dynamics. We implemented novel concentration-dependent rate constants that improved the accuracy of RDME-based simulations, allowing us to efficiently investigate how cluster dynamics might be regulated. Molecular noise was effective in relocating clusters when the clusters contained low numbers of limiting polarity factors, and when Cdc42, the central polarity regulator, exhibited short dwell times at the polarity site. Cluster stabilization occurred when abundances or binding rates were altered to either lengthen dwell times or increase the number of polarity molecules in the cluster. We validated key results using full 3D particle-based simulations. Understanding the mechanisms cells use to regulate the dynamics of polarity clusters should provide insights into how cells dynamically track external directional cues.

Three-dimensional stochastic simulation of chemoattractant-mediated excitability in cells.

During the last decade, a consensus has emerged that the stochastic triggering of an excitable system drives pseudopod formation and subsequent migration of amoeboid cells. The presence of chemoattractant stimuli alters the threshold for triggering this activity and can bias the direction of migration. Though noise plays an important role in these behaviors, mathematical models have typically ignored its origin and merely introduced it as an external signal into a series of reaction-diffusion equations. Here we consider a more realistic description based on a reaction-diffusion master equation formalism to implement these networks. In this scheme, noise arises naturally from a stochastic description of the various reaction and diffusion terms. Working on a three-dimensional geometry in which separate compartments are divided into a tetrahedral mesh, we implement a modular description of the system, consisting of G-protein coupled receptor signaling (GPCR), a local excitation-global inhibition mechanism (LEGI), and signal transduction excitable network (STEN). Our models implement detailed biochemical descriptions whenever this information is available, such as in the GPCR and G-protein interactions. In contrast, where the biochemical entities are less certain, such as the LEGI mechanism, we consider various possible schemes and highlight the differences between them. Our simulations show that even when the LEGI mechanism displays perfect adaptation in terms of the mean level of proteins, the variance shows a dose-dependence. This differs between the various models considered, suggesting a possible means for determining experimentally among the various potential networks. Overall, our simulations recreate temporal and spatial patterns observed experimentally in both wild-type and perturbed cells, providing further evidence for the excitable system paradigm. Moreover, because of the overall importance and ubiquity of the modules we consider, including GPCR signaling and adaptation, our results will be of interest beyond the field of directed migration.

Stochastic reaction and diffusion systems in molecular communications: Recent results and open problems

Chemical reactions and diffusion are two basic mechanisms governing the dynamics of molecules in a fluid. As such, they play a critical role in molecular communication for channel modeling, design of detection rules, implementation of molecular circuits for computation, and modeling interactions with external biochemical systems. For finite numbers of information-carrying molecules, stochastic models naturally arise with the simplest example given by the Wiener process, often known as Brownian motion. Nevertheless, the Wiener process fails to be accurate when external forces, friction, and chemical reactions are present. Recently, there have been several contributions that tailor molecular communication systems to these more challenging channel conditions. In this paper, we first overview a general family of stochastic models of reaction and diffusion systems, including both Langevin diffusion and the reaction-diffusion master equation. These models form a basis for the use of these models as molecular communication channels, from which modulation and detection schemes can be developed. We survey recent results on the design of these schemes, with a focus on a recently developed approach which is robust to a wide range of channel models, known as equilibrium signaling. We then turn to the implementation of these detection schemes and related parameter estimation problems via stochastic molecular circuits, based on stochastic chemical reaction networks. Finally, interactions between molecular communication systems and stochastic biological systems as well as open problems are discussed. Our overarching goal is to highlight how the consideration of general stochastic models of reaction and diffusion can be utilized in order to widen the application of molecular communications both within engineered systems, and also as motivation for advances in the mathematical characterization of these models.

Geant4-DNA Modeling of Water Radiolysis beyond the Microsecond: An On-Lattice Stochastic Approach.

In this work, we use the next sub-volume method (NSM) to investigate the possibility of using the compartment-based (“on-lattice”) model to simulate water radiolysis. We, first, start with a brief description of the reaction-diffusion master equation (RDME) in a spatially discretized simulation volume (“mesh”), which is divided into sub-volumes (or “voxels”). We then discuss the choice of voxel size and merging technique of a given mesh, along with the evolution of the system using the hierarchical algorithm for the RDME (“hRDME”). Since the compartment-based model cannot describe high concentration species of early radiation-induced spurs, we propose a combination of the particle-based step-by-step (“SBS”) Brownian dynamics model and the compartment-based model (“SBS-RDME model”) for the simulation. We, finally, use the particle-based SBS Brownian dynamics model of Geant4-DNA as a reference to test the model implementation through several benchmarks. We find that the compartment-based model can efficiently simulate the system with a large number of species and for longer timescales, beyond the microsecond, with a reasonable computing time. Our aim in developing this model is to study the production and evolution of reactive oxygen species generated under irradiation with different dose rate conditions, such as in FLASH and conventional radiotherapy.

A multiscale compartment-based model of stochastic gene regulatory networks using hitting-time analysis.

Spatial stochastic models of single cell kinetics are capable of capturing both fluctuations in molecular numbers and the spatial dependencies of the key steps of intracellular regulatory networks. The spatial stochastic model can be simulated both on a detailed microscopic level using particle tracking and on a mesoscopic level using the reaction–diffusion master equation. However, despite substantial progress on simulation efficiency for spatial models in the last years, the computational cost quickly becomes prohibitively expensive for tasks that require repeated simulation of thousands or millions of realizations of the model. This limits the use of spatial models in applications such as multicellular simulations, likelihood-free parameter inference, and robustness analysis. Further approximation of the spatial dynamics is needed to accelerate such computational engineering tasks. We here propose a multiscale model where a compartment-based model approximates a detailed spatial stochastic model. The compartment model is constructed via a first-exit time analysis on the spatial model, thus capturing critical spatial aspects of the fine-grained simulations, at a cost close to the simple well-mixed model. We apply the multiscale model to a canonical model of negative-feedback gene regulation, assess its accuracy over a range of parameters, and demonstrate that the approximation can yield substantial speedups for likelihood-free parameter inference.

Effects of different discretisations of the Laplacian upon stochastic simulations of reaction–diffusion systems on both static and growing domains

By discretising space into compartments and letting system dynamics be governed by the reaction–diffusion master equation, it is possible to derive and simulate a stochastic model of reaction and diffusion on an arbitrary domain. However, there are many implementation choices involved in this process, such as the choice of discretisation and method of derivation of the diffusive jump rates, and it is not clear a priori how these affect model predictions. To shed light on this issue, in this work we explore how a variety of discretisations and methods for derivation of the diffusive jump rates affect the outputs of stochastic simulations of reaction–diffusion models, in particular using Turing’s model of pattern formation as a key example. We consider both static and uniformly growing domains and demonstrate that, while only minor differences are observed for simple reaction–diffusion systems, there can be vast differences in model predictions for systems that include complicated reaction kinetics, such as Turing’s diffusion-driven instability model of pattern formation. Our work highlights that care must be taken in using the reaction–diffusion master equation framework to make predictions as to the dynamics of stochastic reaction–diffusion systems.

Stochastic Spatial Simulation of Genetic Information Processes in the Minimal Cell

JCVI-syn3A is a genetically minimal cell whose genome consists of 493 genes, of which 452 are protein-coding. The genetic information processing reactions in JCVI-syn3A were recently solved stochastically using a well-stirred assumption which predicts the times of DNA replication initiation and elongation to be 8 and 50 minutes on average respectively after averaging over 1,000 replicate cells. This model approximately doubles the number of proteins according to genome-wide proteomics data as well as the number of ribosomal components over the course of a cell cycle. We introduce spatial heterogeneity and diffusion by constructing a 3D spatial model of JCVI-syn3A using cryo-electron tomograms to infer placement of large cellular components such as membrane boundaries and ribosomes. The ensemble of DNA configurations for the circular chromosome were generated using a random polymer where self-intersections and intersections with ribosomes were forbidden. Membrane-bound degradosome complexes were randomly distributed across the interior of the cell membrane. We fixed positions of the DNA, degradosomes, and ribosomes while simulating diffusion of RNA and proteins with more than 3,000 reactions including transcription of all 452 protein-coding genes and translation and degradation of all 452 mRNA under a reaction-diffusion master equation formulation using our GPU-based software, Lattice Microbes. Introducing spatial heterogeneity and requiring mRNA to diffuse to ribosomes and degradosomes to be translated or decayed results in a distribution of mRNA half-lives in agreement with measurements in other gram-positive bacteria Mycoplasma gallisepticum and Bacillus subtilis.

Gibbs free energy change of a discrete chemical reaction event.

In modeling the interior of cells by simulating a reaction-diffusion master equation over a grid of compartments, one employs the assumption that the copy numbers of various chemical species are small, discrete quantities. We show that, in this case, textbook expressions for the change in Gibbs free energy accompanying a chemical reaction or diffusion between adjacent compartments are inaccurate. We derive exact expressions for these free energy changes for the case of discrete copy numbers and show how these expressions reduce to traditional expressions under a series of successive approximations leveraging the relative sizes of the stoichiometric coefficients and the copy numbers of the solutes and solvent. Numerical results are presented to corroborate the claim that if the copy numbers are treated as discrete quantities, then only these more accurate expressions lead to correct behavior. Thus, the newly derived expressions are critical for correctly computing entropy production in mesoscopic simulations based on the reaction-diffusion master equation formalism.

A Novel Stochastic Simulation Approach Enables Exploration of Mechanisms to Regulate Polarization Dynamics

Cells polarize in the direction of chemical gradients to orient growth or migration. When budding yeast cells are exposed to a source of mating pheromone, polarity factors assemble into highly dynamic clusters at the plasma membrane, sometimes but not always oriented up the pheromone gradient. The clusters assemble, disassemble, and move between different regions of the cell membrane before eventually forming a more concentrated cluster stably oriented towards the pheromone source. Pathways to generate and stabilize polarity factor clusters have been identified, but clustering dynamics are not well understood. Molecular fluctuations are likely to contribute to clustering dynamics. To gain insight into how they do so, we performed efficient stochastic simulations within the reaction-diffusion master equation framework. While less accurate than particle-based simulations, this approach has been improved by recently derived scale-dependent mesoscopic association rates. A remaining problem is that such mesoscopic rates perform poorly at the high local molecular concentrations typical of polarity systems. We therefore derived and tested novel mesoscopic rates that improve the accuracy of such simulations. We observed that although some existing polarity models failed to produce highly dynamic signaling clusters, a single additional reaction, proposed in another model, enabled dynamic clustering behavior. We investigate which reactions affect cluster dynamics, and validate our predictions using a more realistic particle-based model in a 3D geometry. Overall, our novel stochastic simulation approach allowed us to efficiently explore possible mechanisms by which cells may control polarization dynamics important for gradient detection.

A hybrid method for micro-mesoscopic stochastic simulation of reaction-diffusion systems

The present paper introduces a new micro–meso hybrid algorithm based on the Ghost Cell Method concept in which the microscopic subdomain is governed by the Reactive Multi-Particle Collision (RMPC) dynamics. The mesoscopic subdomain is modeled using the Reaction-Diffusion Master Equation (RDME). The RDME is solved by means of the Inhomogeneous Stochastic Simulation Algorithm. No hybrid algorithm has hitherto used the RMPC dynamics for modeling reactions and the trajectories of each individual particle. The RMPC is faster than other molecular based methods and has the advantage of conserving mass, energy and momentum in the collision and free streaming steps. The new algorithm is tested on three reaction-diffusion systems. In all the systems studied, very good agreement with the deterministic solutions of the corresponding differential equations is obtained. In addition, it has been shown that proper discretization of the computational domain results in significant speed-ups in comparison with the full RMPC algorithm.

Exposing Inter-process Information for Efficient PDES of Spatial Stochastic Systems on Multicores

We present a new approach for efficient process synchronization in parallel discrete event simulation on multicore computers. We aim specifically at simulation of spatially extended stochastic system models where time intervals between successive inter-process events are highly variable and without lower bounds: This includes models governed by the mesoscopic Reaction-Diffusion Master Equation (RDME). A central part of our approach is a mechanism for optimism control, in which each process disseminates accurate information about timestamps of its future outgoing interprocess events to its neighbours. This information gives each process a precise basis for deciding when to pause local processing to reduce the risk of expensive rollbacks caused by future “delayed” incoming events. We apply our approach to a natural parallelization of the Next Subvolume Method (NSM) for simulating systems obeying RDME. Since this natural parallelization does not expose accurate timestamps of future interprocess events, we restructure it to expose such information, resulting in a simulation algorithm called Refined Parallel NSM (Refined PNSM). We have implemented Refined PNSM in a parallel simulator for spatial extended Markovian processes. On 32 cores, it achieves an efficiency ranging between 43--95% for large models, and on average 37% for small models, compared to an efficient sequential simulation without any code for parallelization. It is shown that the gain of restructuring the naive parallelization into Refined PNSM more than outweighs its overhead. We also show that our resulting simulator is superior in performance to existing simulators on multicores for comparable models.

Integration of Reaction-Diffusion Master Equation and Brownian Dynamics Methodologies to Simulate a Minimal Cell

Simulating whole-cells at the genome scale involves the integration of models which operate at incompatible length, time, and concentration scales. Hybrid simulation schemes are necessary to bridge these scales. Progress is presented on a method combining the reaction-diffusion master equation (RDME) with Brownian dynamics (BD) into Lattice Microbes, a GPU-based RDME simulation code. RDME is used for protein diffusion, while BD is used for diffusion of particles larger than the subvolume size (ribosomes). The BD integration is performed on the GPU, allowing both simulations to progress in lockstep. Owing to the dense molecular crowding observed in living cells, accounting for excluded volume interactions is necessary to predict the spatial distribution and diffusion time scale of chemical species. Particles interact through a shared excluded volume field which influences the spatial transition rate of RDME particles as well as the friction coefficient of BD particles. Integration of the methods shows little impact on the time-to-solution compared to a pure RDME simulation. A proof of principle is presented using the particle densities and geometry of a 500-nm diameter minimal cell. This minimal cell, named JCVI-syn3A, with a 543-kbp genome and 493 genes, provides a versatile platform to study the basics of life. This computational methodology will become a major component of the goal to simulate JCVI-syn3A at a spatially resolved, stochastic level.