Green's-function reaction dynamics: a particle-based approach for simulating biochemical networks in time and space

Simulation of Radiation-Induced DNA Damage and Protection by Histones Using the Code RITRACKS

(1) Background: DNA damage is of great importance in the understanding of the effects of ionizing radiation. Various types of DNA damage can result from exposure to ionizing radiation, with clustered types considered the most important for radiobiological effects. (2) Methods: The code RITRACKS (Relativistic Ion Tracks), a program that simulates stochastic radiation track structures, was used to simulate DNA damage by photons and ions spanning a broad range of linear energy transfer (LET) values. To perform these simulations, the transport code was modified to include cross sections for the interactions of ions or electrons with DNA and amino acids for ionizations, dissociative electron attachment, and elastic collisions. The radiochemistry simulations were performed using a step-by-step algorithm that follows the evolution of all particles in time, including reactions between radicals and DNA structures and amino acids. Furthermore, detailed DNA damage events, such as base pair positions, DNA fragment lengths, and fragment yields, were recorded. (3) Results: We report simulation results using photons and the ions 1H+, 4He2+, 12C6+, 16O8+, and 56Fe26+ at various energies, covering LET values from 0.3 to 164 keV/µm, and performed a comparison with other codes and experimental results. The results show evidence of DNA protection from damage at its points of contacts with histone proteins. (4) Conclusions: RITRACKS can provide a framework for studying DNA damage from a variety of ionizing radiation sources with detailed representations of DNA at the atomic scale, DNA-associated proteins, and resulting DNA damage events and statistics, enabling a broader range of future comparisons with experiments such as those based on DNA sequencing.

Computational Approaches to Predict Protein-Protein Interactions in Crowded Cellular Environments

Investigating protein-protein interactions is crucial for understanding cellular biological processes because proteins often function within molecular complexes rather than in isolation. While experimental and computational methods have provided valuable insights into these interactions, they often overlook a critical factor: the crowded cellular environment. This environment significantly impacts protein behavior, including structural stability, diffusion, and ultimately the nature of binding. In this review, we discuss theoretical and computational approaches that allow the modeling of biological systems to guide and complement experiments and can thus significantly advance the investigation, and possibly the predictions, of protein-protein interactions in the crowded environment of cell cytoplasm. We explore topics such as statistical mechanics for lattice simulations, hydrodynamic interactions, diffusion processes in high-viscosity environments, and several methods based on molecular dynamics simulations. By synergistically leveraging methods from biophysics and computational biology, we review the state of the art of computational methods to study the impact of molecular crowding on protein-protein interactions and discuss its potential revolutionizing effects on the characterization of the human interactome.

Technologies for whole-cell modeling: Genome-wide reconstruction of a cell in silico

With advances in high-throughput, large-scale in vivo measurement and genome modification techniques at the single-nucleotide level, there is an increasing demand for the development of new technologies for the flexible design and control of cellular systems. Computer-aided design is a powerful tool to design new cells. Whole-cell modeling aims to integrate various cellular subsystems, determine their interactions and cooperative mechanisms, and predict comprehensive cellular behaviors by computational simulations on a genome-wide scale. It has been applied to prokaryotes, yeasts, and higher eukaryotic cells, and utilized in a wide range of applications, including production of valuable substances, drug discovery, and controlled differentiation. Whole-cell modeling, consisting of several thousand elements with diverse scales and properties, requires innovative model construction, simulation, and analysis techniques. Furthermore, whole-cell modeling has been extended to multiple scales, including high-resolution modeling at the single-nucleotide and single-amino acid levels and multicellular modeling of tissues and organs. This review presents an overview of the current state of whole-cell modeling, discusses the novel computational and experimental technologies driving it, and introduces further developments toward multihierarchical modeling on a whole-genome scale.

Diffusive exit rates through pores in membrane-enclosed structures

The function of many membrane-enclosed intracellular structures relies on release of diffusing particles that exit through narrow pores or channels in the membrane. The rate of release varies with pore size, density, and length of the channel. We propose a simple approximate model, validated with stochastic simulations, for estimating the effective release rate from cylinders, and other simple-shaped domains, as a function of channel parameters. The results demonstrate that, for very small pores, a low density of channels scattered over the boundary is sufficient to achieve substantial rates of particle release. Furthermore, we show that increasing the length of passive channels will both reduce release rates and lead to a less steep dependence on channel density. Our results are compared to previously-measured local calcium release rates from tubules of the endoplasmic reticulum, providing an estimate of the relevant channel density responsible for the observed calcium efflux.

Particle-based model of mechanosensory contractility kit assembly

Cell shape change processes, such as proliferation, polarization, migration, and cancer metastasis, rely on a dynamic network of macromolecules. The proper function of this network enables mechanosensation, the ability of cells to sense and respond to mechanical cues. Myosin II and cortexillin I, critical elements of the cellular mechanosensory machinery, preassemble in the cytoplasm of Dictyostelium cells into complexes that we have termed contractility kits (CKs). Two IQGAP proteins then differentially regulate the mechanoresponsiveness of the cortexillin I-myosin II elements within CKs. To investigate the mechanism of CK self-assembly and gain insight into possible molecular means for IQGAP regulation, we developed a coarse-grained excluded volume molecular model in which all protein polymers are represented by nm-sized spheres connected by spring-like links. The model is parameterized using experimentally measured parameters acquired through fluorescence cross-correlation spectroscopy and fluorescence correlation spectroscopy, which describe the interaction affinities and diffusion coefficients for individual molecular components, and which have also been validated via several orthogonal methods. Simulations of wild-type and null-mutant conditions implied that the temporal order of assembly of these kits is dominated by myosin II dimer formation and that IQGAP proteins mediate cluster growth. In addition, our simulations predicted the existence of "ambiguous" CKs that incorporate both classes of IQGAPs, and we confirmed this experimentally using fluorescence cross-correlation spectroscopy. The model serves to describe the formation of the CKs and how their assembly enables and regulates mechanosensation at the molecular level.

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

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

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.

