Simulation tools for particle-based reaction-diffusion dynamics in continuous space

Interplay between stochastic enzyme activity and microtubule stability drives detyrosination enrichment on microtubule subsets

Microtubules in cells consist of functionally diverse subpopulations carrying distinct post-translational modifications (PTMs). Akin to the histone code, the tubulin code regulates a myriad of microtubule functions, ranging from intracellular transport to chromosome segregation. However, how individual PTMs only occur on subsets of microtubules to contribute to microtubule specialization is not well understood. In particular, microtubule detyrosination, the removal of the C-terminal tyrosine on α-tubulin subunits, marks the stable population of microtubules and modifies how microtubules interact with other microtubule-associated proteins to regulate a wide range of cellular processes. Previously, we found that in certain cell types, only ∼30% of microtubules are highly enriched with the detyrosination mark and that detyrosination spans most of the length of a microtubule, often adjacent to a completely tyrosinated microtubule. How the activity of a cytosolic detyrosinase, vasohibin (VASH), leads to only a small subpopulation of highly detyrosinated microtubules is unclear. Here, using quantitative super-resolution microscopy, we visualized nascent microtubule detyrosination events in cells consisting of 1-3 detyrosinated α-tubulin subunits after nocodazole washout. Microtubule detyrosination accumulates slowly and in a dispersed pattern across the microtubule length. By visualizing single molecules of VASH in live cells, we found that VASH engages with microtubules stochastically on a short timescale, suggesting limited removal of tyrosine per interaction, consistent with the super-resolution results. Combining these quantitative imaging results with simulations incorporating parameters from our experiments, we provide evidence for a stochastic model for cells to establish a subset of detyrosinated microtubules via a detyrosination-stabilization feedback mechanism.

Rapid and Efficient Computation of Cell Paths During Ultrasonic Focusing

This biophysical analysis explores the first-principles physics of movement of white blood cell sized particles, suspended in an aqueous fluid and experiencing progressive or standing waves of acoustic pressure. In many current applications the cells are gradually nudged or herded toward the nodes of the standing wave, providing a degree of acoustic focusing and concentration of the cells in layers perpendicular to the direction of sound propagation. Here the underlying biomechanics of this phenomenon are analyzed specifically for the viscous regime of water and for small diameter microscopic spheroids such as living cells. The resulting mathematical model leads to a single algebraic expression for the creep or drift velocity as a function of sound frequency, amplitude, wavelength, fluid viscosity, boundary dimensions, and boundary reflectivity. This expression can be integrated numerically by a simple and fast computer algorithm to demonstrate net movement of particles as a function of time, providing a guide to optimization in a variety of emerging applications of ultrasonic cell focusing.

Implications of different membrane compartmentalization models in particle-based in silico studies

Studying membrane dynamics is important to understand the cellular response to environmental stimuli. A decisive spatial characteristic of the plasma membrane is its compartmental structure created by the actin-based membrane-skeleton (fences) and anchored transmembrane proteins (pickets). Particle-based reaction-diffusion simulation of the membrane offers a suitable temporal and spatial resolution to analyse its spatially heterogeneous and stochastic dynamics. Fences have been modelled via hop probabilities, potentials or explicit picket fences. Our study analyses the different approaches' constraints and their impact on simulation results and performance. Each of the methods comes with its own constraints; the picket fences require small timesteps, potential fences might induce a bias in diffusion in crowded systems, and probabilistic fences, in addition to carefully scaling the probability with the timesteps, induce higher computational costs for each propagation step.

Protein drift-diffusion dynamics and phase separation in curved cell membranes and dendritic spines: Hybrid discrete-continuum methods

We develop methods for investigating protein drift-diffusion dynamics in heterogeneous cell membranes and the roles played by geometry, diffusion, chemical kinetics, and phase separation. Our hybrid stochastic numerical methods combine discrete particle descriptions with continuum-level models for tracking the individual protein drift-diffusion dynamics when coupled to continuum fields. We show how our approaches can be used to investigate phenomena motivated by protein kinetics within dendritic spines. The spine geometry is hypothesized to play an important biological role regulating synaptic strength, protein kinetics, and self-assembly of clusters. We perform simulation studies for model spine geometries varying the neck size to investigate how phase-separation and protein organization is influenced by different shapes. We also show how our methods can be used to study the roles of geometry in reaction-diffusion systems including Turing instabilities. Our methods provide general approaches for investigating protein kinetics and drift-diffusion dynamics within curved membrane structures.

Thermodynamics and Kinetics of Aggregation of Flexible Peripheral Membrane Proteins

Biomembrane remodeling is essential for cellular trafficking, with membrane-binding peripheral proteins playing a key role in it. Significant membrane remodeling as in endo- and exocytosis is often due to aggregates of many proteins with direct or membrane-mediated interactions. Understanding this process via computer simulations is extremely challenging: protein-membrane systems involve time and length scales that make atomistic simulations impractical, while most coarse-grained models fall short in resolving dynamics and physical effects of protein and membrane flexibility. Here, we develop a coarse-grained model of the bilayer membrane bestrewed with rotationally symmetric flexible proteins, parametrized to reflect local curvatures and lateral dynamics of proteins. We investigate the kinetics, equilibrium distributions, and the free energy landscape governing the formation and breakup of protein clusters on the surface of the membrane. We demonstrate how the flexibility of the proteins as well as their surface concentration play deciding roles in highly selective macroscopic aggregation behavior.

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.

