Stochastic simulation of chemical reactions with spatial resolution and single molecule detail

Transfer of polarity information via diffusion of Wnt ligands in C. elegans embryos

Different signaling mechanisms concur to ensure robust tissue patterning and cell fate instruction during animal development. Most of these mechanisms rely on signaling proteins that are produced, transported, and detected. The spatiotemporal dynamics of signaling molecules are largely unknown, yet they determine signal activity's spatial range and time frame. Here, we use the Caenorhabditis elegans embryo to study how Wnt ligands, an evolutionarily conserved family of signaling proteins, dynamically organize to establish cell polarity in a developing tissue. We identify how Wnt ligands, produced in the posterior half of the embryos, spread extracellularly to transmit information to distant target cells in the anterior half. With quantitative live imaging and fluorescence correlation spectroscopy, we show that Wnt ligands diffuse through the embryo over a timescale shorter than the cell cycle, in the intercellular space, and outside the tissue below the eggshell. We extracted diffusion coefficients of Wnt ligands and their receptor Frizzled and characterized their co-localization. Integrating our different measurements and observations in a simple computational framework, we show how fast diffusion in the embryo can polarize individual cells through a time integration of the arrival of the ligands at the target cells. The polarity established at the tissue level by a posterior Wnt source can be transferred to the cellular level. Our results support a diffusion-based long-range Wnt signaling, which is consistent with the dynamics of developing processes.

Cell4D: a general purpose spatial stochastic simulator for cellular pathways

BACKGROUND

With the generation of vast compendia of biological datasets, the challenge is how best to interpret 'omics data alongside biochemical and other small-scale experiments to gain meaningful biological insights. Key to this challenge are computational methods that enable domain-users to generate novel hypotheses that can be used to guide future experiments. Of particular interest are flexible modeling platforms, capable of simulating a diverse range of biological systems with low barriers of adoption to those with limited computational expertise.

RESULTS

We introduce Cell4D, a spatial-temporal modeling platform combining a robust simulation engine with integrated graphics visualization, a model design editor, and an underlying XML data model capable of capturing a variety of cellular functions. Cell4D provides an interactive visualization mode, allowing intuitive feedback on model behavior and exploration of novel hypotheses, together with a non-graphics mode, compatible with high performance cloud compute solutions, to facilitate generation of statistical data. To demonstrate the flexibility and effectiveness of Cell4D, we investigate the dynamics of CEACAM1 localization in T-cell activation. We confirm the importance of Ca2+ microdomains in activating calmodulin and highlight a key role of activated calmodulin on the surface expression of CEACAM1. We further show how lymphocyte-specific protein tyrosine kinase can help regulate this cell surface expression and exploit spatial modeling features of Cell4D to test the hypothesis that lipid rafts regulate clustering of CEACAM1 to promote trans-binding to neighbouring cells.

CONCLUSIONS

Through demonstrating its ability to test and generate hypotheses, Cell4D represents an effective tool to help integrate knowledge across diverse, large and small-scale datasets.

A probabilistic algorithm for optimising the steady-state diffusional flux into a partially absorbing body

Cells in an aqueous environment absorb diffusing nutrient molecules through nanoscale protein channels in their outer membranes. Assuming that there are constraints on the number of such channels a cell can produce, we ask the question: given a nondepleting source of nutrients, what is the optimal distribution of these channels over the cell surface? We coarse-grain this problem, phrasing it as a diffusion problem with position-dependent Robin boundary conditions on the surface. The aim is to maximize the steady-state total flux through the partially absorbing surface under an integral constraint on the local reactivities. We develop an algorithm to tackle this problem that uses the stored and processed results of a particle-based simulation with reflective boundary conditions to a posteriori estimate absorption flux at essentially negligible additional computational cost. We validate the algorithm against a few cases for which analytical or semi-analytical results are available. We apply it to two examples: a spherical cell in the presence of a point source and a spheroidal cell with an isotropic source at infinity. In the former case, there is a significant gain relative to the homogeneous case, while in the latter case the gain is only [Formula: see text].

Particle-based simulations reveal two positive feedback loops allow relocation and stabilization of the polarity site during yeast mating

Many cells adjust the direction of polarized growth or migration in response to external directional cues. The yeast Saccharomyces cerevisiae orient their cell fronts (also called polarity sites) up pheromone gradients in the course of mating. However, the initial polarity site is often not oriented towards the eventual mating partner, and cells relocate the polarity site in an indecisive manner before developing a stable orientation. During this reorientation phase, the polarity site displays erratic assembly-disassembly behavior and moves around the cell cortex. The mechanisms underlying this dynamic behavior remain poorly understood. Particle-based simulations of the core polarity circuit revealed that molecular-level fluctuations are unlikely to overcome the strong positive feedback required for polarization and generate relocating polarity sites. Surprisingly, inclusion of a second pathway that promotes polarity site orientation generated relocating polarity sites with properties similar to those observed experimentally. This pathway forms a second positive feedback loop involving the recruitment of receptors to the cell membrane and couples polarity establishment to gradient sensing. This second positive feedback loop also allows cells to stabilize their polarity site once the site is aligned with the pheromone gradient.

Close encounters of the sticky kind: Brownian motion at absorbing boundaries

Encounter-based models of diffusion provide a probabilistic framework for analyzing the effects of a partially absorbing reactive surface, in which the probability of absorption depends upon the amount of surface-particle contact time. In this paper we develop a class of encounter-based models that deal with absorption at sticky boundaries. Sticky boundaries occur in a diverse range of applications, including cell biology, colloidal physics, finance, and human crowd dynamics. They also naturally arise in active matter, where confined active particles tend to spontaneously accumulate at boundaries even in the absence of any particle-particle interactions. We begin by constructing a one-dimensional encounter-based model of sticky Brownian motion (BM), which is based on the zero-range limit of nonsticky BM with a short-range attractive potential well near the origin. In this limit, the boundary-contact time is given by the amount of time (occupation time) that the particle spends at the origin. We calculate the joint probability density or propagator for the particle position and the occupation time, and then identify an absorption event as the first time that the occupation time crosses a randomly generated threshold. We illustrate the theory by considering diffusion in a finite interval with a partially absorbing sticky boundary at one end. We show how various quantities, such as the mean first passage time (MFPT) for single-particle absorption and the relaxation to steady state at the multiparticle level, depend on moments of the random threshold distribution. Finally, we determine how sticky BM can be obtained by taking a particular diffusion limit of a sticky run-and-tumble particle (RTP).