Sliding across a surface: Particles with fixed and mobile ligands

A quantitative model of the mobility of ligand-presenting particles at the interface is pivotal to understanding important systems in biology and nanotechnology. In this work, we investigate the emerging dynamics of particles featuring ligands that selectively bind receptors decorating an interface. The formation of a ligand-receptor complex leads to a molecular bridge anchoring the particle to the surface. We consider systems with reversible bridges in which ligand-receptor pairs bind/unbind with finite reaction rates. For a given set of bridges, the particle can explore a tiny fraction of the surface as the extensivity of the bridges is finite. We show how, at timescales longer than the bridges' lifetime, the average position of the particle diffuses away from its initial value. We distill our findings into two analytic equations for the sliding diffusion constant of particles carrying mobile and fixed ligands. We quantitatively validate our theoretical predictions using reaction-diffusion simulations. We compare our findings with results from recent literature studies and discuss the molecular parameters that likely affect the particle's mobility most. Our results, along with recent literature studies, will allow inferring the microscopic parameters at play in complex biological systems from experimental trajectories.

Corrections to reaction-diffusion dynamics above the upper critical dimension

Reaction-diffusion models are common in many areas of statistical physics, where they describe the late-time dynamics of chemical reactions. Using a Bose gas representation, which maps the real-time dynamics of the reactants to the imaginary-time evolution of an interacting Bose gas, we consider corrections to the late-time scaling of k-particle annihilation processes kA→∅ above the upper critical dimension, where mean-field theory sets the leading order. We establish that the leading corrections are not given by a small renormalization of the reaction rate due to k-particle memory effects, but instead set by higher-order correlation functions that capture memory effects of subclusters of reactants. Drawing on methods developed for ultracold quantum gases and nuclear physics, we compute these corrections exactly for various annihilation processes with k>2.

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.

Space-time histories approach to fast stochastic simulation of bimolecular reactions

Computational models of reaction-diffusion systems involving low copy numbers or strongly heterogeneous molecular spatial distributions, such as those frequently found in cellular signaling pathways, require approaches that account for the stochastic dynamics of individual particles, as opposed to approaches representing them through their average concentrations. Efforts to remedy the high computational cost associated with particle-based stochastic approaches by taking advantage of Green's functions are hampered by the need to draw random numbers from complicated, and therefore costly, non-standard probability distributions to update particle positions. Here, we introduce an approach that permits the reconstruction of entire molecular trajectories, including bimolecular encounters, in retrospect, after a simulated time step, while avoiding inefficient draws from non-standard distributions. This means that highly accurate stochastic simulations can be performed for system sizes that would be prohibitively costly to simulate with conventional Green's function based methods. The algorithm applies equally well to one, two, and three dimensional systems and can be readily extended to include deterministic forces specified by an interaction potential, such as the Coulomb potential.

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.

Incorporating domain growth into hybrid methods for reaction-diffusion systems

Reaction-diffusion mechanisms are a robust paradigm that can be used to represent many biological and physical phenomena over multiple spatial scales. Applications include intracellular dynamics, the migration of cells and the patterns formed by vegetation in semi-arid landscapes. Moreover, domain growth is an important process for embryonic growth and wound healing. There are many numerical modelling frameworks capable of simulating such systems on growing domains; however, each of these may be well suited to different spatial scales and particle numbers. Recently, spatially extended hybrid methods on static domains have been produced to bridge the gap between these different modelling paradigms in order to represent multi-scale phenomena. However, such methods have not been developed with domain growth in mind. In this paper, we develop three hybrid methods on growing domains, extending three of the prominent static-domain hybrid methods. We also provide detailed algorithms to allow others to employ them. We demonstrate that the methods are able to accurately model three representative reaction-diffusion systems accurately and without bias.

Concentration sensing in crowded environments

Signal transduction within crowded cellular compartments is essential for the physiological function of cells. Although the accuracy with which receptors can probe the concentration of ligands has been thoroughly investigated in dilute systems, the effect of macromolecular crowding on the inference of concentration remains unclear. In this work, we develop an algorithm to simulate reversible reactions between reacting Brownian particles. Our algorithm facilitates the calculation of reaction rates and correlation times for ligand-receptor systems in the presence of macromolecular crowding. Using this method, we show that it is possible for crowding to increase the accuracy of estimated ligand concentration based on receptor occupancy. In particular, we find that crowding can enhance the effective association rates between small ligands and receptors to a degree sufficient to overcome the increased chance of rebinding due to caging by crowding molecules. For larger ligands, crowding decreases the accuracy of the receptor's estimate primarily by decreasing the microscopic association and dissociation rates.