Hydrodynamic coupling for particle-based solvent-free membrane models

The great challenge with biological membrane systems is the wide range of scales involved, from nanometers and picoseconds for individual lipids to the micrometers and beyond millisecond for cellular signaling processes. While solvent-free coarse-grained membrane models are convenient for large-scale simulations and promising to provide insight into slow processes involving membranes, these models usually have unrealistic kinetics. One major obstacle is the lack of an equally convenient way of introducing hydrodynamic coupling without significantly increasing the computational cost of the model. To address this, we introduce a framework based on anisotropic Langevin dynamics, for which major in-plane and out-of-plane hydrodynamic effects are modeled via friction and diffusion tensors from analytical or semi-analytical solutions to Stokes hydrodynamic equations. Using this framework, in conjunction with our recently developed membrane model, we obtain accurate dispersion relations for planar membrane patches, both free-standing and in the vicinity of a wall. We briefly discuss how non-equilibrium dynamics is affected by hydrodynamic interactions. We also measure the surface viscosity of the model membrane and discuss the affecting dissipative mechanisms.

Electrostatic Reaction Inhibition in Nanoparticle Catalysis

Electrostatic reaction inhibition in heterogeneous catalysis emerges if charged reactants and products with similar charges are adsorbed on the catalyst and thus repel the approaching reactants. In this work, we study the effects of electrostatic inhibition on the reaction rate of unimolecular reactions catalyzed on the surface of a spherical model nanoparticle using particle-based reaction-diffusion simulations. Moreover, we derive closed rate equations based on an approximate Debye-Smoluchowski rate theory, valid for diffusion-controlled reactions, and a modified Langmuir adsorption isotherm, relevant for reaction-controlled reactions, to account for electrostatic inhibition in the Debye-Hückel limit. We study the kinetics of reactions ranging from low to high adsorptions on the nanoparticle surface and from the surface- to diffusion-controlled limits for charge valencies 1 and 2. In the diffusion-controlled limit, electrostatic inhibition drastically slows down the reactions for strong adsorption and low ionic concentration, which is well described by our theory. In particular, the rate decreases with adsorption affinity because, in this case, the inhibiting products are generated at a high rate. In the (slow) reaction-controlled limit, the effect of electrostatic inhibition is much weaker, as semiquantitatively reproduced by our electrostatic-modified Langmuir theory. We finally propose and verify a simple interpolation formula that describes electrostatic inhibition for all reaction speeds ("diffusion-influenced" reactions) in general.

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.

Multiscale simulation approaches to modeling drug-protein binding

Simulations can provide detailed insight into the molecular processes involved in drug action, such as protein-ligand binding, and can therefore be a valuable tool for drug design and development. Processes with a large range of length and timescales may be involved, and understanding these different scales typically requires different types of simulation methodology. Ideally, simulations should be able to connect across scales, to analyze and predict how changes at one scale can influence another. Multiscale simulation methods, which combine different levels of treatment, are an emerging frontier with great potential in this area. Here we review multiscale frameworks of various types, and selected applications to biomolecular systems with a focus on drug-ligand binding.

pSpatiocyte: a high-performance simulator for intracellular reaction-diffusion systems

Background

Studies using quantitative experimental methods have shown that intracellular spatial distribution of molecules plays a central role in many cellular systems. Spatially resolved computer simulations can integrate quantitative data from these experiments to construct physically accurate models of the systems. Although computationally expensive, microscopic resolution reaction-diffusion simulators, such as Spatiocyte can directly capture intracellular effects comprising diffusion-limited reactions and volume exclusion from crowded molecules by explicitly representing individual diffusing molecules in space. To alleviate the steep computational cost typically associated with the simulation of large or crowded intracellular compartments, we present a parallelized Spatiocyte method called pSpatiocyte.

Results

The new high-performance method employs unique parallelization schemes on hexagonal close-packed (HCP) lattice to efficiently exploit the resources of common workstations and large distributed memory parallel computers. We introduce a coordinate system for fast accesses to HCP lattice voxels, a parallelized event scheduler, a parallelized Gillespie’s direct-method for unimolecular reactions, and a parallelized event for diffusion and bimolecular reaction processes. We verified the correctness of pSpatiocyte reaction and diffusion processes by comparison to theory. To evaluate the performance of pSpatiocyte, we performed a series of parallelized diffusion runs on the RIKEN K computer. In the case of fine lattice discretization with low voxel occupancy, pSpatiocyte exhibited 74% parallel efficiency and achieved a speedup of 7686 times with 663552 cores compared to the runtime with 64 cores. In the weak scaling performance, pSpatiocyte obtained efficiencies of at least 60% with up to 663552 cores. When executing the Michaelis-Menten benchmark model on an eight-core workstation, pSpatiocyte required 45- and 55-fold shorter runtimes than Smoldyn and the parallel version of ReaDDy, respectively. As a high-performance application example, we study the dual phosphorylation-dephosphorylation cycle of the MAPK system, a typical reaction network motif in cell signaling pathways.

