Computational methods for diffusion-influenced biochemical reactions

Particle-based MR modeling with diffusion, microstructure, and enzymatic reactions

PURPOSE

To develop a novel particle-based in silico MR model and demonstrate applications of this model to signal mechanisms which are affected by the spatial organization of particles, including metabolic reaction kinetics, microstructural effects on diffusion, and radiofrequency (RF) refocusing effects in gradient-echo sequences.

METHODS

The model was developed by integrating a forward solution of the Bloch equations with a Brownian dynamics simulator. Simulation configurations were then designed to model MR signal dynamics of interest, with a primary focus on hyperpolarized 13C MRI methods. Phantom scans and spectrophotometric assays were conducted to validate model results in vitro.

RESULTS

The model accurately reproduced the reaction kinetics of enzyme-mediated conversion of pyruvate to lactate. When varying proportions of restrictive structure were added to the reaction volume, nonlinear changes in the reaction rate measured in vitro were replicated in silico. Modeling of RF refocusing effects characterized the degree of diffusion-weighted contribution from preserved residual magnetization in nonspoiled gradient-echo sequences.

CONCLUSIONS

These results show accurate reproduction of a range of MR signal mechanisms, establishing the model's capability to investigate the multifactorial signal dynamics such as those underlying hyperpolarized 13C MRI data.

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

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

Stochastic self-tuning hybrid algorithm for reaction-diffusion systems

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

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.

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.

Diffusion-dynamics laws in stochastic reaction networks

Many biological activities are induced by cellular chemical reactions of diffusing reactants. The dynamics of such systems can be captured by stochastic reaction networks. A recent numerical study has shown that diffusion can significantly enhance the fluctuations in gene regulatory networks. However, the universal relation between diffusion and stochastic system dynamics remains veiled. Within the approximation of reaction-diffusion master equation (RDME), we find general relation that the steady-state distribution in complex balanced networks is diffusion-independent. Here, complex balance is the nonequilibrium generalization of detailed balance. We also find that for a diffusion-included network with a Poisson-like steady-state distribution, the diffusion can be ignored at steady state. We then derive a necessary and sufficient condition for networks holding such steady-state distributions. Moreover, we show that for linear reaction networks the RDME reduces to the chemical master equation, which implies that the stochastic dynamics of networks is unaffected by diffusion at any arbitrary time. Our findings shed light on the fundamental question of when diffusion can be neglected, or (if nonnegligible) its effects on the stochastic dynamics of the reaction network.

Stochastic Simulation of Pattern Formation in Growing Tissue: A Multilevel Approach

We take up the challenge of designing realistic computational models of large interacting cell populations. The goal is essentially to bring Gillespie's celebrated stochastic methodology to the level of an interacting population of cells. Specifically, we are interested in how the gold standard of single-cell computational modeling, here taken to be spatial stochastic reaction-diffusion models, may be efficiently coupled with a similar approach at the cell population level. Concretely, we target a recently proposed set of pathways for pattern formation involving Notch-Delta signaling mechanisms. These involve cell-to-cell communication as mediated both via direct membrane contact sites and via cellular protrusions. We explain how to simulate the process in growing tissue using a multilevel approach and we discuss implications for future development of the associated computational methods.

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.

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.

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.

Looking "Under the Hood" of Cellular Mechanotransduction with Computational Tools: A Systems Biomechanics Approach across Multiple Scales

Signal modulation has been developed in living cells throughout evolution to promote utilizing the same machinery for multiple cellular functions. Chemical and mechanical modules of signal transmission and transduction are interconnected and necessary for organ development and growth. However, due to the high complexity of the intercommunication of physical intracellular connections with biochemical pathways, there are many missing details in our overall understanding of mechanotransduction processes, i.e., the process by which mechanical signals are converted to biochemical cascades. Cell-matrix adhesions are mechanically coupled to the nucleus through the cytoskeleton. This modulated and tightly integrated network mediates the transmission of mechanochemical signals from the extracellular matrix to the nucleus. Various experimental and computational techniques have been utilized to understand the basic mechanisms of mechanotransduction, yet many aspects have remained elusive. Recently, in silico experiments have made important contributions to the field of mechanobiology. Herein, computational modeling efforts devoted to understanding integrin-mediated mechanotransduction pathways are reviewed, and an outlook is presented for future directions toward using suitable computational approaches and developing novel techniques for addressing important questions in the field of mechanotransduction.

Macromolecular Crowding Regulates the Gene Expression Profile by Limiting Diffusion

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