Stochastic simulation algorithms for Interacting Particle Systems

Interacting Particle Systems (IPSs) are used to model spatio-temporal stochastic systems in many disparate areas of science. We design an algorithmic framework that reduces IPS simulation to simulation of well-mixed Chemical Reaction Networks (CRNs). This framework minimizes the number of associated reaction channels and decouples the computational cost of the simulations from the size of the lattice. Decoupling allows our software to make use of a wide class of techniques typically reserved for well-mixed CRNs. We implement the direct stochastic simulation algorithm in the open source programming language Julia. We also apply our algorithms to several complex spatial stochastic phenomena. including a rock-paper-scissors game, cancer growth in response to immunotherapy, and lipid oxidation dynamics. Our approach aids in standardizing mathematical models and in generating hypotheses based on concrete mechanistic behavior across a wide range of observed spatial phenomena.

A review of simulation codes and approaches for radiation chemistry

Historically, the field of radiation chemistry began shortly after the discovery of radioactivity, and its development has been closely related to discoveries in other related fields such as radiation and nuclear physics. Radiolysis of water and radiation chemistry have been very important in elucidating how radiation affects living matter and how it induces DNA damage. Nowadays, we recognize the importance of chemistry to understanding the effects of radiation on cells; however, it took several decades to obtain this insight, and much is still unknown. The radiolysis of water and aqueous solutions have been the subject of much experimental and theoretical research for many decades. One important concept closely related to radiation chemistry is radiation track structure. Track structure results from early physical and physicochemical events that lead to a highly non-homogenous distribution of radiolytic species. Because ionizing radiation creates unstable species that are distributed non-homogenously, the use of conventional reaction kinetics methods does not describe this chemistry well. In recent years, several methods have been developed for simulating radiation chemistry. In this review, we give a brief history of the field and the development of the simulation codes. We review the current methods used to simulate radiolysis of water and radiation chemistry, and we describe several radiation chemistry codes and their applications.

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.

How Anionic Lipids Affect Spatiotemporal Properties of KRAS4B on Model Membranes

RAS proteins are small membrane-anchored GTPases that regulate key cellular signaling networks. It has been recently shown that different anionic lipid types can affect the spatiotemporal properties of RAS through dimerization/clustering and signaling fidelity. To understand the effects of anionic lipids on key spatiotemporal properties of RAS, we dissected 1 ms of data from all-atom molecular dynamics simulations for KRAS4B on two model anionic lipid membranes that have 30% of POPS mixed with neutral POPC and 8% of PIP2 mixed with POPC. We unveiled the orientation space of KRAS4B, whose kinetics were slower and more distinguishable on the membrane containing PIP2 than the membrane containing POPS. Particularly, the PIP2-mixed membrane can differentiate a third kinetic orientation state from the other two known orientation states. We observed that each orientation state may yield different binding modes with an RAF kinase, which is required for activating the MAPK/ERK signaling pathway. However, an overall occluded probability, for which RAF kinases cannot bind KRAS4B, remains unchanged on the two different membranes. We identified rare fast diffusion modes of KRAS4B that appear coupled with orientations exposed to cytosolic RAF. Particularly, on the membrane having PIP2, we found nonlinear correlations between the orientation states and the conformations of the cationic farnesylated hypervariable region, which acts as an anchor in the membrane. Using diffusion coefficients estimated from the all-atom simulations, we quantified the effect of PIP2 and POPS on the KRAS4B dimerization via Green's function reaction dynamics simulations, in which the averaged dimerization rate is 12.5% slower on PIP2-mixed membranes.

Presence or Absence of Ras Dimerization Shows Distinct Kinetic Signature in Ras-Raf Interaction

Initiations of cell signaling pathways often occur through the formation of multiprotein complexes that form through protein-protein interactions. Therefore, detecting their presence is central to understanding the function of a cell signaling pathway, aberration of which often leads to fatal diseases, including cancers. However, the multiprotein complexes are often difficult to detect using microscopes due to their small sizes. Therefore, currently, their presence can be only detected through indirect means. In this article, we propose to investigate the presence or absence of protein complexes through some easily measurable kinetic parameters, such as activation rates. As a proof of concept, we investigate the Ras-Raf system, a well-characterized cell signaling system. It has been hypothesized that Ras dimerization is necessary to create activated Raf dimers. Although there are circumstantial evidences supporting the Ras dimerization hypothesis, direct proof of Ras dimerization is still inconclusive. In the absence of conclusive direct experimental proof, this hypothesis can only be examined through indirect evidences of Ras dimerization. In this article, using a multiscale simulation technique, we provide multiple criteria that distinguishes an activation mechanism involving Ras dimerization from another mechanism that does not involve Ras dimerization. The provided criteria will be useful in the investigation of not only Ras-Raf interaction but also other two-protein interactions.

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.

An implicit lipid model for efficient reaction-diffusion simulations of protein binding to surfaces of arbitrary topology