Conclusions

pSpatiocyte demonstrates good accuracies, fast runtimes and a significant performance advantage over well-known microscopic particle methods in large-scale simulations of intracellular reaction-diffusion systems. The source code of pSpatiocyte is available at https://spatiocyte.org.

Computational resources and strategies to assess single-molecule dynamics of the translation process in S. cerevisiae

This work provides a systematic and comprehensive overview of available resources for the molecular-scale modelling of the translation process through agent-based modelling. The case study is the translation in Saccharomyces cerevisiae, one of the most studied yeasts. The data curation workflow encompassed structural information about the yeast (i.e. the simulation environment), and the proteins, ribonucleic acids and other types of molecules involved in the process (i.e. the agents). Moreover, it covers the main process events, such as diffusion (i.e. motion of molecules in the environment) and collision efficiency (i.e. interaction between molecules). Data previously determined by wet-lab techniques were preferred, resorting to computational predictions/extrapolations only when strictly necessary. The computational modelling of the translation processes is of added industrial interest, since it may bring forward knowledge on how to control such phenomena and enhance the production of proteins of interest in a faster and more efficient manner.

Potential based, spatial simulation of dynamically nested particles

BACKGROUND

To study cell biological phenomena which depend on diffusion, active transport processes, or the locations of species, modeling and simulation studies need to take space into account. To describe the system as a collection of discrete objects moving and interacting in continuous space, various particle-based reaction diffusion simulators for cell-biological system have been developed. So far the focus has been on particles as solid spheres or points. However, spatial dynamics might happen at different organizational levels, such as proteins, vesicles or cells with interrelated dynamics which requires spatial approaches that take this multi-levelness of cell biological systems into account.

RESULTS

Based on the perception of particles forming hollow spheres, ML-Force contributes to the family of particle-based simulation approaches: in addition to excluded volumes and forces, it also supports compartmental dynamics and relating dynamics between different organizational levels explicitly. Thereby, compartmental dynamics, e.g., particles entering and leaving other particles, and bimolecular reactions are modeled using pair-wise potentials (forces) and the Langevin equation. In addition, forces that act independently of other particles can be applied to direct the movement of particles. Attributes and the possibility to define arbitrary functions on particles, their attributes and content, to determine the results and kinetics of reactions add to the expressiveness of ML-Force. Its implementation comprises a rudimentary rule-based embedded domain-specific modeling language for specifying models and a simulator for executing models continuously. Applications inspired by cell biological models from literature, such as vesicle transport or yeast growth, show the value of the realized features. They facilitate capturing more complex spatial dynamics, such as the fission of compartments or the directed movement of particles, and enable the integration of non-spatial intra-compartmental dynamics as stochastic events.

CONCLUSIONS

By handling all dynamics based on potentials (forces) and the Langevin equation, compartmental dynamics, such as dynamic nesting, fusion and fission of compartmental structures are handled continuously and are seamlessly integrated with traditional particle-based reaction-diffusion dynamics within the cell. Thereby, attributes and arbitrary functions allow to flexibly describe diverse spatial phenomena, and relate dynamics across organizational levels. Also they prove crucial in modeling intra-cellular or intra-compartmental dynamics in a non-spatial manner, and, thus, to abstract from spatial dynamics, on demand which increases the range of multi-compartmental processes that can be captured.

Diffusion-influenced reaction rates in the presence of pair interactions

The kinetics of bimolecular reactions in solution depends, among other factors, on intermolecular forces such as steric repulsion or electrostatic interaction. Microscopically, a pair of molecules first has to meet by diffusion before the reaction can take place. In this work, we establish an extension of Doi's volume reaction model to molecules interacting via pair potentials, which is a key ingredient for interacting-particle-based reaction-diffusion (iPRD) simulations. As a central result, we relate model parameters and macroscopic reaction rate constants in this situation. We solve the corresponding reaction-diffusion equation in the steady state and derive semi-analytical expressions for the reaction rate constant and the local concentration profiles. Our results apply to the full spectrum from well-mixed to diffusion-limited kinetics. For limiting cases, we give explicit formulas, and we provide a computationally inexpensive numerical scheme for the general case, including the intermediate, diffusion-influenced regime. The obtained rate constants decompose uniquely into encounter and formation rates, and we discuss the effect of the potential on both subprocesses, exemplified for a soft harmonic repulsion and a Lennard-Jones potential. The analysis is complemented by extensive stochastic iPRD simulations, and we find excellent agreement with the theoretical predictions.

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.

Surface reaction-diffusion kinetics on lattice at the microscopic scale

