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.

Risk of Harm Associated With Using Rapid Sequence Induction Intubation and Positive Pressure Ventilation in Patients With Hemorrhagic Shock.

Based on limited published evidence, physiological principles, clinical experience, and expertise, the author group has developed a consensus statement on the potential for iatrogenic harm with rapid sequence induction (RSI) intubation and positive-pressure ventilation (PPV) on patients in hemorrhagic shock. "In hemorrhagic shock, or any low flow (central hypovolemic) state, it should be noted that RSI and PPV are likely to cause iatrogenic harm by decreasing cardiac output." The use of RSI and PPV leads to an increased burden of shock due to a decreased cardiac output (CO)2 which is one of the primary determinants of oxygen delivery (DO2). The diminishing DO2 creates a state of systemic hypoxia, the severity of which will determine the magnitude of the shock (shock dose) and a growing deficit of oxygen, referred to as oxygen debt. Rapid accumulation of critical levels of oxygen debt results in coagulopathy and organ dysfunction and failure. Spontaneous respiration induced negative intrathoracic pressure (ITP) provides the pressure differential driving venous return. PPV subsequently increases ITP and thus right atrial pressure. The loss in pressure differential directly decreases CO and DO2 with a resultant increase in systemic hypoxia. If RSI and PPV are deemed necessary, prior or parallel resuscitation with blood products is required to mitigate post intervention reduction of DO2 and the potential for inducing cardiac arrest in the critically shocked patient.

Hierarchical algorithm for the reaction-diffusion master equation.

We have developed an algorithm coupling mesoscopic simulations on different levels in a hierarchy of Cartesian meshes. Based on the multiscale nature of the chemical reactions, some molecules in the system will live on a fine-grained mesh, while others live on a coarse-grained mesh. By allowing molecules to transfer from the fine levels to the coarse levels when appropriate, we show that we can save up to three orders of magnitude of computational time compared to microscopic simulations or highly resolved mesoscopic simulations, without losing significant accuracy. We demonstrate this in several numerical examples with systems that cannot be accurately simulated with a coarse-grained mesoscopic model.

Control over single-cell distribution of G1 lengths by WNT governs pluripotency.

The link between single-cell variation and population-level fate choices lacks a mechanistic explanation despite extensive observations of gene expression and epigenetic variation among individual cells. Here, we found that single human embryonic stem cells (hESCs) have different and biased differentiation potentials toward either neuroectoderm or mesendoderm depending on their G1 lengths before the onset of differentiation. Single-cell variation in G1 length operates in a dynamic equilibrium that establishes a G1 length probability distribution for a population of hESCs and predicts differentiation outcome toward neuroectoderm or mesendoderm lineages. Although sister stem cells generally share G1 lengths, a variable proportion of cells have asymmetric G1 lengths, which maintains the population dispersion. Environmental Wingless-INT (WNT) levels can control the G1 length distribution, apparently as a means of priming the fate of hESC populations once they undergo differentiation. As a downstream mechanism, global 5-hydroxymethylcytosine levels are regulated by G1 length and thereby link G1 length to differentiation outcomes of hESCs. Overall, our findings suggest that intrapopulation heterogeneity in G1 length underlies the pluripotent differentiation potential of stem cell populations.

A 3D Multiscale Model to Explore the Role of EGFR Overexpression in Tumourigenesis.

The epidermal growth factor receptor (EGFR) signalling cascade is one of the main pathways that regulate the survival and division of mammalian cells. It is also one of the most altered transduction pathways in cancer. Acquired mutations in the EGFR/ERK pathway can cause the overexpression of EGFR on the surface of the cell, while others downregulate the inactivation of switched on intracellular proteins such as Ras and Raf. This upregulates the activity of ERK and promotes cell division. We develop a 3D multiscale model to explore the role of EGFR overexpression on tumour initiation. In this model, cells are described as individual objects that move, interact, divide, proliferate, and die by apoptosis. We use Brownian Dynamics to describe the extracellular and intracellular regulations of cells as well as the spatial and stochastic effects influencing them. The fate of each cell depends on the number of active transcription factors in the nucleus. We use numerical simulations to investigate the individual and combined effects of mutations on the intracellular regulation of individual cells. Next, we show that the distance between active receptors increase the level of EGFR/ERK signalling. We demonstrate the usefulness of the model by quantifying the impact of mutational alterations in the EGFR/ERK pathway on the growth rate of in silico tumours.