Localization of proteins to a membrane is an essential step in a broad range of biological processes such as signaling, virion formation, and clathrin-mediated endocytosis. The strength and specificity of proteins binding to a membrane depend on the lipid composition. Single-particle reaction-diffusion methods offer a powerful tool for capturing lipid-specific binding to membrane surfaces by treating lipids explicitly as individual diffusible binding sites. However, modeling lipid particle populations is expensive. Here, we present an algorithm for reversible binding of proteins to continuum surfaces with implicit lipids, providing dramatic speed-ups to many body simulations. Our algorithm can be readily integrated into most reaction-diffusion software packages. We characterize changes to kinetics that emerge from explicit vs implicit lipids as well as surface adsorption models, showing excellent agreement between our method and the full explicit lipid model. Compared to models of surface adsorption, which couple together binding affinity and lipid concentration, our implicit lipid model decouples them to provide more flexibility for controlling surface binding properties and lipid inhomogeneity, thus reproducing binding kinetics and equilibria. Crucially, we demonstrate our method's application to membranes of arbitrary curvature and topology, modeled via a subdivision limit surface, again showing excellent agreement with explicit lipid simulations. Unlike adsorption models, our method retains the ability to bind lipids after proteins are localized to the surface (through, e.g., a protein-protein interaction), which can greatly increase the stability of multiprotein complexes on the surface. Our method will enable efficient cell-scale simulations involving proteins localizing to realistic membrane models, which is a critical step for predictive modeling and quantification of in vitro and in vivo dynamics.

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.

Particle-Based Simulation Reveals Macromolecular Crowding Effects on the Michaelis-Menten Mechanism

Many computational models for analyzing and predicting cell physiology rely on in vitro data collected in dilute and controlled buffer solutions. However, this can mislead models because up to 40% of the intracellular volume—depending on the organism, the physiology, and the cellular compartment—is occupied by a dense mixture of proteins, lipids, polysaccharides, RNA, and DNA. These intracellular macromolecules interfere with the interactions of enzymes and their reactants and thus affect the kinetics of biochemical reactions, making in vivo reactions considerably more complex than the in vitro data indicates. In this work, we present a new, to our knowledge, type of kinetics that captures and quantifies the effect of volume exclusion and other spatial phenomena on the kinetics of elementary reactions. We further developed a framework that allows for the efficient parameterization of these kinetics using particle simulations. Our formulation, entitled generalized elementary kinetics, can be used to analyze and predict the effect of intracellular crowding on enzymatic reactions and was herein applied to investigate the influence of crowding on phosphoglycerate mutase in Escherichia coli, which exhibits prototypical reversible Michaelis-Menten kinetics. Current research indicates that many enzymes are reaction limited and not diffusion limited, and our results suggest that the influence of fractal diffusion is minimal for these reaction-limited enzymes. Instead, increased association rates and decreased dissociation rates lead to a strong decrease in the effective maximal velocities V_max and the effective Michaelis-Menten constants K_M under physiologically relevant volume occupancies. Finally, the effects of crowding were explored in the context of a linear pathway, with the finding that crowding can have a redistributing effect on the effective flux responses in the case of twofold enzyme overexpression. We suggest that this framework, in combination with detailed kinetics models, will improve our understanding of enzyme reaction networks under nonideal conditions.

Minimal coarse-grained models for molecular self-organisation in biology

The molecular machinery of life is largely created via self-organisation of individual molecules into functional assemblies. Minimal coarse-grained models, in which a whole macromolecule is represented by a small number of particles, can be of great value in identifying the main driving forces behind self-organisation in cell biology. Such models can incorporate data from both molecular and continuum scales, and their results can be directly compared to experiments. Here we review the state of the art of models for studying the formation and biological function of macromolecular assemblies in living organisms. We outline the key ingredients of each model and their main findings. We illustrate the contribution of this class of simulations to identifying the physical mechanisms behind life and diseases, and discuss their future developments.

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 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.

eGFRD in all dimensions

Biochemical reactions often occur at low copy numbers but at once in crowded and diverse environments. Space and stochasticity therefore play an essential role in biochemical networks. Spatial-stochastic simulations have become a prominent tool for understanding how stochasticity at the microscopic level influences the macroscopic behavior of such systems. While particle-based models guarantee the level of detail necessary to accurately describe the microscopic dynamics at very low copy numbers, the algorithms used to simulate them typically imply trade-offs between computational efficiency and biochemical accuracy. eGFRD (enhanced Green's Function Reaction Dynamics) is an exact algorithm that evades such trade-offs by partitioning the N-particle system into M ≤ N analytically tractable one- and two-particle systems; the analytical solutions (Green's functions) then are used to implement an event-driven particle-based scheme that allows particles to make large jumps in time and space while retaining access to their state variables at arbitrary simulation times. Here we present "eGFRD2," a new eGFRD version that implements the principle of eGFRD in all dimensions, thus enabling efficient particle-based simulation of biochemical reaction-diffusion processes in the 3D cytoplasm, on 2D planes representing membranes, and on 1D elongated cylinders representative of, e.g., cytoskeletal tracks or DNA; in 1D, it also incorporates convective motion used to model active transport. We find that, for low particle densities, eGFRD2 is up to 6 orders of magnitude faster than conventional Brownian dynamics. We exemplify the capabilities of eGFRD2 by simulating an idealized model of Pom1 gradient formation, which involves 3D diffusion, active transport on microtubules, and autophosphorylation on the membrane, confirming recent experimental and theoretical results on this system to hold under genuinely stochastic conditions.

ReaDDy 2: Fast and flexible software framework for interacting-particle reaction dynamics