Microscopic models of reaction-diffusion processes on the cell membrane can link local spatiotemporal effects to macroscopic self-organized patterns often observed on the membrane. Simulation schemes based on the microscopic lattice method (MLM) can model these processes at the microscopic scale by tracking individual molecules, represented as hard spheres, on fine lattice voxels. Although MLM is simple to implement and is generally less computationally demanding than off-lattice approaches, its accuracy and consistency in modeling surface reactions have not been fully verified. Using the Spatiocyte scheme, we study the accuracy of MLM in diffusion-influenced surface reactions. We derive the lattice-based bimolecular association rates for two-dimensional (2D) surface-surface reaction and one-dimensional (1D) volume-surface adsorption according to the Smoluchowski-Collins-Kimball model and random walk theory. We match the time-dependent rates on lattice with off-lattice counterparts to obtain the correct expressions for MLM parameters in terms of physical constants. The expressions indicate that the voxel size needs to be at least 0.6% larger than the molecule to accurately simulate surface reactions on triangular lattice. On square lattice, the minimum voxel size should be even larger, at 5%. We also demonstrate the ability of MLM-based schemes such as Spatiocyte to simulate a reaction-diffusion model that involves all dimensions: three-dimensional (3D) diffusion in the cytoplasm, 2D diffusion on the cell membrane, and 1D cytoplasm-membrane adsorption. With the model, we examine the contribution of the 2D reaction pathway to the overall reaction rate at different reactant diffusivity, reactivity, and concentrations.

Dynamics of proteins in solution

The dynamics of proteins in solution includes a variety of processes, such as backbone and side-chain fluctuations, interdomain motions, as well as global rotational and translational (i.e. center of mass) diffusion. Since protein dynamics is related to protein function and essential transport processes, a detailed mechanistic understanding and monitoring of protein dynamics in solution is highly desirable. The hierarchical character of protein dynamics requires experimental tools addressing a broad range of time- and length scales. We discuss how different techniques contribute to a comprehensive picture of protein dynamics, and focus in particular on results from neutron spectroscopy. We outline the underlying principles and review available instrumentation as well as related analysis frameworks.

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.

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.

A spatio-temporal model reveals self-limiting FcɛRI cross-linking by multivalent antigens

Aggregation of cell surface receptor proteins by multivalent antigens is an essential early step for immune cell signalling. A number of experimental and modelling studies in the past have investigated multivalent ligand-mediated aggregation of IgE receptors (FcɛRI) in the plasma membrane of mast cells. However, understanding of the mechanisms of FcɛRI aggregation remains incomplete. Experimental reports indicate that FcɛRI forms relatively small and finite-sized clusters when stimulated by a multivalent ligand. By contrast, modelling studies have shown that receptor cross-linking by a trivalent ligand may lead to the formation of large receptor superaggregates that may potentially give rise to hyperactive cellular responses. In this work, we have developed a Brownian dynamics-based spatio-temporal model to analyse FcɛRI aggregation by a trivalent antigen. Unlike the existing models, which implemented non-spatial simulation approaches, our model explicitly accounts for the coarse-grained site-specific features of the multivalent species (molecules and complexes). The model incorporates membrane diffusion, steric collisions and sub-nanometre-scale site-specific interaction of the time-evolving species of arbitrary structures. Using the model, we investigated temporal evolution of the species and their diffusivities. Consistent with a recent experimental report, our model predicted sharp decay in species mobility in the plasma membrane in response receptor cross-linking by a multivalent antigen. We show that, due to such decay in the species mobility, post-stimulation receptor aggregation may become self-limiting. Our analysis reveals a potential regulatory mechanism suppressing hyperactivation of immune cells in response to multivalent antigens.

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.

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.

Particle-based membrane model for mesoscopic simulation of cellular dynamics

We present a simple and computationally efficient coarse-grained and solvent-free model for simulating lipid bilayer membranes. In order to be used in concert with particle-based reaction-diffusion simulations, the model is purely based on interacting and reacting particles, each representing a coarse patch of a lipid monolayer. Particle interactions include nearest-neighbor bond-stretching and angle-bending and are parameterized so as to reproduce the local membrane mechanics given by the Helfrich energy density over a range of relevant curvatures. In-plane fluidity is implemented with Monte Carlo bond-flipping moves. The physical accuracy of the model is verified by five tests: (i) Power spectrum analysis of equilibrium thermal undulations is used to verify that the particle-based representation correctly captures the dynamics predicted by the continuum model of fluid membranes. (ii) It is verified that the input bending stiffness, against which the potential parameters are optimized, is accurately recovered. (iii) Isothermal area compressibility modulus of the membrane is calculated and is shown to be tunable to reproduce available values for different lipid bilayers, independent of the bending rigidity. (iv) Simulation of two-dimensional shear flow under a gravity force is employed to measure the effective in-plane viscosity of the membrane model and show the possibility of modeling membranes with specified viscosities. (v) Interaction of the bilayer membrane with a spherical nanoparticle is modeled as a test case for large membrane deformations and budding involved in cellular processes such as endocytosis. The results are shown to coincide well with the predicted behavior of continuum models, and the membrane model successfully mimics the expected budding behavior. We expect our model to be of high practical usability for ultra coarse-grained molecular dynamics or particle-based reaction-diffusion simulations of biological systems.

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.