The cellular nucleus is packed with macromolecules such as DNAs and proteins, which leaves limited space for other molecules to move around. Recent experimental results by C. Tan et al. have shown that macromolecular crowding can regulate gene expression, resulting in a more homogenous cell population. We introduce a computational model to uncover the mechanism by which macromolecular crowding functions. Our results suggest that macromolecular crowding limits the diffusion of the transcription factors and attenuates the transcriptional bursting, which leads to a more homogenous cell population. Regulation of gene expression noise by macromolecules depends on the size of the crowders, i.e. larger macromolecules can reduce the noise more effectively than smaller macromolecules. We also demonstrate that local variation of chromatin density can affect the noise of gene expression. This shows the importance of the chromatin structure in gene expression regulation.

Quantitative and logic modelling of molecular and gene networks

Behaviours of complex biomolecular systems are often irreducible to the elementary properties of their individual components. Explanatory and predictive mathematical models are therefore useful for fully understanding and precisely engineering cellular functions. The development and analyses of these models require their adaptation to the problems that need to be solved and the type and amount of available genetic or molecular data. Quantitative and logic modelling are among the main methods currently used to model molecular and gene networks. Each approach comes with inherent advantages and weaknesses. Recent developments show that hybrid approaches will become essential for further progress in synthetic biology and in the development of virtual organisms.

Methodologies for the modeling and simulation of biochemical networks, illustrated for signal transduction pathways: a primer

Biochemical networks depict the chemical interactions that take place among elements of living cells. They aim to elucidate how cellular behavior and functional properties of the cell emerge from the relationships between its components, i.e. molecules. Biochemical networks are largely characterized by dynamic behavior, and exhibit high degrees of complexity. Hence, the interest in such networks is growing and they have been the target of several recent modeling efforts. Signal transduction pathways (STPs) constitute a class of biochemical networks that receive, process, and respond to stimuli from the environment, as well as stimuli that are internal to the organism. An STP consists of a chain of intracellular signaling processes that ultimately result in generating different cellular responses. This primer presents the methodologies used for the modeling and simulation of biochemical networks, illustrated for STPs. These methodologies range from qualitative to quantitative, and include structural as well as dynamic analysis techniques. We describe the different methodologies, outline their underlying assumptions, and provide an assessment of their advantages and disadvantages. Moreover, publicly and/or commercially available implementations of these methodologies are listed as appropriate. In particular, this primer aims to provide a clear introduction and comprehensive coverage of biochemical modeling and simulation methodologies for the non-expert, with specific focus on relevant literature of STPs.

Immobilized IL-8 Triggers Phagocytosis and Dynamic Changes in Membrane Microtopology in Human Neutrophils

The interaction of leukocytes with surface bound ligands can be limited by the location of the molecules relative to the surface topology of the cell. In this report, we examine the dynamic response of neutrophils to IL-8-fractalkine chimera immobilized on bead surfaces, taking into account changes in receptor occupancy resulting from changes in surface topography. As a readout for receptor signaling, we observe the dynamics of calcium release in neutrophils following contact with the IL-8 coated surface. After a delay that depended on the initial area of contact and the surface density of IL-8, the cell began to phagocytose the IL-8 coated bead. This appeared to be a pre-requisite for release of calcium, which typically followed shortly after the initiation of phagocytosis. In separate experiments, effective kinetic coefficients for the formation of bonds between immobilized IL-8 and receptors on the cell surface were determined. Using these coefficients, we were able to estimate the number of bound receptors in the nascent contact zone. Kinetic modeling of the signaling response predicted that cell spreading and a concomitant increase in the density of occupied receptors would be required for the experimentally observed calcium dynamics. Postulating that there is an increase in receptor occupancy resulting from smoothing of the cell surface as it is stretched over the bead enabled us to obtain model predictions consistent with experimental observations. This study reveals the likely importance of membrane microtopology as a rate-limiting property and potential means of regulation of cell responses stimulated by two-dimensional surface interactions.

FOUNDATIONS FOR MODELING THE DYNAMICS OF GENE REGULATORY NETWORKS: A MULTILEVEL-PERSPECTIVE REVIEW

A promising alternative for unraveling the principles under which the dynamic interactions among genes lead to cellular phenotypes relies on mathematical and computational models at different levels of abstraction, from the molecular level of protein-DNA interactions to the system level of functional relationships among genes. This review article presents, under a bottom–up perspective, a hierarchy of approaches to modeling gene regulatory network dynamics, from microscopic descriptions at the single-molecule level in the spatial context of an individual cell to macroscopic models providing phenomenological descriptions at the population-average level. The reviewed modeling approaches include Molecular Dynamics, Particle-Based Brownian Dynamics, the Master Equation approach, Ordinary Differential Equations, and the Boolean logic abstraction. Each of these frameworks is motivated by a particular biological context and the nature of the insight being pursued. The setting of gene network dynamic models from such frameworks involves assumptions and mathematical artifacts often ignored by the non-specialist. This article aims at providing an entry point for biologists new to the field and computer scientists not acquainted with some recent biophysically-inspired models of gene regulation. The connections promoting intuition between different abstraction levels and the role that approximations play in the modeling process are highlighted throughout the paper.