Mesoscopic-microscopic spatial stochastic simulation with automatic system partitioning.

The reaction-diffusion master equation (RDME) is a model that allows for efficient on-lattice simulation of spatially resolved stochastic chemical kinetics. Compared to off-lattice hard-sphere simulations with Brownian dynamics or Green's function reaction dynamics, the RDME can be orders of magnitude faster if the lattice spacing can be chosen coarse enough. However, strongly diffusion-controlled reactions mandate a very fine mesh resolution for acceptable accuracy. It is common that reactions in the same model differ in their degree of diffusion control and therefore require different degrees of mesh resolution. This renders mesoscopic simulation inefficient for systems with multiscale properties. Mesoscopic-microscopic hybrid methods address this problem by resolving the most challenging reactions with a microscale, off-lattice simulation. However, all methods to date require manual partitioning of a system, effectively limiting their usefulness as "black-box" simulation codes. In this paper, we propose a hybrid simulation algorithm with automatic system partitioning based on indirect a priori error estimates. We demonstrate the accuracy and efficiency of the method on models of diffusion-controlled networks in 3D.

The Role of Chromatin Density in Cell Population Heterogeneity during Stem Cell Differentiation.

We incorporate three-dimensional (3D) conformation of chromosome (Hi-C) and single-cell RNA sequencing data together with discrete stochastic simulation, to explore the role of chromatin reorganization in determining gene expression heterogeneity during development. While previous research has emphasized the importance of chromatin architecture on activation and suppression of certain regulatory genes and gene networks, our study demonstrates how chromatin remodeling can dictate gene expression distribution by folding into distinct topological domains. We hypothesize that the local DNA density during differentiation accentuate transcriptional bursting due to the crowding effect of chromatin. This phenomenon yields a heterogeneous cell population, thereby increasing the potential of differentiation of the stem cells.

Reaction rates for reaction-diffusion kinetics on unstructured meshes.

The reaction-diffusion master equation is a stochastic model often utilized in the study of biochemical reaction networks in living cells. It is applied when the spatial distribution of molecules is important to the dynamics of the system. A viable approach to resolve the complex geometry of cells accurately is to discretize space with an unstructured mesh. Diffusion is modeled as discrete jumps between nodes on the mesh, and the diffusion jump rates can be obtained through a discretization of the diffusion equation on the mesh. Reactions can occur when molecules occupy the same voxel. In this paper, we develop a method for computing accurate reaction rates between molecules occupying the same voxel in an unstructured mesh. For large voxels, these rates are known to be well approximated by the reaction rates derived by Collins and Kimball, but as the mesh is refined, no analytical expression for the rates exists. We reduce the problem of computing accurate reaction rates to a pure preprocessing step, depending only on the mesh and not on the model parameters, and we devise an efficient numerical scheme to estimate them to high accuracy. We show in several numerical examples that as we refine the mesh, the results obtained with the reaction-diffusion master equation approach those of a more fine-grained Smoluchowski particle-tracking model.

Stochastic Simulation Service: Bridging the Gap between the Computational Expert and the Biologist.

We present StochSS: Stochastic Simulation as a Service, an integrated development environment for modeling and simulation of both deterministic and discrete stochastic biochemical systems in up to three dimensions. An easy to use graphical user interface enables researchers to quickly develop and simulate a biological model on a desktop or laptop, which can then be expanded to incorporate increasing levels of complexity. StochSS features state-of-the-art simulation engines. As the demand for computational power increases, StochSS can seamlessly scale computing resources in the cloud. In addition, StochSS can be deployed as a multi-user software environment where collaborators share computational resources and exchange models via a public model repository. We demonstrate the capabilities and ease of use of StochSS with an example of model development and simulation at increasing levels of complexity.