Challenges in Reproducibility, Replicability, and Comparability of Computational Models and Tools for Neuronal and Glial Networks, Cells, and Subcellular Structures

The possibility to replicate and reproduce published research results is one of the biggest challenges in all areas of science. In computational neuroscience, there are thousands of models available. However, it is rarely possible to reimplement the models based on the information in the original publication, let alone rerun the models just because the model implementations have not been made publicly available. We evaluate and discuss the comparability of a versatile choice of simulation tools: tools for biochemical reactions and spiking neuronal networks, and relatively new tools for growth in cell cultures. The replicability and reproducibility issues are considered for computational models that are equally diverse, including the models for intracellular signal transduction of neurons and glial cells, in addition to single glial cells, neuron-glia interactions, and selected examples of spiking neuronal networks. We also address the comparability of the simulation results with one another to comprehend if the studied models can be used to answer similar research questions. In addition to presenting the challenges in reproducibility and replicability of published results in computational neuroscience, we highlight the need for developing recommendations and good practices for publishing simulation tools and computational models. Model validation and flexible model description must be an integral part of the tool used to simulate and develop computational models. Constant improvement on experimental techniques and recording protocols leads to increasing knowledge about the biophysical mechanisms in neural systems. This poses new challenges for computational neuroscience: extended or completely new computational methods and models may be required. Careful evaluation and categorization of the existing models and tools provide a foundation for these future needs, for constructing multiscale models or extending the models to incorporate additional or more detailed biophysical mechanisms. Improving the quality of publications in computational neuroscience, enabling progressive building of advanced computational models and tools, can be achieved only through adopting publishing standards which underline replicability and reproducibility of research results.

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.

Fluctuating Finite Element Analysis (FFEA): A continuum mechanics software tool for mesoscale simulation of biomolecules

Fluctuating Finite Element Analysis (FFEA) is a software package designed to perform continuum mechanics simulations of proteins and other globular macromolecules. It combines conventional finite element methods with stochastic thermal noise, and is appropriate for simulations of large proteins and protein complexes at the mesoscale (length-scales in the range of 5 nm to 1 μm), where there is currently a paucity of modelling tools. It requires 3D volumetric information as input, which can be low resolution structural information such as cryo-electron tomography (cryo-ET) maps or much higher resolution atomistic co-ordinates from which volumetric information can be extracted. In this article we introduce our open source software package for performing FFEA simulations which we have released under a GPLv3 license. The software package includes a C ++ implementation of FFEA, together with tools to assist the user to set up the system from Electron Microscopy Data Bank (EMDB) or Protein Data Bank (PDB) data files. We also provide a PyMOL plugin to perform basic visualisation and additional Python tools for the analysis of FFEA simulation trajectories. This manuscript provides a basic background to the FFEA method, describing the implementation of the core mechanical model and how intermolecular interactions and the solvent environment are included within this framework. We provide prospective FFEA users with a practical overview of how to set up an FFEA simulation with reference to our publicly available online tutorials and manuals that accompany this first release of the package.

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.

Bridging intracellular scales by mechanistic computational models

The impact of intracellular spatial organization beyond classical compartments on processes such as cell signaling is increasingly recognized. A quantitative, mechanistic understanding of cellular systems therefore needs to account for different scales in at least three coordinates: time, molecular abundances, and space. Mechanistic mathematical models may span all these scales, but corresponding multi-scale models need to resolve mechanistic details on small scales while maintaining computational tractability for larger ones. This typically results in models that combine different levels of description: from a microscopic representation of chemical reactions up to continuum dynamics in space and time. We highlight recent progress in bridging these model classes and outline current challenges in multi-scale models such as active transport and dynamic geometries.

Multi-Algorithm Particle Simulations with Spatiocyte

As quantitative biologists get more measurements of spatially regulated systems such as cell division and polarization, simulation of reaction and diffusion of proteins using the data is becoming increasingly relevant to uncover the mechanisms underlying the systems. Spatiocyte is a lattice-based stochastic particle simulator for biochemical reaction and diffusion processes. Simulations can be performed at single molecule and compartment spatial scales simultaneously. Molecules can diffuse and react in 1D (filament), 2D (membrane), and 3D (cytosol) compartments. The implications of crowded regions in the cell can be investigated because each diffusing molecule has spatial dimensions. Spatiocyte adopts multi-algorithm and multi-timescale frameworks to simulate models that simultaneously employ deterministic, stochastic, and particle reaction-diffusion algorithms. Comparison of light microscopy images to simulation snapshots is supported by Spatiocyte microscopy visualization and molecule tagging features. Spatiocyte is open-source software and is freely available at http://spatiocyte.org .

Smoldyn: particle-based simulation with rule-based modeling, improved molecular interaction and a library interface

Motivation

Smoldyn is a spatial and stochastic biochemical simulator. It treats each molecule of interest as an individual particle in continuous space, simulating molecular diffusion, molecule-membrane interactions and chemical reactions, all with good accuracy. This article presents several new features.