Interacting-particle reaction dynamics (iPRD) combines the simulation of dynamical trajectories of interacting particles as in molecular dynamics (MD) simulations with reaction kinetics, in which particles appear, disappear, or change their type and interactions based on a set of reaction rules. This combination facilitates the simulation of reaction kinetics in crowded environments, involving complex molecular geometries such as polymers, and employing complex reaction mechanisms such as breaking and fusion of polymers. iPRD simulations are ideal to simulate the detailed spatiotemporal reaction mechanism in complex and dense environments, such as in signalling processes at cellular membranes, or in nano- to microscale chemical reactors. Here we introduce the iPRD software ReaDDy 2, which provides a Python interface in which the simulation environment, particle interactions and reaction rules can be conveniently defined and the simulation can be run, stored and analyzed. A C++ interface is available to enable deeper and more flexible interactions with the framework. The main computational work of ReaDDy 2 is done in hardware-specific simulation kernels. While the version introduced here provides single- and multi-threading CPU kernels, the architecture is ready to implement GPU and multi-node kernels. We demonstrate the efficiency and validity of ReaDDy 2 using several benchmark examples. ReaDDy 2 is available at the https://readdy.github.io/ website.

Rate constants for proteins binding to substrates with multiple binding sites using a generalized forward flux sampling expression

To predict the response of a biochemical system, knowledge of the intrinsic and effective rate constants of proteins is crucial. The experimentally accessible effective rate constant for association can be decomposed in a diffusion-limited rate at which proteins come into contact and an intrinsic association rate at which the proteins in contact truly bind. Reversely, when dissociating, bound proteins first separate into a contact pair with an intrinsic dissociation rate, before moving away by diffusion. While microscopic expressions exist that enable the calculation of the intrinsic and effective rate constants by conducting a single rare event simulation of the protein dissociation reaction, these expressions are only valid when the substrate has just one binding site. If the substrate has multiple binding sites, a bound enzyme can, besides dissociating into the bulk, also hop to another binding site. Calculating transition rate constants between multiple states with forward flux sampling requires a generalized rate expression. We present this expression here and use it to derive explicit expressions for all intrinsic and effective rate constants involving binding to multiple states, including rebinding. We illustrate our approach by computing the intrinsic and effective association, dissociation, and hopping rate constants for a system in which a patchy particle model enzyme binds to a substrate with two binding sites. We find that these rate constants increase as a function of the rotational diffusion constant of the particles. The hopping rate constant decreases as a function of the distance between the binding sites. Finally, we find that blocking one of the binding sites enhances both association and dissociation rate constants. Our approach and results are important for understanding and modeling association reactions in enzyme-substrate systems and other patchy particle systems and open the way for large multiscale simulations of such systems.

Modeling the Self-Assembly of Protein Complexes through a Rigid-Body Rotational Reaction-Diffusion Algorithm

The reaction-diffusion equations provide a powerful framework for modeling nonequilibrium, cell-scale dynamics over the long time scales that are inaccessible by traditional molecular modeling approaches. Single-particle reaction-diffusion offers the highest resolution technique for tracking such dynamics, but it has not been applied to the study of protein self-assembly due to its treatment of reactive species as single-point particles. Here, we develop a relatively simple but accurate approach for building rigid structure and rotation into single-particle reaction-diffusion methods, providing a rate-based method for studying protein self-assembly. Our simplifying assumption is that reactive collisions can be evaluated purely on the basis of the separations between the sites, and not their orientations. The challenge of evaluating reaction probabilities can then be performed using well-known equations based on translational diffusion in both 3D and 2D, by employing an effective diffusion constant we derive here. We show how our approach reproduces both the kinetics of association, which is altered by rotational diffusion, and the equilibrium of reversible association, which is not. Importantly, the macroscopic kinetics of association can be predicted on the basis of the microscopic parameters of our structurally resolved model, allowing for critical comparisons with theory and other rate-based simulations. We demonstrate this method for efficient, rate-based simulations of self-assembly of clathrin trimers, highlighting how formation of regular lattices impacts the kinetics of association.

Reversible Interacting-Particle Reaction Dynamics

Interacting-particle reaction dynamics (iPRD) simulates the spatiotemporal evolution of particles that experience interaction forces and can react with one another. The combination of interaction forces and reactions enables a wide range of complex reactive systems in biology and chemistry to be simulated, but gives rise to new questions such as how to evolve the dynamical equations in a computationally efficient and statistically correct manner. Here we consider reversible reactions such as A + B ⇄ C with interacting particles and derive expressions for the microscopic iPRD simulation parameters such that desired values for the equilibrium constant and the dissociation rate are obtained in the dilute limit. We then introduce a Monte Carlo algorithm that ensures detailed balance in the iPRD time-evolution (iPRD-DB). iPRD-DB guarantees the correct thermodynamics at all concentrations and maintains the desired kinetics in the dilute limit, where chemical rates are well-defined and kinetic measurement experiments usually operate. We show that in dense particle systems, the incorporation of detailed balance is essential to obtain physically realistic solutions. iPRD-DB is implemented in ReaDDy 2.

Grand canonical diffusion-influenced reactions: A stochastic theory with applications to multiscale reaction-diffusion simulations