A framework for discrete stochastic simulation on 3D moving boundary domains.

We have developed a method for modeling spatial stochastic biochemical reactions in complex, three-dimensional, and time-dependent domains using the reaction-diffusion master equation formalism. In particular, we look to address the fully coupled problems that arise in systems biology where the shape and mechanical properties of a cell are determined by the state of the biochemistry and vice versa. To validate our method and characterize the error involved, we compare our results for a carefully constructed test problem to those of a microscale implementation. We demonstrate the effectiveness of our method by simulating a model of polarization and shmoo formation during the mating of yeast. The method is generally applicable to problems in systems biology where biochemistry and mechanics are coupled, and spatial stochastic effects are critical.

Macromolecular Crowding Regulates the Gene Expression Profile by Limiting Diffusion.

We seek to elucidate the role of macromolecular crowding in transcription and translation. It is well known that stochasticity in gene expression can lead to differential gene expression and heterogeneity in a cell population. Recent experimental observations by Tan et al. have improved our understanding of the functional role of macromolecular crowding. It can be inferred from their observations that macromolecular crowding can lead to robustness in gene expression, resulting in a more homogeneous cell population. We introduce a spatial stochastic model to provide insight into this process. Our results show that macromolecular crowding reduces noise (as measured by the kurtosis of the mRNA distribution) in a cell population by limiting the diffusion of transcription factors (i.e. removing the unstable intermediate states), and that crowding by large molecules reduces noise more efficiently than crowding by small molecules. Finally, our simulation results provide evidence that the local variation in chromatin density as well as the total volume exclusion of the chromatin in the nucleus can induce a homogenous cell population.

Reaction rates for a generalized reaction-diffusion master equation.

It has been established that there is an inherent limit to the accuracy of the reaction-diffusion master equation. Specifically, there exists a fundamental lower bound on the mesh size, below which the accuracy deteriorates as the mesh is refined further. In this paper we extend the standard reaction-diffusion master equation to allow molecules occupying neighboring voxels to react, in contrast to the traditional approach, in which molecules react only when occupying the same voxel. We derive reaction rates, in two dimensions as well as three dimensions, to obtain an optimal match to the more fine-grained Smoluchowski model and show in two numerical examples that the extended algorithm is accurate for a wide range of mesh sizes, allowing us to simulate systems that are intractable with the standard reaction-diffusion master equation. In addition, we show that for mesh sizes above the fundamental lower limit of the standard algorithm, the generalized algorithm reduces to the standard algorithm. We derive a lower limit for the generalized algorithm which, in both two dimensions and three dimensions, is of the order of the reaction radius of a reacting pair of molecules.

Convergence of methods for coupling of microscopic and mesoscopic reaction-diffusion simulations.

In this paper, three multiscale methods for coupling of mesoscopic (compartment-based) and microscopic (molecular-based) stochastic reaction-diffusion simulations are investigated. Two of the three methods that will be discussed in detail have been previously reported in the literature; the two-regime method (TRM) and the compartment-placement method (CPM). The third method that is introduced and analysed in this paper is called the ghost cell method (GCM), since it works by constructing a "ghost cell" in which molecules can disappear and jump into the compartment-based simulation. Presented is a comparison of sources of error. The convergent properties of this error are studied as the time step Δt (for updating the molecular-based part of the model) approaches zero. It is found that the error behaviour depends on another fundamental computational parameter h, the compartment size in the mesoscopic part of the model. Two important limiting cases, which appear in applications, are considered: (i) Δt → 0 and h is fixed; (ii) Δt → 0 and h → 0 such that √Δt/h is fixed. The error for previously developed approaches (the TRM and CPM) converges to zero only in the limiting case (ii), but not in case (i). It is shown that the error of the GCM converges in the limiting case (i). Thus the GCM is superior to previous coupling techniques if the mesoscopic description is much coarser than the microscopic part of the model.