Results

Smoldyn now supports two types of rule-based modeling. These are a wildcard method, which is very convenient, and the BioNetGen package with extensions for spatial simulation, which is better for complicated models. Smoldyn also includes new algorithms for simulating the diffusion of surface-bound molecules and molecules with excluded volume. Both are exact in the limit of short time steps and reasonably good with longer steps. In addition, Smoldyn supports single-molecule tracking simulations. Finally, the Smoldyn source code can be accessed through a C/C ++ language library interface.

Availability and Implementation

Smoldyn software, documentation, code, and examples are at http://www.smoldyn.org .

Contact

steven.s.andrews@gmail.com.

Molecular simulations of cellular processes

It is, nowadays, possible to simulate biological processes in conditions that mimic the different cellular compartments. Several groups have performed these calculations using molecular models that vary in performance and accuracy. In many cases, the atomistic degrees of freedom have been eliminated, sacrificing both structural complexity and chemical specificity to be able to explore slow processes. In this review, we will discuss the insights gained from computer simulations on macromolecule diffusion, nuclear body formation, and processes involving the genetic material inside cell-mimicking spaces. We will also discuss the challenges to generate new models suitable for the simulations of biological processes on a cell scale and for cell-cycle-long times, including non-equilibrium events such as the co-translational folding, misfolding, and aggregation of proteins. A prominent role will be played by the wise choice of the structural simplifications and, simultaneously, of a relatively complex energetic description. These challenging tasks will rely on the integration of experimental and computational methods, achieved through the application of efficient algorithms. Graphical abstract.

Protein-peptide association kinetics beyond the seconds timescale from atomistic simulations

Understanding and control of structures and rates involved in protein ligand binding are essential for drug design. Unfortunately, atomistic molecular dynamics (MD) simulations cannot directly sample the excessively long residence and rearrangement times of tightly binding complexes. Here we exploit the recently developed multi-ensemble Markov model framework to compute full protein-peptide kinetics of the oncoprotein fragment 25-109Mdm2 and the nano-molar inhibitor peptide PMI. Using this system, we report, for the first time, direct estimates of kinetics beyond the seconds timescale using simulations of an all-atom MD model, with high accuracy and precision. These results only require explicit simulations on the sub-milliseconds timescale and are tested against existing mutagenesis data and our own experimental measurements of the dissociation and association rates. The full kinetic model reveals an overall downhill but rugged binding funnel with multiple pathways. The overall strong binding arises from a variety of conformations with different hydrophobic contact surfaces that interconvert on the milliseconds timescale.

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.

Lipid-mediated PX-BAR domain recruitment couples local membrane constriction to endocytic vesicle fission

Clathrin-mediated endocytosis (CME) involves membrane-associated scaffolds of the bin-amphiphysin-rvs (BAR) domain protein family as well as the GTPase dynamin, and is accompanied and perhaps triggered by changes in local lipid composition. How protein recruitment, scaffold assembly and membrane deformation is spatiotemporally controlled and coupled to fission is poorly understood. We show by computational modelling and super-resolution imaging that phosphatidylinositol 3,4-bisphosphate [PI(3,4)P₂] synthesis within the clathrin-coated area of endocytic intermediates triggers selective recruitment of the PX-BAR domain protein SNX9, as a result of complex interactions of endocytic proteins competing for phospholipids. The specific architecture induces positioning of SNX9 at the invagination neck where its self-assembly regulates membrane constriction, thereby providing a template for dynamin fission. These data explain how lipid conversion at endocytic pits couples local membrane constriction to fission. Our work demonstrates how computational modelling and super-resolution imaging can be combined to unravel function and mechanisms of complex cellular processes.

The spatiotemporal regulation of membrane scaffolds recruitment and coupling between membrane deformation and fission in endocytosis are unclear. Here the authors show that lipid conversion at endocytic pits recruits SNX9, which couples local membrane constriction to fission in endocytosis.

Excluded volume effects in on- and off-lattice reaction-diffusion models

Mathematical models are important tools to study the excluded volume effects on reaction-diffusion systems, which are known to play an important role inside living cells. Detailed microscopic simulations with off-lattice Brownian dynamics become computationally expensive in crowded environments. In this study, the authors therefore investigate to which extent on-lattice approximations, the so-called cellular automata models, can be used to simulate reactions and diffusion in the presence of crowding molecules. They show that the diffusion is most severely slowed down in the off-lattice model, since randomly distributed obstacles effectively exclude more volume than those ordered on an artificial grid. Crowded reaction rates can be both increased and decreased by the grid structure and it proves important to model the molecules with realistic sizes when excluded volume is taken into account. The grid artefacts increase with increasing crowder density and they conclude that the computationally more efficient on-lattice simulations are accurate approximations only for low crowder densities.

SpringSaLaD: A Spatial, Particle-Based Biochemical Simulation Platform with Excluded Volume