Dissecting the general mechanisms of protein cage self-assembly by coarse-grained simulations

The development of artificial protein cages has recently gained massive attention due to their promising application prospect as novel delivery vehicles for therapeutics. These nanoparticles are formed through a process called self-assembly, in which individual subunits spontaneously arrange into highly ordered patterns via non-covalent but specific interactions. Therefore, the first step toward the design of novel engineered protein cages is to understand the general mechanisms of their self-assembling dynamics. Here we have developed a new computational method to tackle this problem. Our method is based on a coarse-grained model and a diffusion-reaction simulation algorithm. Using a tetrahedral cage as test model, we showed that self-assembly of protein cage requires of a seeding process in which specific configurations of kinetic intermediate states are identified. We further found that there is a critical concentration to trigger self-assembly of protein cages. This critical concentration allows that cages can only be successfully assembled under a persistently high concentration. Additionally, phase diagram of self-assembly has been constructed by systematically testing the model across a wide range of binding parameters. Finally, our simulations demonstrated the importance of protein's structural flexibility in regulating the dynamics of cage assembly. In summary, this study throws lights on the general principles underlying self-assembly of large cage-like protein complexes and thus provides insights to design new nanomaterials.

Spatial Stochastic Model of the Pre-B Cell Receptor

Survival and proliferation of immature B lymphocytes requires expression and tonic signaling of the pre-B cell receptor (pre-BCR). This low level, ligand-independent signaling is likely achieved through frequent, but short-lived, homo interactions. Tonic signaling is also central in the pathology of precursor B acute lymphoblastic leukemia (B-ALL). In order to understand how repeated, transient events can lead to sustained signaling and to assess the impact of receptor accumulation induced by the membrane landscape, we developed a spatial stochastic model of receptor aggregation and downstream signaling events. Our rule- and agent-based model builds on previous mature BCR signaling models and incorporates novel parameters derived from single particle tracking of pre-BCR on surfaces of two different B-ALL cell lines, 697 and Nalm6. Live cell tracking of receptors on the two cell lines revealed characteristic differences in their dimer dissociation rates and diffusion coefficients. We report here that these differences affect pre-BCR aggregation and consequent signal initiation events. Receptors on Nalm6 cells, which have a lower off-rate and lower diffusion coefficient, more frequently form higher order oligomers than pre-BCR on 697 cells, resulting in higher levels of downstream phosphorylation in the Nalm6 cell line.

Rapid Particle-Based Simulations of Cellular Signalling with the FLAME-Accelerated Signalling Tool (FaST) and GPUs

Cellular signalling is a vital process in living organisms for coordinating highly diverse responses to various stimuli. Particle-based modelling excels in its ability to model complex features of cellular signalling pathways including stochasticity, spatial effects, and heterogeneity, thus improving our understanding of critical decision processes in biology. Yet, particle-based modelling is computationally prohibitive to implement. We recently developed FaST (FLAME-accelerated signalling tool), a software tool that harnesses the power of high-performance computation to reduce the computational burden of particle-based modelling. In particular, employing the unique massively parallel architecture of graphic processing units (GPUs) provided extreme speed ups of simulations by >650-fold. In this chapter, we provide a step-by-step walkthrough of how to use FaST to create GPU-accelerated simulations of a simple cellular signalling network. We further explore how the flexibility of FaST can be used to implement entirely customized simulations while still including the intrinsic speed up advantages of GPU-based parallelization.

Close agreement between deterministic versus stochastic modeling of first-passage time to vesicle fusion

Ca2+-dependent cell processes, such as neurotransmitter or endocrine vesicle fusion, are inherently stochastic due to large fluctuations in Ca2+ channel gating, Ca2+ diffusion, and Ca2+ binding to buffers and target sensors. However, previous studies revealed closer-than-expected agreement between deterministic and stochastic simulations of Ca2+ diffusion, buffering, and sensing if Ca2+ channel gating is not Ca2+ dependent. To understand this result more fully, we present a comparative study complementing previous work, focusing on Ca2+ dynamics downstream of Ca2+ channel gating. Specifically, we compare deterministic (mean-field/mass-action) and stochastic simulations of vesicle exocytosis latency, quantified by the probability density of the first-passage time (FPT) to the Ca2+-bound state of a vesicle fusion sensor, following a brief Ca2+ current pulse. We show that under physiological constraints, the discrepancy between FPT densities obtained using the two approaches remains small even if as few as ∼50 Ca2+ ions enter per single channel-vesicle release unit. Using a reduced two-compartment model for ease of analysis, we illustrate how this close agreement arises from the smallness of correlations between fluctuations of the reactant molecule numbers, despite the large magnitude of fluctuation amplitudes. This holds if all relevant reactions are heteroreaction between molecules of different species, as is the case for bimolecular Ca2+ binding to buffers and downstream sensor targets. In this case, diffusion and buffering effectively decorrelate the state of the Ca2+ sensor from local Ca2+ fluctuations. Thus, fluctuations in the Ca2+ sensor's state underlying the FPT distribution are only weakly affected by the fluctuations in the local Ca2+ concentration around its average, deterministically computable value.

STEPS 4.0: Fast and memory-efficient molecular simulations of neurons at the nanoscale