Smoluchowski-type models for diffusion-influenced reactions (A + B → C) can be formulated within two frameworks: the probabilistic-based approach for a pair A, B of reacting particles and the concentration-based approach for systems in contact with a bath that generates a concentration gradient of B particles that interact with A. Although these two approaches are mathematically similar, it is not straightforward to establish a precise mathematical relationship between them. Determining this relationship is essential to derive particle-based numerical methods that are quantitatively consistent with bulk concentration dynamics. In this work, we determine the relationship between the two approaches by introducing the grand canonical Smoluchowski master equation (GC-SME), which consists of a continuous-time Markov chain that models an arbitrary number of B particles, each one of them following Smoluchowski's probabilistic dynamics. We show that the GC-SME recovers the concentration-based approach by taking either the hydrodynamic or the large copy number limit. In addition, we show that the GC-SME provides a clear statistical mechanical interpretation of the concentration-based approach and yields an emergent chemical potential for nonequilibrium spatially inhomogeneous reaction processes. We further exploit the GC-SME robust framework to accurately derive multiscale/hybrid numerical methods that couple particle-based reaction-diffusion simulations with bulk concentration descriptions, as described in detail through two computational implementations.

The auxiliary region method: a hybrid method for coupling PDE- and Brownian-based dynamics for reaction-diffusion systems

Reaction-diffusion systems are used to represent many biological and physical phenomena. They model the random motion of particles (diffusion) and interactions between them (reactions). Such systems can be modelled at multiple scales with varying degrees of accuracy and computational efficiency. When representing genuinely multiscale phenomena, fine-scale models can be prohibitively expensive, whereas coarser models, although cheaper, often lack sufficient detail to accurately represent the phenomenon at hand. Spatial hybrid methods couple two or more of these representations in order to improve efficiency without compromising accuracy. In this paper, we present a novel spatial hybrid method, which we call the auxiliary region method (ARM), which couples PDE- and Brownian-based representations of reaction-diffusion systems. Numerical PDE solutions on one side of an interface are coupled to Brownian-based dynamics on the other side using compartment-based 'auxiliary regions'. We demonstrate that the hybrid method is able to simulate reaction-diffusion dynamics for a number of different test problems with high accuracy. Furthermore, we undertake error analysis on the ARM which demonstrates that it is robust to changes in the free parameters in the model, where previous coupling algorithms are not. In particular, we envisage that the method will be applicable for a wide range of spatial multi-scales problems including filopodial dynamics, intracellular signalling, embryogenesis and travelling wave phenomena.

The intrinsic rate constants in diffusion-influenced reactions

Intrinsic rate constants play a dominant role in the theory of diffusion-influenced reactions, but usually as abstract quantities that are implicitly assumed to be known. However, recently it has become clear that modeling complex processes requires explicit knowledge of these intrinsic rates. In this paper we provide microscopic expressions for the intrinsic rate constants for association and dissociation processes of isotropically interacting particles and illustrate how these rates can be computed efficiently using rare event simulations techniques. In addition, we address the role of the orientational dynamics, for particles interacting via anisotropic potentials.

Multiscale simulations of anisotropic particles combining molecular dynamics and Green's function reaction dynamics

The modeling of complex reaction-diffusion processes in, for instance, cellular biochemical networks or self-assembling soft matter can be tremendously sped up by employing a multiscale algorithm which combines the mesoscopic Green's Function Reaction Dynamics (GFRD) method with explicit stochastic Brownian, Langevin, or deterministic molecular dynamics to treat reactants at the microscopic scale [A. Vijaykumar, P. G. Bolhuis, and P. R. ten Wolde, J. Chem. Phys. 143, 214102 (2015)]. Here we extend this multiscale MD-GFRD approach to include the orientational dynamics that is crucial to describe the anisotropic interactions often prevalent in biomolecular systems. We present the novel algorithm focusing on Brownian dynamics only, although the methodology is generic. We illustrate the novel algorithm using a simple patchy particle model. After validation of the algorithm, we discuss its performance. The rotational Brownian dynamics MD-GFRD multiscale method will open up the possibility for large scale simulations of protein signalling networks.

MSM/RD: Coupling Markov state models of molecular kinetics with reaction-diffusion simulations

Molecular dynamics (MD) simulations can model the interactions between macromolecules with high spatiotemporal resolution but at a high computational cost. By combining high-throughput MD with Markov state models (MSMs), it is now possible to obtain long time-scale behavior of small to intermediate biomolecules and complexes. To model the interactions of many molecules at large length scales, particle-based reaction-diffusion (RD) simulations are more suitable but lack molecular detail. Thus, coupling MSMs and RD simulations (MSM/RD) would be highly desirable, as they could efficiently produce simulations at large time and length scales, while still conserving the characteristic features of the interactions observed at atomic detail. While such a coupling seems straightforward, fundamental questions are still open: Which definition of MSM states is suitable? Which protocol to merge and split RD particles in an association/dissociation reaction will conserve the correct bimolecular kinetics and thermodynamics? In this paper, we make the first step toward MSM/RD by laying out a general theory of coupling and proposing a first implementation for association/dissociation of a protein with a small ligand (A + B ⇌ C). Applications on a toy model and CO diffusion into the heme cavity of myoglobin are reported.

Efficient reactive Brownian dynamics