We introduce Springs, Sites, and Langevin Dynamics (SpringSaLaD), a comprehensive software platform for spatial, stochastic, particle-based modeling of biochemical systems. SpringSaLaD models biomolecules in a coarse-grained manner as a group of linked spherical sites with excluded volume. This mesoscopic approach bridges the gap between highly detailed molecular dynamics simulations and the various methods used to study network kinetics and diffusion at the cellular level. SpringSaLaD is a standalone tool that supports model building, simulation, visualization, and data analysis, all through a user-friendly graphical user interface that should make it more accessible than tools built into more comprehensive molecular dynamics infrastructures. Importantly, for bimolecular reactions we derive an exact expression relating the macroscopic on-rate to the various microscopic parameters with the inclusion of excluded volume; this makes SpringSaLaD more accurate than other tools, which rely on approximate relationships between these parameters.

Reaction-diffusion basis of retroviral infectivity

Retrovirus particle (virion) infectivity requires diffusion and clustering of multiple transmembrane envelope proteins (Env₃) on the virion exterior, yet is triggered by protease-dependent degradation of a partially occluding, membrane-bound Gag polyprotein lattice on the virion interior. The physical mechanism underlying such coupling is unclear and only indirectly accessible via experiment. Modelling stands to provide insight but the required spatio-temporal range far exceeds current accessibility by all-atom or even coarse-grained molecular dynamics simulations. Nor do such approaches account for chemical reactions, while conversely, reaction kinetics approaches handle neither diffusion nor clustering. Here, a recently developed multiscale approach is considered that applies an ultra-coarse-graining scheme to treat entire proteins at near-single particle resolution, but which also couples chemical reactions with diffusion and interactions. A model is developed of Env₃ molecules embedded in a truncated Gag lattice composed of membrane-bound matrix proteins linked to capsid subunits, with freely diffusing protease molecules. Simulations suggest that in the presence of Gag but in the absence of lateral lattice-forming interactions, Env₃ diffuses comparably to Gag-absent Env₃ Initial immobility of Env₃ is conferred through lateral caging by matrix trimers vertically coupled to the underlying hexameric capsid layer. Gag cleavage by protease vertically decouples the matrix and capsid layers, induces both matrix and Env₃ diffusion, and permits Env₃ clustering. Spreading across the entire membrane surface reduces crowding, in turn, enhancing the effect and promoting infectivity.This article is part of the themed issue 'Multiscale modelling at the physics-chemistry-biology interface'.

Collective ratchet effects and reversals for active matter particles on quasi-one-dimensional asymmetric substrates

Using computer simulations, we study a two-dimensional system of sterically interacting self-mobile run-and-tumble disk-shaped particles with an underlying periodic quasi-one-dimensional asymmetric substrate, and show that a rich variety of collective active ratchet behaviors arise as a function of particle density, activity, substrate period, and the maximum force exerted by the substrate. The net dc drift, or ratchet transport flux, is nonmonotonic since it increases with increased activity but is diminished by the onset of self-clustering of the active particles. Increasing the particle density decreases the ratchet transport flux for shallow substrates but increases the ratchet transport flux for deep substrates due to collective hopping events. At the highest particle densities, the ratchet motion is destroyed by a self-jamming effect. We show that it is possible to realize reversals of the direction of the net dc drift in the deep substrate limit when multiple rows of active particles can be confined in each substrate minimum, permitting emergent particle-like excitations to appear that experience an inverted effective substrate potential. We map out a phase diagram of the forward and reverse ratchet effects as a function of the particle density, activity, and substrate properties.

Computational support for a scaffolding mechanism of centriole assembly

Centrioles are essential for forming cilia, flagella and centrosomes. Successful centriole assembly requires proteins of the SAS-6 family, which can form oligomeric ring structures with ninefold symmetry in vitro. While important progress has been made in understanding SAS-6 protein biophysics, the mechanisms enabling ring formation in vivo remain elusive. Likewise, the mechanisms by which a nascent centriole forms near-orthogonal to an existing one are not known. Here, we investigate possible mechanisms of centriole assembly using coarse-grained Brownian dynamics computer simulations in combination with a rate equation approach. Our results suggest that without any external factors, strong stabilization associated with ring closure would be needed to enable efficient ring formation. Strikingly, our simulations reveal that a scaffold-assisted assembly mechanism can trigger robust ring formation owing to local cooperativity, and that this mechanism can also impart orthogonalilty to centriole assembly. Overall, our findings provide novel insights into the organizing principles governing the assembly of this important organelle.

Non-Brownian diffusion in lipid membranes: Experiments and simulations