Recent advances in computational neuroscience have demonstrated the usefulness and importance of stochastic, spatial reaction-diffusion simulations. However, ever increasing model complexity renders traditional serial solvers, as well as naive parallel implementations, inadequate. This paper introduces a new generation of the STochastic Engine for Pathway Simulation (STEPS) project (http://steps.sourceforge.net/), denominated STEPS 4.0, and its core components which have been designed for improved scalability, performance, and memory efficiency. STEPS 4.0 aims to enable novel scientific studies of macroscopic systems such as whole cells while capturing their nanoscale details. This class of models is out of reach for serial solvers due to the vast quantity of computation in such detailed models, and also out of reach for naive parallel solvers due to the large memory footprint. Based on a distributed mesh solution, we introduce a new parallel stochastic reaction-diffusion solver and a deterministic membrane potential solver in STEPS 4.0. The distributed mesh, together with improved data layout and algorithm designs, significantly reduces the memory footprint of parallel simulations in STEPS 4.0. This enables massively parallel simulations on modern HPC clusters and overcomes the limitations of the previous parallel STEPS implementation. Current and future improvements to the solver are not sustainable without following proper software engineering principles. For this reason, we also give an overview of how the STEPS codebase and the development environment have been updated to follow modern software development practices. We benchmark performance improvement and memory footprint on three published models with different complexities, from a simple spatial stochastic reaction-diffusion model, to a more complex one that is coupled to a deterministic membrane potential solver to simulate the calcium burst activity of a Purkinje neuron. Simulation results of these models suggest that the new solution dramatically reduces the per-core memory consumption by more than a factor of 30, while maintaining similar or better performance and scalability.

Stochastic model of ERK-mediated progesterone receptor translocation, clustering and transcriptional activity

Progesterone receptor (PR) transcriptional activity is a key factor in the differentiation of the uterine endometrium. By consequence, progestin has been identified as an important treatment modality for endometrial cancer. PR transcriptional activity is controlled by extracellular-signal-regulated kinase (ERK) mediated phosphorylation, downstream of growth factor receptors such as EGFR. However, phosphorylation of PR also targets it for ubiquitination and destruction in the proteasome. Quantitative studies of these opposing roles are much needed toward validation of potential new progestin-based therapeutics. In this work, we propose a spatial stochastic model to study the effects of the opposing roles for PR phosphorylation on the levels of active transcription factor. Our numerical simulations confirm earlier in vitro experiments in endometrial cancer cell lines, identifying clustering as a mechanism that amplifies the ability of progesterone receptors to influence gene transcription. We additionally show the usefulness of a statistical method we developed to quantify and control variations in stochastic simulations in general biochemical systems, assisting modelers in defining minimal but meaningful numbers of simulations while guaranteeing outputs remain within a pre-defined confidence level.

Detailed balance for particle models of reversible reactions in bounded domains

In particle-based stochastic reaction–diffusion models, reaction rates and placement kernels are used to decide the probability per time a reaction can occur between reactant particles and to decide where product particles should be placed. When choosing kernels to use in reversible reactions, a key constraint is to ensure that detailed balance of spatial reaction fluxes holds at all points at equilibrium. In this work, we formulate a general partial-integral differential equation model that encompasses several of the commonly used contact reactivity (e.g., Smoluchowski-Collins-Kimball) and volume reactivity (e.g., Doi) particle models. From these equations, we derive a detailed balance condition for the reversible A + B ⇆ C reaction. In bounded domains with no-flux boundary conditions, when choosing unbinding kernels consistent with several commonly used binding kernels, we show that preserving detailed balance of spatial reaction fluxes at all points requires spatially varying unbinding rate functions near the domain boundary. Brownian dynamics simulation algorithms can realize such varying rates through ignoring domain boundaries during unbinding and rejecting unbinding events that result in product particles being placed outside the domain.

Hybrid modelling of biological systems: current progress and future prospects

Integrated modelling of biological systems is becoming a necessity for constructing models containing the major biochemical processes of such systems in order to obtain a holistic understanding of their dynamics and to elucidate emergent behaviours. Hybrid modelling methods are crucial to achieve integrated modelling of biological systems. This paper reviews currently popular hybrid modelling methods, developed for systems biology, mainly revealing why they are proposed, how they are formed from single modelling formalisms and how to simulate them. By doing this, we identify future research requirements regarding hybrid approaches for further promoting integrated modelling of biological systems.

Understanding the functional role of membrane confinements in TNF-mediated signaling by multiscale simulations

The interaction between TNFα and TNFR1 is essential in maintaining tissue development and immune responses. While TNFR1 is a cell surface receptor, TNFα exists in both soluble and membrane-bound forms. Interestingly, it was found that the activation of TNFR1-mediated signaling pathways is preferentially through the soluble form of TNFα, which can also induce the clustering of TNFR1 on plasma membrane of living cells. We developed a multiscale simulation framework to compare receptor clustering induced by soluble and membrane-bound ligands. Comparing with the freely diffusive soluble ligands, we hypothesize that the conformational dynamics of membrane-bound ligands are restricted, which affects the clustering of ligand-receptor complexes at cell-cell interfaces. Our simulation revealed that only small clusters can form if TNFα is bound on cell surface. In contrast, the clustering triggered by soluble TNFα is more dynamic, and the size of clusters is statistically larger. We therefore demonstrated the impact of membrane-bound ligand on dynamics of receptor clustering. Moreover, considering that larger TNFα-TNFR1 clusters is more likely to provide spatial platform for downstream signaling pathway, our studies offer new mechanistic insights about why the activation of TNFR1-mediated signaling pathways is not preferred by membrane-bound form of TNFα.

Narrow capture problem: An encounter-based approach to partially reactive targets

A general topic of current interest is the analysis of diffusion problems in singularly perturbed domains with small interior targets or traps (the narrow capture problem). One major application is to intracellular diffusion, where the targets typically represent some form of reactive biochemical substrate. Most studies of the narrow capture problem treat the target boundaries as totally absorbing (Dirichlet), that is, the chemical reaction occurs immediately on first encounter between particle and target surface. In this paper, we analyze the three-dimensional narrow capture problem in the more realistic case of partially reactive target boundaries. We begin by considering classical Robin boundary conditions. Matching inner and outer solutions of the single-particle probability density, we derive an asymptotic expansion of the Laplace transformed flux into each reactive surface in powers of ε, where ερ is a given target size. In turn, the fluxes determine the splitting probabilities for target absorption. We then extend our analysis to more general types of reactive targets by combining matched asymptotic analysis with an encounter-based formulation of diffusion-mediated surface reactions. That is, we derive an asymptotic expansion of the joint probability density for particle position and the so-called boundary local time, which characterizes the amount of time that a Brownian particle spends in the neighborhood of a point on a totally reflecting boundary. The effects of surface reactions are then incorporated via an appropriate stopping condition for the boundary local time. Robin boundary conditions are recovered in the special case of an exponential law for the stopping local times. Finally, we illustrate the theory by exploring how the leading-order contributions to the splitting probabilities depend on the choice of surface reactions. In particular, we show that there is an effective renormalization of the target radius of the form ρ→ρ-Ψ[over ̃](1/ρ), where Ψ[over ̃] is the Laplace transform of the stopping local time distribution.

Simulation of receptor triggering by kinetic segregation shows role of oligomers and close-contacts

The activation of T cells, key players of the immune system, involves local evacuation of phosphatase CD45 from a region of the T cell's surface, segregating it from the T cell receptor. What drives this evacuation? In the presence of antigen, what ensures evacuation happens in the subsecond timescales necessary to initiate signaling? In the absence of antigen, what mechanisms ensure that evacuation does not happen spontaneously, which could cause signaling errors? Phenomena known to influence spatial organization of CD45 or similar surface molecules include diffusive motion in the lipid bilayer, oligomerization reactions, and mechanical compression against a nearby surface, such as that of the cell presenting the antigen. Computer simulations can investigate hypothesized spatiotemporal mechanisms of T cell signaling. The challenge to computational studies of evacuation is that the base process, spontaneous evacuation by simple diffusion, is in the extreme rare event limit, meaning direct stochastic simulation is unfeasible. Here, we combine particle-based spatial stochastic simulation with the weighted ensemble method for rare events to compute the mean first passage time for cell surface availability by surface reorganization of CD45. We confirm mathematical estimates that, at physiological concentrations, spontaneous evacuation is extremely rare, roughly 300 years. We find that dimerization decreases the time required for evacuation. A weak bimolecular interaction (dissociation constant estimate 460 μM) is sufficient for an order of magnitude reduction of spontaneous evacuation times, and oligomerization to hexamers reduces times to below 1 s. This introduces a mechanism whereby externally induced CD45 oligomerization could significantly modify T cell function. For large regions of close contact, such as those induced by large microvilli, molecular size and compressibility imply a nonzero reentry probability of 60%, decreasing evacuation times. Simulations show that these reduced evacuation times are still unrealistically long (even with a fourfold variation centered around previous estimates of parameters), suggesting that a yet-to-be-described mechanism, besides compressional exclusion at a close contact, drives evacuation.

Contact statistics in populations of noninteracting random walkers in two dimensions

The interaction between individuals in biological populations, dilute components of chemical systems, or particles transported by turbulent flows depends critically on their contact statistics. This work clarifies those statistics under the simplifying assumptions that the underlying motions approximate a Brownian random walk and that the particles can be treated as noninteracting. We measure the contact-interval (also called the waiting-time or interarrival-time), contact-count, and contact-duration distributions in populations of individuals undergoing noninteracting continuous-space-time random walks on a periodic two-dimensional plane (a torus) as functions of the population number density, walker radius, and random-walk step size. The contact-interval is exponentially distributed for times longer than the ballistic mean-free-collision time but not for times shorter than that, and the contact duration distribution is strongly peaked at the ballistic-crossing time for head-on collisions when the ballistic-crossing time is short compared to the mean step duration. While successive contacts between individuals are independent, the probability of repeat contact decreases with time after a previous contact. This leads to a negative duration dependence of the waiting-time interval and overdispersion of the contact-count probability density function for all time intervals. The paper demonstrates that for populations of small particles (with a walker radius that is small compared to the mean-separation or random-walk step size), the ballistic mean-free-collision interval, the ballistic-crossing time, and the random-walk-step duration can be used to construct temporal scalings which allow for common waiting-time, contact-count, and contact-duration distributions across different populations. Semi-analytic approximations for both the waiting-time and contact-duration distributions are also presented.

Multiscale molecular kinetics by coupling Markov state models and reaction-diffusion dynamics

A novel approach to simulate simple protein-ligand systems at large time and length scales is to couple Markov state models (MSMs) of molecular kinetics with particle-based reaction-diffusion (RD) simulations, MSM/RD. Currently, MSM/RD lacks a mathematical framework to derive coupling schemes, is limited to isotropic ligands in a single conformational state, and lacks multiparticle extensions. In this work, we address these needs by developing a general MSM/RD framework by coarse-graining molecular dynamics into hybrid switching diffusion processes. Given enough data to parameterize the model, it is capable of modeling protein-protein interactions over large time and length scales, and it can be extended to handle multiple molecules. We derive the MSM/RD framework, and we implement and verify it for two protein-protein benchmark systems and one multiparticle implementation to model the formation of pentameric ring molecules. To enable reproducibility, we have published our code in the MSM/RD software package.

Cell patterning by secretion-induced plasma membrane flows

Cells self-organize using reaction-diffusion and fluid-flow principles. Whether bulk membrane flows contribute to cell patterning has not been established. Here, using mathematical modeling, optogenetics, and synthetic probes, we show that polarized exocytosis causes lateral membrane flows away from regions of membrane insertion. Plasma membrane–associated proteins with sufficiently low diffusion and/or detachment rates couple to the flows and deplete from areas of exocytosis. In rod-shaped fission yeast cells, zones of Cdc42 GTPase activity driving polarized exocytosis are limited by GTPase activating proteins (GAPs). We show that membrane flows pattern the GAP Rga4 distribution and that coupling of a synthetic GAP to membrane flows is sufficient to establish the rod shape. Thus, membrane flows induced by Cdc42-dependent exocytosis form a negative feedback restricting the zone of Cdc42 activity.

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.

Cells localize polarity molecules in a small region of the plasma membrane forming a polarity cluster that directs functions such as migration, reproduction, and growth. Guided by external signals, these clusters move across the membrane allowing cells to reorient growth or motion. The polarity molecules continuously and randomly shuttle between the cluster and the cell cytosol and, as a result, the number and distribution of molecules at the cluster constantly changes. Here we present an improved stochastic simulation algorithm to investigate how such molecular-scale fluctuations induce cluster movement across the cell membrane. Unexpectedly, cluster mobility does not correlate with variations in total molecule abundance within the cluster, but rather with changes in the spatial distribution of molecules that form the cluster. Cluster motion is faster when polarity molecules are scarce and when they shuttle rapidly between the cluster and the cytosol. Our results suggest that cells control cluster mobility by regulating the abundance of polarity molecules and biochemical reactions that affect the time molecules spend at the cluster. We provide insights into how cells harness random molecular behavior to perform functions important for survival, such as detecting the direction of external signals.

Diffusive search and trajectories on tubular networks: a propagator approach

Several organelles in eukaryotic cells, including mitochondria and the endoplasmic reticulum, form interconnected tubule networks extending throughout the cell. These tubular networks host many biochemical pathways that rely on proteins diffusively searching through the network to encounter binding partners or localized target regions. Predicting the behavior of such pathways requires a quantitative understanding of how confinement to a reticulated structure modulates reaction kinetics. In this work, we develop both exact analytical methods to compute mean first passage times and efficient kinetic Monte Carlo algorithms to simulate trajectories of particles diffusing in a tubular network. Our approach leverages exact propagator functions for the distribution of transition times between network nodes and allows large simulation time steps determined by the network structure. The methodology is applied to both synthetic planar networks and organelle network structures, demonstrating key general features such as the heterogeneity of search times in different network regions and the functional advantage of broadly distributing target sites throughout the network. The proposed algorithms pave the way for future exploration of the interrelationship between tubular network structure and biomolecular reaction kinetics.

Mechanistic dissection of spatial organization in NF-κB signaling pathways by hybrid simulations

The nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) is one of the most important transcription factors involved in the regulation of inflammatory signaling pathways. Inappropriate activation of these pathways has been linked to autoimmunity and cancers. Emerging experimental evidences have been showing the existence of elaborate spatial organizations for various molecular components in the pathways. One example is the scaffold protein tumor necrosis factor receptor associated factor (TRAF). While most TRAF proteins form trimeric quaternary structure through their coiled-coil regions, the N-terminal region of some members in the family can further be dimerized. This dimerization of TRAF trimers can drive them into higher-order clusters as a response to receptor stimulation, which functions as a spatial platform to mediate the downstream poly-ubiquitination. However, the molecular mechanism underlying the TRAF protein clustering and its functional impacts are not well-understood. In this article, we developed a hybrid simulation method to tackle this problem. The assembly of TRAF-based signaling platform at the membrane-proximal region is modeled with spatial resolution, while the dynamics of downstream signaling network, including the negative feedbacks through various signaling inhibitors, is simulated as stochastic chemical reactions. These two algorithms are further synchronized under a multiscale simulation framework. Using this computational model, we illustrated that the formation of TRAF signaling platform can trigger an oscillatory NF-κB response. We further demonstrated that the temporal patterns of downstream signal oscillations are closely regulated by the spatial factors of TRAF clustering, such as the geometry and energy of dimerization between TRAF trimers. In general, our study sheds light on the basic mechanism of NF-κB signaling pathway and highlights the functional importance of spatial regulation within the pathway. The simulation framework also showcases its potential of application to other signaling pathways in cells.