We develop a Split Reactive Brownian Dynamics (SRBD) algorithm for particle simulations of reaction-diffusion systems based on the Doi or volume reactivity model, in which pairs of particles react with a specified Poisson rate if they are closer than a chosen reactive distance. In our Doi model, we ensure that the microscopic reaction rules for various association and dissociation reactions are consistent with detailed balance (time reversibility) at thermodynamic equilibrium. The SRBD algorithm uses Strang splitting in time to separate reaction and diffusion and solves both the diffusion-only and reaction-only subproblems exactly, even at high packing densities. To efficiently process reactions without uncontrolled approximations, SRBD employs an event-driven algorithm that processes reactions in a time-ordered sequence over the duration of the time step. A grid of cells with size larger than all of the reactive distances is used to schedule and process the reactions, but unlike traditional grid-based methods such as reaction-diffusion master equation algorithms, the results of SRBD are statistically independent of the size of the grid used to accelerate the processing of reactions. We use the SRBD algorithm to compute the effective macroscopic reaction rate for both reaction-limited and diffusion-limited irreversible association in three dimensions and compare to existing theoretical predictions at low and moderate densities. We also study long-time tails in the time correlation functions for reversible association at thermodynamic equilibrium and compare to recent theoretical predictions. Finally, we compare different particle and continuum methods on a model exhibiting a Turing-like instability and pattern formation. Our studies reinforce the common finding that microscopic mechanisms and correlations matter for diffusion-limited systems, making continuum and even mesoscopic modeling of such systems difficult or impossible. We also find that for models in which particles diffuse off lattice, such as the Doi model, reactions lead to a spurious enhancement of the effective diffusion coefficients.

An efficient multi-scale Green's function reaction dynamics scheme

Molecular Dynamics-Green's Function Reaction Dynamics (MD-GFRD) is a multiscale simulation method for particle dynamics or particle-based reaction-diffusion dynamics that is suited for systems involving low particle densities. Particles in a low-density region are just diffusing and not interacting. In this case, one can avoid the costly integration of microscopic equations of motion, such as molecular dynamics (MD), and instead turn to an event-based scheme in which the times to the next particle interaction and the new particle positions at that time can be sampled. At high (local) concentrations, however, e.g., when particles are interacting in a nontrivial way, particle positions must still be updated with small time steps of the microscopic dynamical equations. The efficiency of a multi-scale simulation that uses these two schemes largely depends on the coupling between them and the decisions when to switch between the two scales. Here we present an efficient scheme for multi-scale MD-GFRD simulations. It has been shown that MD-GFRD schemes are more efficient than brute-force molecular dynamics simulations up to a molar concentration of 10² μM. In this paper, we show that the choice of the propagation domains has a relevant impact on the computational performance. Domains are constructed using a local optimization of their sizes and a minimal domain size is proposed. The algorithm is shown to be more efficient than brute-force Brownian dynamics simulations up to a molar concentration of 10³ μM and is up to an order of magnitude more efficient compared with previous MD-GFRD schemes.

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.

Simulating biological processes: stochastic physics from whole cells to colonies

The last few decades have revealed the living cell to be a crowded spatially heterogeneous space teeming with biomolecules whose concentrations and activities are governed by intrinsically random forces. It is from this randomness, however, that a vast array of precisely timed and intricately coordinated biological functions emerge that give rise to the complex forms and behaviors we see in the biosphere around us. This seemingly paradoxical nature of life has drawn the interest of an increasing number of physicists, and recent years have seen stochastic modeling grow into a major subdiscipline within biological physics. Here we review some of the major advances that have shaped our understanding of stochasticity in biology. We begin with some historical context, outlining a string of important experimental results that motivated the development of stochastic modeling. We then embark upon a fairly rigorous treatment of the simulation methods that are currently available for the treatment of stochastic biological models, with an eye toward comparing and contrasting their realms of applicability, and the care that must be taken when parameterizing them. Following that, we describe how stochasticity impacts several key biological functions, including transcription, translation, ribosome biogenesis, chromosome replication, and metabolism, before considering how the functions may be coupled into a comprehensive model of a 'minimal cell'. Finally, we close with our expectation for the future of the field, focusing on how mesoscopic stochastic methods may be augmented with atomic-scale molecular modeling approaches in order to understand life across a range of length and time scales.

Particle-based simulations of polarity establishment reveal stochastic promotion of Turing pattern formation

Polarity establishment, the spontaneous generation of asymmetric molecular distributions, is a crucial component of many cellular functions. Saccharomyces cerevisiae (yeast) undergoes directed growth during budding and mating, and is an ideal model organism for studying polarization. In yeast and many other cell types, the Rho GTPase Cdc42 is the key molecular player in polarity establishment. During yeast polarization, multiple patches of Cdc42 initially form, then resolve into a single front. Because polarization relies on strong positive feedback, it is likely that the amplification of molecular-level fluctuations underlies the generation of multiple nascent patches. In the absence of spatial cues, these fluctuations may be key to driving polarization. Here we used particle-based simulations to investigate the role of stochastic effects in a Turing-type model of yeast polarity establishment. In the model, reactions take place either between two molecules on the membrane, or between a cytosolic and a membrane-bound molecule. Thus, we developed a computational platform that explicitly simulates molecules at and near the cell membrane, and implicitly handles molecules away from the membrane. To evaluate stochastic effects, we compared particle simulations to deterministic reaction-diffusion equation simulations. Defining macroscopic rate constants that are consistent with the microscopic parameters for this system is challenging, because diffusion occurs in two dimensions and particles exchange between the membrane and cytoplasm. We address this problem by empirically estimating macroscopic rate constants from appropriately designed particle-based simulations. Ultimately, we find that stochastic fluctuations speed polarity establishment and permit polarization in parameter regions predicted to be Turing stable. These effects can operate at Cdc42 abundances expected of yeast cells, and promote polarization on timescales consistent with experimental results. To our knowledge, our work represents the first particle-based simulations of a model for yeast polarization that is based on a Turing mechanism.