Stochastic operator-splitting method for reaction-diffusion systems

Many biochemical processes at the sub-cellular level involve a small number of molecules. The local numbers of these molecules vary in space and time, and exhibit random fluctuations that can only be captured with stochastic simulations. We present a novel stochastic operator-splitting algorithm to model such reaction-diffusion phenomena. The reaction and diffusion steps employ stochastic simulation algorithms and Brownian dynamics, respectively. Through theoretical analysis, we have developed an algorithm to identify if the system is reaction-controlled, diffusion-controlled or is in an intermediate regime. The time-step size is chosen accordingly at each step of the simulation. We have used three examples to demonstrate the accuracy and robustness of the proposed algorithm. The first example deals with diffusion of two chemical species undergoing an irreversible bimolecular reaction. It is used to validate our algorithm by comparing its results with the solution obtained from a corresponding deterministic partial differential equation at low and high number of molecules. In this example, we also compare the results from our method to those obtained using a Gillespie multi-particle (GMP) method. The second example, which models simplified RNA synthesis, is used to study the performance of our algorithm in reaction- and diffusion-controlled regimes and to investigate the effects of local inhomogeneity. The third example models reaction-diffusion of CheY molecules through the cytoplasm of Escherichia coli during chemotaxis. It is used to compare the algorithm's performance against the GMP method. Our analysis demonstrates that the proposed algorithm enables accurate simulation of the kinetics of complex and spatially heterogeneous systems. It is also computationally more efficient than commonly used alternatives, such as the GMP method.

A computational approach to increase time scales in Brownian dynamics-based reaction-diffusion modeling

Particle-based Brownian dynamics simulations offer the opportunity to not only simulate diffusion of particles but also the reactions between them. They therefore provide an opportunity to integrate varied biological data into spatially explicit models of biological processes, such as signal transduction or mitosis. However, particle based reaction-diffusion methods often are hampered by the relatively small time step needed for accurate description of the reaction-diffusion framework. Such small time steps often prevent simulation times that are relevant for biological processes. It is therefore of great importance to develop reaction-diffusion methods that tolerate larger time steps while maintaining relatively high accuracy. Here, we provide an algorithm, which detects potential particle collisions prior to a BD-based particle displacement and at the same time rigorously obeys the detailed balance rule of equilibrium reactions. We can show that for reaction-diffusion processes of particles mimicking proteins, the method can increase the typical BD time step by an order of magnitude while maintaining similar accuracy in the reaction diffusion modelling.

Lattice Microbes: high-performance stochastic simulation method for the reaction-diffusion master equation

Spatial stochastic simulation is a valuable technique for studying reactions in biological systems. With the availability of high-performance computing (HPC), the method is poised to allow integration of data from structural, single-molecule and biochemical studies into coherent computational models of cells. Here, we introduce the Lattice Microbes software package for simulating such cell models on HPC systems. The software performs either well-stirred or spatially resolved stochastic simulations with approximated cytoplasmic crowding in a fast and efficient manner. Our new algorithm efficiently samples the reaction-diffusion master equation using NVIDIA graphics processing units and is shown to be two orders of magnitude faster than exact sampling for large systems while maintaining an accuracy of !0.1%. Display of cell models and animation of reaction trajectories involving millions of molecules is facilitated using a plug-in to the popular VMD visualization platform. The Lattice Microbes software is open source and available for download at http://www.scs.illinois.edu/schulten/lm

A new method for choosing the computational cell in stochastic reaction-diffusion systems

How to choose the computational compartment or cell size for the stochastic simulation of a reaction-diffusion system is still an open problem, and a number of criteria have been suggested. A generalized measure of the noise for finite-dimensional systems based on the largest eigenvalue of the covariance matrix of the number of molecules of all species has been suggested as a measure of the overall fluctuations in a multivariate system, and we apply it here to a discretized reaction-diffusion system. We show that for a broad class of first-order reaction networks this measure converges to the square root of the reciprocal of the smallest mean species number in a compartment at the steady state. We show that a suitably re-normalized measure stabilizes as the volume of a cell approaches zero, which leads to a criterion for the maximum volume of the compartments in a computational grid. We then derive a new criterion based on the sensitivity of the entire network, not just of the fastest step, that predicts a grid size that assures that the concentrations of all species converge to a spatially-uniform solution. This criterion applies for all orders of reactions and for reaction rate functions derived from singular perturbation or other reduction methods, and encompasses both diffusing and non-diffusing species. We show that this predicts the maximal allowable volume found in a linear problem, and we illustrate our results with an example motivated by anterior-posterior pattern formation in Drosophila, and with several other examples.

Accounting for diffusion in agent based models of reaction-diffusion systems with application to cytoskeletal diffusion