A multiscale study on the mechanisms of spatial organization in ligand-receptor interactions on cell surfaces

The binding of cell surface receptors with extracellular ligands triggers distinctive signaling pathways, leading into the corresponding phenotypic variation of cells. It has been found that in many systems, these ligand-receptor complexes can further oligomerize into higher-order structures. This ligand-induced oligomerization of receptors on cell surfaces plays an important role in regulating the functions of cell signaling. The underlying mechanism, however, is not well understood. One typical example is proteins that belong to the tumor necrosis factor receptor (TNFR) superfamily. Using a generic multiscale simulation platform that spans from atomic to subcellular levels, we compared the detailed physical process of ligand-receptor oligomerization for two specific members in the TNFR superfamily: the complex formed between ligand TNFα and receptor TNFR1 versus the complex formed between ligand TNFβ and receptor TNFR2. Interestingly, although these two systems share high similarity on the tertiary and quaternary structural levels, our results indicate that their oligomers are formed with very different dynamic properties and spatial patterns. We demonstrated that the changes of receptor’s conformational fluctuations due to the membrane confinements are closely related to such difference. Consistent to previous experiments, our simulations also showed that TNFR can preassemble into dimers prior to ligand binding, while the introduction of TNF ligands induced higher-order oligomerization due to a multivalent effect. This study, therefore, provides the molecular basis to TNFR oligomerization and reveals new insights to TNFR-mediated signal transduction. Moreover, our multiscale simulation framework serves as a prototype that paves the way to study higher-order assembly of cell surface receptors in many other bio-systems.