Reaction rates for mesoscopic reaction-diffusion kinetics.

The mesoscopic reaction-diffusion master equation (RDME) is a popular modeling framework frequently applied to stochastic reaction-diffusion kinetics in systems biology. The RDME is derived from assumptions about the underlying physical properties of the system, and it may produce unphysical results for models where those assumptions fail. In that case, other more comprehensive models are better suited, such as hard-sphere Brownian dynamics (BD). Although the RDME is a model in its own right, and not inferred from any specific microscale model, it proves useful to attempt to approximate a microscale model by a specific choice of mesoscopic reaction rates. In this paper we derive mesoscopic scale-dependent reaction rates by matching certain statistics of the RDME solution to statistics of the solution of a widely used microscopic BD model: the Smoluchowski model with a Robin boundary condition at the reaction radius of two molecules. We also establish fundamental limits on the range of mesh resolutions for which this approach yields accurate results and show both theoretically and in numerical examples that as we approach the lower fundamental limit, the mesoscopic dynamics approach the microscopic dynamics. We show that for mesh sizes below the fundamental lower limit, results are less accurate. Thus, the lower limit determines the mesh size for which we obtain the most accurate results.

Stochastic reaction-diffusion processes with embedded lower-dimensional structures.

Small copy numbers of many molecular species in biological cells require stochastic models of the chemical reactions between the molecules and their motion. Important reactions often take place on one-dimensional structures embedded in three dimensions with molecules migrating between the dimensions. Examples of polymer structures in cells are DNA, microtubules, and actin filaments. An algorithm for simulation of such systems is developed at a mesoscopic level of approximation. An arbitrarily shaped polymer is coupled to a background Cartesian mesh in three dimensions. The realization of the system is made with a stochastic simulation algorithm in the spirit of Gillespie. The method is applied to model problems for verification and two more detailed models of transcription factor interaction with the DNA.

Single molecule simulations in complex geometries with embedded dynamic one-dimensional structures.

Stochastic models of reaction-diffusion systems are important for the study of biochemical reaction networks where species are present in low copy numbers or if reactions are highly diffusion limited. In living cells many such systems include reactions and transport on one-dimensional structures, such as DNA and microtubules. The cytoskeleton is a dynamic structure where individual fibers move, grow, and shrink. In this paper we present a simulation algorithm that combines single molecule simulations in three-dimensional space with single molecule simulations on one-dimensional structures of arbitrary shape. Molecules diffuse and react with each other in space, they associate with and dissociate from one-dimensional structures as well as diffuse and react with each other on the one-dimensional structure. A general curve embedded in space can be approximated by a piecewise linear curve to arbitrary accuracy. The resulting algorithm is hence very flexible. Molecules bound to a curve can move by pure diffusion or via active transport, and the curve can move in space as well as grow and shrink. The flexibility and accuracy of the algorithm is demonstrated in five numerical examples.

Reaction-diffusion master equation in the microscopic limit.

Stochastic modeling of reaction-diffusion kinetics has emerged as a powerful theoretical tool in the study of biochemical reaction networks. Two frequently employed models are the particle-tracking Smoluchowski framework and the on-lattice reaction-diffusion master equation (RDME) framework. As the mesh size goes from coarse to fine, the RDME initially becomes more accurate. However, recent developments have shown that it will become increasingly inaccurate compared to the Smoluchowski model as the lattice spacing becomes very fine. Here we give a general and simple argument for why the RDME breaks down. Our analysis reveals a hard limit on the voxel size for which no local RDME can agree with the Smoluchowski model and lets us quantify this limit in two and three dimensions. In this light we review and discuss recent work in which the RDME has been modified in different ways in order to better agree with the microscale model for very small voxel sizes.