Diffusion plays a key role in many biochemical reaction systems seen in nature. Scenarios where diffusion behavior is critical can be seen in the cell and subcellular compartments where molecular crowding limits the interaction between particles. We investigate the application of a computational method for modeling the diffusion of molecules and macromolecules in three-dimensional solutions using agent based modeling. This method allows for realistic modeling of a system of particles with different properties such as size, diffusion coefficients, and affinity as well as the environment properties such as viscosity and geometry. Simulations using these movement probabilities yield behavior that mimics natural diffusion. Using this modeling framework, we simulate the effects of molecular crowding on effective diffusion and have validated the results of our model using Langevin dynamics simulations and note that they are in good agreement with previous experimental data. Furthermore, we investigate an extension of this framework where single discrete cells can contain multiple particles of varying size in an effort to highlight errors that can arise from discretization that lead to the unnatural behavior of particles undergoing diffusion. Subsequently, we explore various algorithms that differ in how they handle the movement of multiple particles per cell and suggest an algorithm that properly accommodates multiple particles of various sizes per cell that can replicate the natural behavior of these particles diffusing. Finally, we use the present modeling framework to investigate the effect of structural geometry on the directionality of diffusion in the cell cytoskeleton with the observation that parallel orientation in the structural geometry of actin filaments of filopodia and the branched structure of lamellipodia can give directionality to diffusion at the filopodia-lamellipodia interface.

Inertial effects in diffusion-limited reactions

Diffusion-limited reactions are commonly found in biochemical processes such as enzyme catalysis, colloid and protein aggregation and binding between different macromolecules in cells. Usually, such reactions are modeled within the Smoluchowski framework by considering purely diffusive boundary problems. However, inertial effects are not always negligible in real biological or physical media on typical observation time frames. This is all the more so for non-bulk phenomena involving physical boundaries, that introduce additional time and space constraints. In this paper, we present and test a novel numerical scheme, based on event-driven Brownian dynamics, that allows us to explore a wide range of velocity relaxation times, from the purely diffusive case to the underdamped regime. We show that our algorithm perfectly reproduces the solution of the Fokker-Planck problem with absorbing boundary conditions in all the regimes considered and is thus a good tool for studying diffusion-guided reactions in complex biological environments.

Detailed simulations of cell biology with Smoldyn 2.1

Most cellular processes depend on intracellular locations and random collisions of individual protein molecules. To model these processes, we developed algorithms to simulate the diffusion, membrane interactions, and reactions of individual molecules, and implemented these in the Smoldyn program. Compared to the popular MCell and ChemCell simulators, we found that Smoldyn was in many cases more accurate, more computationally efficient, and easier to use. Using Smoldyn, we modeled pheromone response system signaling among yeast cells of opposite mating type. This model showed that secreted Bar1 protease might help a cell identify the fittest mating partner by sharpening the pheromone concentration gradient. This model involved about 200,000 protein molecules, about 7000 cubic microns of volume, and about 75 minutes of simulated time; it took about 10 hours to run. Over the next several years, as faster computers become available, Smoldyn will allow researchers to model and explore systems the size of entire bacterial and smaller eukaryotic cells.

Stochastic off-lattice modeling of molecular self-assembly in crowded environments by Green's function reaction dynamics

The environment inside a living cell is dramatically different from that found in in vitro models, presenting a problem for computational models of biochemistry that are only beginning to capture these differences. This deviation between idealized in vitro models and more realistic intracellular conditions is particularly problematic for models of molecular self-assembly, but also specifically hard to address because the large sizes and long assembly times of biological self-assembly systems force the use of highly simplified models. We have developed a prototype of a molecular self-assembly simulator based on the Green's function reaction dynamics (GFRD) model to achieve more realistic models of assembly in the crowded conditions of the cell without unduly sacrificing tractability. We tested the model on a simple representation of dimer assembly in a two-dimensional space. Our simulations verify that the model is computationally efficient, provides a realistic quantitative model of reaction kinetics in uncrowded conditions, and exhibits expected excluded volume effects under conditions of high crowding. This work confirms the effectiveness of the GFRD technique for more realistic coarse-grained modeling of self-assembly in crowded conditions and helps lay the groundwork for exploring the effects of in vivo crowding on more complex assembly systems.

mlegp: statistical analysis for computer models of biological systems using R

Gaussian processes (GPs) are flexible statistical models commonly used for predicting output from complex computer codes. As such, GPs are well suited for the analysis of computer models of biological systems, which have been traditionally difficult to analyze due to their high-dimensional, non-linear and resource-intensive nature. We describe an R package, mlegp, that fits GPs to computer model outputs and performs sensitivity analysis to identify and characterize the effects of important model inputs.

AVAILABILITY

http://www.biomath.org/mlegp