Exploratory polarization facilitates mating partner selection in Saccharomyces cerevisiae

Yeast decode pheromone gradients to locate mating partners, providing a model for chemotropism. How yeast polarize toward a single partner in crowded environments is unclear. Initially, cells often polarize in unproductive directions, but then they relocate the polarity site until two partners’ polarity sites align, whereupon the cells “commit” to each other by stabilizing polarity to promote fusion. Here we address the role of the early mobile polarity sites. We found that commitment by either partner failed if just one partner was defective in generating, orienting, or stabilizing its mobile polarity sites. Mobile polarity sites were enriched for pheromone receptors and G proteins, and we suggest that such sites engage in an exploratory search of the local pheromone landscape, stabilizing only when they detect elevated pheromone levels. Mobile polarity sites were also enriched for pheromone secretion factors, and simulations suggest that only focal secretion at polarity sites would produce high pheromone concentrations at the partner’s polarity site, triggering commitment.

High-Performance Whole-Cell Simulation Exploiting Modular Cell Biology Principles

One of the grand challenges of this century is modeling and simulating a whole cell. Extreme regulation of an extensive quantity of model and simulation data during whole-cell modeling and simulation renders it a computationally expensive research problem in systems biology. In this article, we present a high-performance whole-cell simulation exploiting modular cell biology principles. We prepare the simulation by dividing the unicellular bacterium, Escherichia coli (E. coli), into subcells utilizing the spatially localized densely connected protein clusters/modules. We set up a Brownian dynamics-based parallel whole-cell simulation framework by utilizing the Hamiltonian mechanics-based equations of motion. Though the velocity Verlet integration algorithm possesses the capability of solving the equations of motion, it lacks the ability to capture and deal with particle-collision scenarios. Hence, we propose an algorithm for detecting and resolving both elastic and inelastic collisions and subsequently modify the velocity Verlet integrator by incorporating our algorithm into it. Also, we address the boundary conditions to arrest the molecules' motion outside the subcell. For efficiency, we define one hashing-based data structure called the cellular dictionary to store all of the subcell-related information. A benchmark analysis of our CUDA C/C++ simulation code when tested on E. coli using the CPU-GPU cluster indicates that the computational time requirement decreases with the increase in the number of computing cores and becomes stable at around 128 cores. Additional testing on higher organisms such as rats and humans informs us that our proposed work can be extended to any organism and is scalable for high-end CPU-GPU clusters.

A computational study of co-inhibitory immune complex assembly at the interface between T cells and antigen presenting cells

The activation and differentiation of T-cells are mainly directly by their co-regulatory receptors. T lymphocyte-associated protein-4 (CTLA-4) and programed cell death-1 (PD-1) are two of the most important co-regulatory receptors. Binding of PD-1 and CTLA-4 with their corresponding ligands programed cell death-ligand 1 (PD-L1) and B7 on the antigen presenting cells (APC) activates two central co-inhibitory signaling pathways to suppress T cell functions. Interestingly, recent experiments have identified a new cis-interaction between PD-L1 and B7, suggesting that a crosstalk exists between two co-inhibitory receptors and the two pairs of ligand-receptor complexes can undergo dynamic oligomerization. Inspired by these experimental evidences, we developed a coarse-grained model to characterize the assembling of an immune complex consisting of CLTA-4, B7, PD-L1 and PD-1. These four proteins and their interactions form a small network motif. The temporal dynamics and spatial pattern formation of this network was simulated by a diffusion-reaction algorithm. Our simulation method incorporates the membrane confinement of cell surface proteins and geometric arrangement of different binding interfaces between these proteins. A wide range of binding constants was tested for the interactions involved in the network. Interestingly, we show that the CTLA-4/B7 ligand-receptor complexes can first form linear oligomers, while these oligomers further align together into two-dimensional clusters. Similar phenomenon has also been observed in other systems of cell surface proteins. Our test results further indicate that both co-inhibitory signaling pathways activated by B7 and PD-L1 can be down-regulated by the new cis-interaction between these two ligands, consistent with previous experimental evidences. Finally, the simulations also suggest that the dynamic and the spatial properties of the immune complex assembly are highly determined by the energetics of molecular interactions in the network. Our study, therefore, brings new insights to the co-regulatory mechanisms of T cell activation.