Many cells need to generate and maintain biochemical signals in specific subcellular regions. This phenomenon is broadly called polarity establishment, and is important in fundamental processes such as cell migration and differentiation. A key polarity factor found in diverse organisms, including yeast and humans, is the protein Cdc42. In yeast, Cdc42-dependent polarization occurs through a self-reinforcing biochemical signaling loop. Directional cues can guide polarity establishment, but interestingly, yeast can polarize in the absence of such a cue. The mechanism thought to underlie this symmetry breaking involves the amplification of inhomogeneities in molecular distributions that arise from molecular-level fluctuations. We investigated the effects of random fluctuations on polarization by performing particle-based simulations of the Cdc42 signaling network. We found that fluctuations can facilitate polarization, allowing faster polarization, and polarization over a broader range of concentrations. Our observations may help understand how polarity works in other systems.

Mixed analytical-stochastic simulation method for the recovery of a Brownian gradient source from probability fluxes to small windows

Is it possible to recover the position of a source from the steady-state fluxes of Brownian particles to small absorbing windows located on the boundary of a domain? To address this question, we develop a numerical procedure to avoid tracking Brownian trajectories in the entire infinite space. Instead, we generate particles near the absorbing windows, computed from the analytical expression of the exit probability. When the Brownian particles are generated by a steady-state gradient at a single point, we compute asymptotically the fluxes to small absorbing holes distributed on the boundary of half-space and on a disk in two dimensions, which agree with stochastic simulations. We also derive an expression for the splitting probability between small windows using the matched asymptotic method. Finally, when there are more than two small absorbing windows, we show how to reconstruct the position of the source from the diffusion fluxes. The present approach provides a computational first principle for the mechanism of sensing a gradient of diffusing particles, a ubiquitous problem in cell biology.

Spatially extended hybrid methods: a review

Many biological and physical systems exhibit behaviour at multiple spatial, temporal or population scales. Multiscale processes provide challenges when they are to be simulated using numerical techniques. While coarser methods such as partial differential equations are typically fast to simulate, they lack the individual-level detail that may be required in regions of low concentration or small spatial scale. However, to simulate at such an individual level throughout a domain and in regions where concentrations are high can be computationally expensive. Spatially coupled hybrid methods provide a bridge, allowing for multiple representations of the same species in one spatial domain by partitioning space into distinct modelling subdomains. Over the past 20 years, such hybrid methods have risen to prominence, leading to what is now a very active research area across multiple disciplines including chemistry, physics and mathematics. There are three main motivations for undertaking this review. Firstly, we have collated a large number of spatially extended hybrid methods and presented them in a single coherent document, while comparing and contrasting them, so that anyone who requires a multiscale hybrid method will be able to find the most appropriate one for their need. Secondly, we have provided canonical examples with algorithms and accompanying code, serving to demonstrate how these types of methods work in practice. Finally, we have presented papers that employ these methods on real biological and physical problems, demonstrating their utility. We also consider some open research questions in the area of hybrid method development and the future directions for the field.

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 magnitude of the intrinsic rate constant: How deep can association reactions be in the diffusion limited regime?

Intrinsic and effective rate constants have an important role in the theory of diffusion-limited reactions. In a previous paper, we provide detailed microscopic expressions for these intrinsic rates [A. Vijaykumar, P. G. Bolhuis, and P. R. ten Wolde, Faraday Discuss. 195, 421 (2016)], which are usually considered as abstract quantities and assumed to be implicitly known. Using these microscopic expressions, we investigate how the rate of association depends on the strength and the range of the isotropic potential and the strength of the non-specific attraction in case of the anisotropic potential. In addition, we determine the location of the interface where these expressions become valid for anisotropic potentials. In particular, by investigating the particles' orientational distributions, we verify whether the interface at which these distributions become isotropic agrees with the interface predicted by the effective association rate constant. Finally, we discuss how large the intrinsic association rate can become, and what are the consequences for the existence of the diffusion limited regime.

Multiscale Modeling of Diffusion in a Crowded Environment

We present a multiscale approach to model diffusion in a crowded environment and its effect on the reaction rates. Diffusion in biological systems is often modeled by a discrete space jump process in order to capture the inherent noise of biological systems, which becomes important in the low copy number regime. To model diffusion in the crowded cell environment efficiently, we compute the jump rates in this mesoscopic model from local first exit times, which account for the microscopic positions of the crowding molecules, while the diffusing molecules jump on a coarser Cartesian grid. We then extract a macroscopic description from the resulting jump rates, where the excluded volume effect is modeled by a diffusion equation with space-dependent diffusion coefficient. The crowding molecules can be of arbitrary shape and size, and numerical experiments demonstrate that those factors together with the size of the diffusing molecule play a crucial role on the magnitude of the decrease in diffusive motion. When correcting the reaction rates for the altered diffusion we can show that molecular crowding either enhances or inhibits chemical reactions depending on local fluctuations of the obstacle density.