The dynamics of constituents and the surface response of cellular membranes-also in connection to the binding of various particles and macromolecules to the membrane-are still a matter of controversy in the membrane biophysics community, particularly with respect to crowded membranes of living biological cells. We here put into perspective recent single particle tracking experiments in the plasma membranes of living cells and supercomputing studies of lipid bilayer model membranes with and without protein crowding. Special emphasis is put on the observation of anomalous, non-Brownian diffusion of both lipid molecules and proteins embedded in the lipid bilayer. While single component, pure lipid bilayers in simulations exhibit only transient anomalous diffusion of lipid molecules on nanosecond time scales, the persistence of anomalous diffusion becomes significantly longer ranged on the addition of disorder-through the addition of cholesterol or proteins-and on passing of the membrane lipids to the gel phase. Concurrently, experiments demonstrate the anomalous diffusion of membrane embedded proteins up to macroscopic time scales in the minute time range. Particular emphasis will be put on the physical character of the anomalous diffusion, in particular, the occurrence of ageing observed in the experiments-the effective diffusivity of the measured particles is a decreasing function of time. Moreover, we present results for the time dependent local scaling exponent of the mean squared displacement of the monitored particles. Recent results finding deviations from the commonly assumed Gaussian diffusion patterns in protein crowded membranes are reported. The properties of the displacement autocorrelation function of the lipid molecules are discussed in the light of their appropriate physical anomalous diffusion models, both for non-crowded and crowded membranes. In the last part of this review we address the upcoming field of membrane distortion by elongated membrane-binding particles. We discuss how membrane compartmentalisation and the particle-membrane binding energy may impact the dynamics and response of lipid membranes. This article is part of a Special Issue entitled: Biosimulations edited by Ilpo Vattulainen and Tomasz Róg.

ANALYSIS AND DESIGN OF JUMP COEFFICIENTS IN DISCRETE STOCHASTIC DIFFUSION MODELS

In computational systems biology, the mesoscopic model of reaction-diffusion kinetics is described by a continuous time, discrete space Markov process. To simulate diffusion stochastically, the jump coefficients are obtained by a discretization of the diffusion equation. Using unstructured meshes to represent complicated geometries may lead to negative coefficients when using piecewise linear finite elements. Several methods have been proposed to modify the coefficients to enforce the nonnegativity needed in the stochastic setting. In this paper, we present a method to quantify the error introduced by that change. We interpret the modified discretization matrix as the exact finite element discretization of a perturbed equation. The forward error, the error between the analytical solutions to the original and the perturbed equations, is bounded by the backward error, the error between the diffusion of the two equations. We present a backward analysis algorithm to compute the diffusion coefficient from a given discretization matrix. The analysis suggests a new way of deriving nonnegative jump coefficients that minimizes the backward error. The theory is tested in numerical experiments indicating that the new method is superior and also minimizes the forward error.

MOLNs: A CLOUD PLATFORM FOR INTERACTIVE, REPRODUCIBLE, AND SCALABLE SPATIAL STOCHASTIC COMPUTATIONAL EXPERIMENTS IN SYSTEMS BIOLOGY USING PyURDME

Computational experiments using spatial stochastic simulations have led to important new biological insights, but they require specialized tools and a complex software stack, as well as large and scalable compute and data analysis resources due to the large computational cost associated with Monte Carlo computational workflows. The complexity of setting up and managing a large-scale distributed computation environment to support productive and reproducible modeling can be prohibitive for practitioners in systems biology. This results in a barrier to the adoption of spatial stochastic simulation tools, effectively limiting the type of biological questions addressed by quantitative modeling. In this paper, we present PyURDME, a new, user-friendly spatial modeling and simulation package, and MOLNs, a cloud computing appliance for distributed simulation of stochastic reaction-diffusion models. MOLNs is based on IPython and provides an interactive programming platform for development of sharable and reproducible distributed parallel computational experiments.

ReaDDyMM: Fast interacting particle reaction-diffusion simulations using graphical processing units

ReaDDy is a modular particle simulation package combining off-lattice reaction kinetics with arbitrary particle interaction forces. Here we present a graphical processing unit implementation of ReaDDy that employs the fast multiplatform molecular dynamics package OpenMM. A speedup of up to two orders of magnitude is demonstrated, giving us access to timescales of multiple seconds on single graphical processing units. This opens up the possibility of simulating cellular signal transduction events while resolving all protein copies.

Dynamical Organization of Syntaxin-1A at the Presynaptic Active Zone

Synaptic vesicle fusion is mediated by SNARE proteins forming in between synaptic vesicle (v-SNARE) and plasma membrane (t-SNARE), one of which is Syntaxin-1A. Although exocytosis mainly occurs at active zones, Syntaxin-1A appears to cover the entire neuronal membrane. By using STED super-resolution light microscopy and image analysis of Drosophila neuro-muscular junctions, we show that Syntaxin-1A clusters are more abundant and have an increased size at active zones. A computational particle-based model of syntaxin cluster formation and dynamics is developed. The model is parametrized to reproduce Syntaxin cluster-size distributions found by STED analysis, and successfully reproduces existing FRAP results. The model shows that the neuronal membrane is adjusted in a way to strike a balance between having most syntaxins stored in large clusters, while still keeping a mobile fraction of syntaxins free or in small clusters that can efficiently search the membrane or be traded between clusters. This balance is subtle and can be shifted toward almost no clustering and almost complete clustering by modifying the syntaxin interaction energy on the order of only 1 kBT. This capability appears to be exploited at active zones. The larger active-zone syntaxin clusters are more stable and provide regions of high docking and fusion capability, whereas the smaller clusters outside may serve as flexible reserve pool or sites of spontaneous ectopic release.