The activation of a T cell can be regulated by the receptors on its surface, such as CTLA-4 and PD-1. People used to think that these two receptors inhibit T cell activation through distinct pathways. However, recent experiments discovered that the ligands of these two receptors, B7 and PD-L1, can interact with each other on the same surface of antigen presenting cells. Here we utilized computational simulations to investigate functional roles of this newly discovered interaction in T cell coregulation. The specific environment of interface between T cell and antigen presenting cell has been taken into account of our model. Ligand and receptors randomly diffuse within this interface area. They further involve in different types of interactions, with each other from the same side or the opposite side of cell surface. Using this method, we found ligands and receptors can not only form complexes, but also aggregate into large-scale clusters. We also demonstrated that the engagement between B7 and PD-L1 can reduce the interactions with their corresponding receptors. This study, therefore, offers new insights to our understanding of signal regulation in T cells.

Quantifying the roles of space and stochasticity in computer simulations for cell biology and cellular biochemistry

Most of the fascinating phenomena studied in cell biology emerge from interactions among highly organized multimolecular structures embedded into complex and frequently dynamic cellular morphologies. For the exploration of such systems, computer simulation has proved to be an invaluable tool, and many researchers in this field have developed sophisticated computational models for application to specific cell biological questions. However, it is often difficult to reconcile conflicting computational results that use different approaches to describe the same phenomenon. To address this issue systematically, we have defined a series of computational test cases ranging from very simple to moderately complex, varying key features of dimensionality, reaction type, reaction speed, crowding, and cell size. We then quantified how explicit spatial and/or stochastic implementations alter outcomes, even when all methods use the same reaction network, rates, and concentrations. For simple cases, we generally find minor differences in solutions of the same problem. However, we observe increasing discordance as the effects of localization, dimensionality reduction, and irreversible enzymatic reactions are combined. We discuss the strengths and limitations of commonly used computational approaches for exploring cell biological questions and provide a framework for decision making by researchers developing new models. As computational power and speed continue to increase at a remarkable rate, the dream of a fully comprehensive computational model of a living cell may be drawing closer to reality, but our analysis demonstrates that it will be crucial to evaluate the accuracy of such models critically and systematically.

High-Throughput Measurements of Intra-Cellular and Secreted Cytokine from Single Spheroids Using Anchored Microfluidic Droplets

While many single-cell approaches have been developed to measure secretions from anchorage-independent cells, these protocols cannot be applied to adherent cells, especially when these cells require to be cultured in 3D formats. Here, a platform to measure secretions from individual spheroids of human mesenchymal stem cells, cultured within microfluidic droplets is introduced. The platform allows to quantify the secretions from hundreds of individual spheroids in each device, by using a secondary droplet to bring functionalized micro-beads in proximity to each spheroid. Vascular endothelial growth factor (VEGF-A) is measured on and a broad distribution of secretion levels within the population of spheroids is observed. The intra-cellular level of VEGF-A on each spheroid, measured through immuno-staining, correlates well with the extra-cellular measurement, indicating that the heterogeneities observed at the spheroid level result from variations at the intra-cellular level. Further, the molecular accumulation within the droplets is modeled and it is found that physical confinement is crucial for measurements of protein secretions. The model predicts that the time to achieve a measurement scales with droplet volume. These first measurements of secretions from individual spheroids provide several new biological and technological insights.

The blending region hybrid framework for the simulation of stochastic reaction-diffusion processes

The simulation of stochastic reaction–diffusion systems using fine-grained representations can become computationally prohibitive when particle numbers become large. If particle numbers are sufficiently high then it may be possible to ignore stochastic fluctuations and use a more efficient coarse-grained simulation approach. Nevertheless, for multiscale systems which exhibit significant spatial variation in concentration, a coarse-grained approach may not be appropriate throughout the simulation domain. Such scenarios suggest a hybrid paradigm in which a computationally cheap, coarse-grained model is coupled to a more expensive, but more detailed fine-grained model, enabling the accurate simulation of the fine-scale dynamics at a reasonable computational cost. In this paper, in order to couple two representations of reaction–diffusion at distinct spatial scales, we allow them to overlap in a ‘blending region’. Both modelling paradigms provide a valid representation of the particle density in this region. From one end of the blending region to the other, control of the implementation of diffusion is passed from one modelling paradigm to another through the use of complementary ‘blending functions’ which scale up or down the contribution of each model to the overall diffusion. We establish the reliability of our novel hybrid paradigm by demonstrating its simulation on four exemplar reaction–diffusion scenarios.

The FLAME-accelerated signalling tool (FaST) for facile parallelisation of flexible agent-based models of cell signalling

Agent-based modelling is particularly adept at modelling complex features of cell signalling pathways, where heterogeneity, stochastic and spatial effects are important, thus increasing our understanding of decision processes in biology in such scenarios. However, agent-based modelling often is computationally prohibitive to implement. Parallel computing, either on central processing units (CPUs) or graphical processing units (GPUs), can provide a means to improve computational feasibility of agent-based applications but generally requires specialist coding knowledge and extensive optimisation. In this paper, we address these challenges through the development and implementation of the FLAME-accelerated signalling tool (FaST), a software that permits easy creation and parallelisation of agent-based models of cell signalling, on CPUs or GPUs. FaST incorporates validated new agent-based methods, for accurate modelling of reaction kinetics and, as proof of concept, successfully converted an ordinary differential equation (ODE) model of apoptosis execution into an agent-based model. We finally parallelised this model through FaST on CPUs and GPUs resulting in an increase in performance of 5.8× (16 CPUs) and 53.9×, respectively. The FaST takes advantage of the communicating X-machine approach used by FLAME and FLAME GPU to allow easy alteration or addition of functionality to parallel applications, but still includes inherent parallelisation optimisation. The FaST, therefore, represents a new and innovative tool to easily create and parallelise bespoke, robust, agent-based models of cell signalling.

3D mesh processing using GAMer 2 to enable reaction-diffusion simulations in realistic cellular geometries

Recent advances in electron microscopy have enabled the imaging of single cells in 3D at nanometer length scale resolutions. An uncharted frontier for in silico biology is the ability to simulate cellular processes using these observed geometries. Enabling such simulations requires watertight meshing of electron micrograph images into 3D volume meshes, which can then form the basis of computer simulations of such processes using numerical techniques such as the finite element method. In this paper, we describe the use of our recently rewritten mesh processing software, GAMer 2, to bridge the gap between poorly conditioned meshes generated from segmented micrographs and boundary marked tetrahedral meshes which are compatible with simulation. We demonstrate the application of a workflow using GAMer 2 to a series of electron micrographs of neuronal dendrite morphology explored at three different length scales and show that the resulting meshes are suitable for finite element simulations. This work is an important step towards making physical simulations of biological processes in realistic geometries routine. Innovations in algorithms to reconstruct and simulate cellular length scale phenomena based on emerging structural data will enable realistic physical models and advance discovery at the interface of geometry and cellular processes. We posit that a new frontier at the intersection of computational technologies and single cell biology is now open.

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.

A computational model for understanding the oligomerization mechanisms of TNF receptor superfamily

By recognizing members in the tumor necrosis factor (TNF) receptor superfamily, TNF ligand proteins function as extracellular cytokines to activate various signaling pathways involved in inflammation, proliferation, and apoptosis. Most ligands in TNF superfamily are trimeric and can simultaneously bind to three receptors on cell surfaces. It has been experimentally observed that the formation of these molecular complexes further triggers the oligomerization of TNF receptors, which in turn regulate the intracellular signaling processes by providing transient compartmentalization in the membrane proximal regions of cytoplasm. In order to decode the molecular mechanisms of oligomerization in TNF receptor superfamily, we developed a new computational method that can physically simulate the spatial-temporal process of binding between TNF ligands and their receptors. The simulations show that the TNF receptors can be organized into hexagonal oligomers. The formation of this spatial pattern is highly dependent not only on the molecular properties such as the affinities of trans and cis binding, but also on the cellular factors such as the concentration of TNF ligands in the extracellular area or the density of TNF receptors on cell surfaces. Moreover, our model suggests that if TNF receptors are pre-organized into dimers before ligand binding, these lateral interactions between receptor monomers can play a positive role in stabilizing the ligand-receptor interactions, as well as in regulating the kinetics of receptor oligomerization. Altogether, this method throws lights on the mechanisms of TNF ligand-receptor interactions in cellular environments.

Applications and Challenges of Machine Learning to Enable Realistic Cellular Simulations

In this perspective, we examine three key aspects of an end-to-end pipeline for realistic cellular simulations: reconstruction and segmentation of cellular structures; generation of cellular structures; and mesh generation, simulation, and data analysis. We highlight some of the relevant prior work in these distinct but overlapping areas, with a particular emphasis on current use of machine learning technologies, as well as on future opportunities.

Stochastic self-tuning hybrid algorithm for reaction-diffusion systems

Many biochemical phenomena involve reactants with vastly different concentrations, some of which are amenable to continuum-level descriptions, while the others are not. We present a hybrid self-tuning algorithm to model such systems. The method combines microscopic (Brownian) dynamics for diffusion with mesoscopic (Gillespie-type) methods for reactions and remains efficient in a wide range of regimes and scenarios with large variations of concentrations. Its accuracy, robustness, and versatility are balanced by redefining propensities and optimizing the mesh size and time step. We use a bimolecular reaction to demonstrate the potential of our method in a broad spectrum of scenarios: from almost completely reaction-dominated systems to cases where reactions rarely occur or take place very slowly. The simulation results show that the number of particles present in the system does not degrade the performance of our method. This makes it an accurate and computationally efficient tool to model complex multireaction systems.

Multiscale simulation unravel the kinetic mechanisms of inflammasome assembly

In the innate immune system, the host defense from the invasion of external pathogens triggers the inflammatory responses. Proteins involved in the inflammatory pathways were often found to aggregate into supramolecular oligomers, called 'inflammasome', mostly through the homotypic interaction between their domains that belong to the death domain superfamily. Although much has been known about the formation of these helical molecular machineries, the detailed correlation between the dynamics of their assembly and the structure of each domain is still not well understood. Using the filament formed by the PYD domains of adaptor molecule ASC as a test system, we constructed a new multiscale simulation framework to study the kinetics of inflammasome assembly. We found that the filament assembly is a multi-step, but highly cooperative process. Moreover, there are three types of binding interfaces between domain subunits in the ASC^PYD filament. The multiscale simulation results suggest that dynamics of domain assembly are rooted in the primary protein sequence which defines the energetics of molecular recognition through three binding interfaces. Interface I plays a more regulatory role than the other two in mediating both the kinetics and the thermodynamics of assembly. Finally, the efficiency of our computational framework allows us to design mutants on a systematic scale and predict their impacts on filament assembly. In summary, this is, to the best of our knowledge, the first simulation method to model the spatial-temporal process of inflammasome assembly. Our work is a useful addition to a suite of existing experimental techniques to study the functions of inflammasome in innate immune system.

Ratiometric GPCR signaling enables directional sensing in yeast

Accurate detection of extracellular chemical gradients is essential for many cellular behaviors. Gradient sensing is challenging for small cells, which can experience little difference in ligand concentrations on the up-gradient and down-gradient sides of the cell. Nevertheless, the tiny cells of the yeast Saccharomyces cerevisiae reliably decode gradients of extracellular pheromones to find their mates. By imaging the behavior of polarity factors and pheromone receptors, we quantified the accuracy of initial polarization during mating encounters. We found that cells bias the orientation of initial polarity up-gradient, even though they have unevenly distributed receptors. Uneven receptor density means that the gradient of ligand-bound receptors does not accurately reflect the external pheromone gradient. Nevertheless, yeast cells appear to avoid being misled by responding to the fraction of occupied receptors rather than simply the concentration of ligand-bound receptors. Such ratiometric sensing also serves to amplify the gradient of active G protein. However, this process is quite error-prone, and initial errors are corrected during a subsequent indecisive phase in which polarity clusters exhibit erratic mobile behavior.

Cells use surface receptors to decode spatial information from chemical gradients, but accurate decoding is hampered by small cell size and the presence of molecular noise. This study shows that yeast cells decode pheromone gradients by measuring the local ratio of bound to unbound receptors. This mechanism corrects for uneven receptor density at the surface and amplifies the gradient transmitted to downstream components.

Accuracy of position determination in Ca^{2+} signaling

A living cell senses its environment and responds to external signals. In this paper, we study theoretically the precision at which cells can determine the position of a spatially localized transient extracellular signal. To this end, we focus on the case where the stimulus is converted into the release of a small molecule that acts as a second messenger, for example, Ca^{2+}, and activates kinases that change the activity of enzymes by phosphorylating them. We analyze the spatial distribution of phosphorylation events using stochastic simulations as well as a mean-field approach. Kinases that need to bind to the cell membrane for getting activated provide more accurate estimates than cytosolic kinases. Our results could explain why the rate of Ca^{2+} detachment from the membrane-binding conventional protein kinase Cα is larger than its phosphorylation rate.

Motion of a Brownian particle in the presence of reactive boundaries

We study the one-dimensional motion of a Brownian particle inside a confinement described by two reactive boundaries which can partially reflect or absorb the particle. Understanding the effects of such boundaries is important in physics, chemistry, and biology. We compute the probability density of the particle displacement exactly, from which we derive expressions for the survival probability and the mean absorption time as a function of the reactive coefficients. Furthermore, using the Feynman-Kac formalism, we investigate the local time profile, which is the fluctuating time spent by the particle at a given location, both till a fixed observation time and till the absorption time. Our analytical results are compared to numerical simulations, showing perfect agreement.

FISIK: Framework for the Inference of In Situ Interaction Kinetics from Single-Molecule Imaging Data

Recent experimental and computational developments have been pushing the limits of live-cell single-molecule imaging, enabling the monitoring of intermolecular interactions in their native environment with high spatiotemporal resolution. However, interactions are captured only for the labeled subset of molecules, which tends to be a small fraction. As a result, it has remained a challenge to calculate molecular interaction kinetics, in particular association rates, from live-cell single-molecule tracking data. To overcome this challenge, we developed a mathematical modeling-based Framework for the Inference of in Situ Interaction Kinetics (FISIK) from single-molecule imaging data with substoichiometric labeling. FISIK consists of (I) devising a mathematical model of molecular movement and interactions, mimicking the biological system and data-acquisition setup, and (II) estimating the unknown model parameters, including molecular association and dissociation rates, by fitting the model to experimental single-molecule data. Due to the stochastic nature of the model and data, we adapted the method of indirect inference for model calibration. We validated FISIK using a series of tests in which we simulated trajectories of diffusing molecules that interact with each other, considering a wide range of model parameters, and including resolution limitations, tracking errors, and mismatches between the model and the biological system it mimics. We found that FISIK has the sensitivity to determine association and dissociation rates, with accuracy and precision depending on the labeled fraction of molecules and the extent of molecule tracking errors. For cases where the labeled fraction is too low (e.g., to afford accurate tracking), combining dynamic but sparse single-molecule imaging data with almost-whole population oligomer distribution data improves FISIK's performance. All in all, FISIK is a promising approach for the derivation of molecular interaction kinetics in their native environment from single-molecule imaging data with substoichiometric labeling.

Positional Information Readout in Ca^{2+} Signaling

Living cells respond to spatially confined signals. Intracellular signal transmission often involves the release of second messengers like Ca^{2+}. They eventually trigger a physiological response, for example, by activating kinases that in turn activate target proteins through phosphorylation. Here, we investigate theoretically how positional information can be accurately read out by protein phosphorylation in spite of rapid second messenger diffusion. We find that accuracy is increased by binding of kinases to the cell membrane prior to phosphorylation and by increasing the rate of Ca^{2+} loss from the cell interior. These findings could explain some salient features of the conventional protein kinase Cα.

Spatial Stochastic Intracellular Kinetics: A Review of Modelling Approaches

Models of chemical kinetics that incorporate both stochasticity and diffusion are an increasingly common tool for studying biology. The variety of competing models is vast, but two stand out by virtue of their popularity: the reaction-diffusion master equation and Brownian dynamics. In this review, we critically address a number of open questions surrounding these models: How can they be justified physically? How do they relate to each other? How do they fit into the wider landscape of chemical models, ranging from the rate equations to molecular dynamics? This review assumes no prior knowledge of modelling chemical kinetics and should be accessible to a wide range of readers.

Rule-Based Modeling Using Wildcards in the Smoldyn Simulator

Many biological molecules exist in multiple variants, such as proteins with different posttranslational modifications, DNAs with different sequences, and phospholipids with different chain lengths. Representing these variants as distinct species, as most biochemical simulators do, leads to the problem that the number of species, and chemical reactions that interconvert them, typically increase combinatorially with the number of ways that the molecules can vary. This can be alleviated by "rule-based modeling methods," in which software generates the chemical reaction network from relatively simple "rules." This chapter presents a new approach to rule-based modeling. It is based on wildcards that match to species names, much as wildcards can match to file names in computer operating systems. It is much simpler to use than the formal rule-based modeling approaches developed previously but can lead to unintended consequences if not used carefully. This chapter demonstrates rule-based modeling with wildcards through examples for signaling systems, protein complexation, polymerization, nucleic acid sequence copying and mutation, the "SMILES" chemical notation, and others. The method is implemented in Smoldyn, a spatial and stochastic biochemical simulator, for both generate-first and on-the-fly expansion, meaning whether the reaction network is generated before or during the simulation.

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.

A Novel A Priori Simulation Algorithm for Absorbing Receivers in Diffusion-Based Molecular Communication Systems

A novel a priori Monte Carlo (APMC) algorithm is proposed to accurately simulate the molecules absorbed at spherical receiver(s) with low-computational complexity in diffusion-based molecular communication (MC) systems. It is demonstrated that the APMC algorithm achieves high-simulation efficiency since by using this algorithm, the fraction of molecules absorbed for a relatively large time step length precisely matches the analytical result. Therefore, the APMC algorithm overcomes the shortcoming of the existing refined Monte Carlo (RMC) algorithm, which enables accurate simulation for a relatively small time step length only. Moreover, for the RMC algorithm, an expression is proposed to quickly predict the simulation accuracy as a function of the time step length and system parameters, which facilitates the choice of simulation time step for a given system. Furthermore, a likelihood threshold is proposed for both the RMC and APMC algorithms to significantly save computational complexity while causing an extremely small loss in accuracy.

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.