Stochastic Reaction-diffusion Model Research Articles

Resilience for stochastic systems interacting via a quasi-degenerate network.

A stochastic reaction-diffusion model is studied on a networked support. In each patch of the network, two species are assumed to interact following a non-normal reaction scheme. When the interaction unit is replicated on a directed linear lattice, noise gets amplified via a self-consistent process, which we trace back to the degenerate spectrum of the embedding support. The same phenomenon holds when the system is bound to explore a quasidegenerate network. In this case, the eigenvalues of the Laplacian operator, which governs species diffusion, accumulate over a limited portion of the complex plane. The larger the network, the more pronounced the amplification. Beyond a critical network size, a system deemed deterministically stable, hence resilient, can develop seemingly regular patterns in the concentration amount. Non-normality and quasidegenerate networks may, therefore, amplify the inherent stochasticity and so contribute to altering the perception of resilience, as quantified via conventional deterministic methods.

Turing-like mechanism in a stochastic reaction-diffusion model recreates three dimensional vascular patterning of plant stems.

Vascular tissue in plants provides a resource distribution network for water and nutrients that exhibits remarkable diversity in patterning among different species. In many succulent plants, the vascular network includes longitudinally-oriented supplemental vascular bundles (SVBs) in the central core of the succulent stems and roots in addition to the more typical zone of vascular tissue development (vascular cambium) in a cylinder at the periphery of the succulent organ. Plant SVBs evolved in over 38 plant families often in tandem with evolutionary increases in stem and root parenchyma storage tissue, so it is of interest to understand the evolutionary-developmental processes responsible for their recurrent evolution and patterning. Previous mathematical models have successfully recreated the two-dimensional vascular patterns in stem and root cross sections, but such models have yet to recreate three-dimensional vascular patterning. Here, a stochastic reaction-diffusion model of plant vascular bundle patterning is developed in an effort to highlight a potential mechanism of three dimensional patterning–Turing pattern formation coupled with longitudinal efflux of a regulatory molecule. A relatively simple model of four or five molecules recreated empirical SVB patterns and many other common vascular arrangements. SVBs failed to develop below a threshold width of parenchymatous tissues, suggesting a mechanism of evolutionary character loss due to changes in the spatial context in which development takes place. Altered diffusion rates of the modeled activator and substrate molecules affected the number and size of the simulated SVBs. This work provides a first mathematical model employing a stochastic Turing-type mechanism that recreates three dimensional vascular patterns seen in plant stems. The model offers predictions that can be tested using molecular-genetic approaches. Evolutionary-developmental ramifications concerning evolution of diffusion rates, organ size and geometry are discussed.

Estimation of stochastic behaviour in cardiac myocytes: I. Ca2+ movements inside the cytosol and sarcoplasmic reticulum on curvilinear domains.

BackgroundSince the discovery of Ca2+ sparks and their stochastic behaviour in cardiac myocytes, models have focused on the inclusion of stochasticity in their studies. While most models pay much attention to the stochastic modelling of cytosolic Ca2+ concentration the coupling of Ca2+ sparks and blinks in a stochastic model has not been explored fully. The cell morphology in in silico studies in the past is assumed to be Cartesian, spherical or cylindrical. The application on curvilinear grids can easily address certain restrictions posed by such grid set up and provide more realistic cell morphology. In this paper, we present a stochastic reaction–diffusion model that couples Ca2+ sparks and blinks in realistic shapes of cells in curvilinear domains.Methodology: Transformation of the model was performed to the curvilinear coordinate system. The set of equations is used to produce Ca2+ waves initiated from sparks and blinks. A non-buffered and non-dyed version as well as a buffered and dyed version of these equations were studied in light of observing the dynamics on the two different systems. For comparison, results for both the Cartesian and curvilinear grids are provided.Results and conclusions: A successful demonstration of the application of curvilinear grids serving as basis for future developments.

Swelling and Diffusion during Methanol Sorption into Hydrated Nafion

Diffusion within polymer electrolyte membranes is often coincident with time-dependent processes such as swelling and polymer relaxation, which are factors that limit their ability to block molecular crossover during use. The solution-diffusion model of membrane permeation, which is the accepted theory for dense polymers, applies only to steady-state processes and does not address dynamic internal structural changes that can accompany permeation. To begin discovery of how such changes can be coupled to the permeation process, we have constructed a stochastic multiscale reaction-diffusion model that examines time-dependent methanol uptake into and swelling of hydrated Nafion. Several potential mechanisms of diffusion and polymer response are tested. The simulation predictions are compared to real-time Fourier transform infrared attenuated total reflectance spectroscopy (FTIR-ATR) absorbance reported in the literature [ Hallinan , D. T. , Jr. ; Elabd , Y. A. J. Phys. Chem. B 2007 , 111 , 13221 - 13230 ]. Of the proposed polymer response mechanisms, only one, a reaction-limited, local response to increasing methanol concentration that takes the entire experimental time frame of 600 s, produces simulated FTIR-ATR data consistent with experiment. The simulations show that water diffusion out of the membrane is minimal during methanol sorption and that changes in the measured infrared absorbances are due primarily to the increase in methanol concentration accompanied by dilution of water during swelling. Swelling involves densification of the polymer structure even as there is an overall volume expansion of the film. Potential connections between the polymer densification and molecular-level structural changes of Nafion in methanol are discussed. These results indicate that the interaction between methanol and Nafion serves to increase Nafion's capacity to accommodate large volumes of methanol-water solutions, facilitating increased permeation across the membrane relative to pure water.

An unstructured mesh convergent reaction–diffusion master equation for reversible reactions

The convergent reaction–diffusion master equation (CRDME) was recently developed to provide a lattice particle-based stochastic reaction–diffusion model that is a convergent approximation in the lattice spacing to an underlying spatially-continuous particle dynamics model. The CRDME was designed to be identical to the popular lattice reaction–diffusion master equation (RDME) model for systems with only linear reactions, while overcoming the RDME's loss of bimolecular reaction effects as the lattice spacing is taken to zero. In our original work we developed the CRDME to handle bimolecular association reactions on Cartesian grids. In this work we develop several extensions to the CRDME to facilitate the modeling of cellular processes within realistic biological domains. Foremost, we extend the CRDME to handle reversible bimolecular reactions on unstructured grids. Here we develop a generalized CRDME through discretization of the spatially continuous volume reactivity model, extending the CRDME to encompass a larger variety of particle–particle interactions. Finally, we conclude by examining several numerical examples to demonstrate the convergence and accuracy of the CRDME in approximating the volume reactivity model.

Weak universality for a class of 3d stochastic reaction–diffusion models

We establish the large scale convergence of a class of stochastic weakly nonlinear reaction-diffusion models on a three dimensional periodic domain to the dynamic Phi^3_4 model within the framework of paracontrolled distributions. Our work extends previous results of Hairer and Xu to nonlinearities with a finite amount of smoothness (in particular C^9 is enough). We use the Malliavin calculus to perform a partial chaos expansion of the stochastic terms and control their L^p norms in terms of the graphs of the standard Phi^3_4 stochastic terms.

A stochastic reaction-diffusion model for protein aggregation on DNA

Vital functions of DNA, such as transcription and packaging, depend on the proper clustering of proteins on the double strand. The present study investigates how the interplay between DNA allostery and electrostatic interactions affects protein clustering. The statistical analysis of a simple but transparent computational model reveals two major consequences of this interplay. First, depending on the protein and salt concentration, protein filaments exhibit a bimodal DNA stiffening and softening behavior. Second, within a certain domain of the control parameters, electrostatic interactions can cause energetic frustration that forces proteins to assemble in rigid spiral configurations. Such spiral filaments might trigger both positive and negative supercoiling, which can ultimately promote gene compaction and regulate the promoter. It has been experimentally shown that bacterial histone-like proteins assemble in similar spiral patterns and/or exhibit the same bimodal behavior. The proposed model can, thus, provide computational insights into the physical mechanisms used by proteins to control the mechanical properties of the DNA.

A novel stochastic reaction-diffusion model of Ca2+ blink in cardiac myocytes

A novel stochastic reaction-diffusion model of Ca2+ blink in cardiac myocytes

Catalytic conversion in nanoporous materials: Concentration oscillations and spatial correlations due to inhibited transport and intermolecular interactions.

We show that steady-state catalytic conversion in nanoporous materials can occur in a quasi-counter-diffusion mode with the reactant (product) concentration strongly decaying (growing) into the pore, but also with oscillations in the total concentration. These oscillations reflect the response of the fluid to the transition from an extended to a confined environment near the pore opening. We focus on the regime of strongly inhibited transport in narrow pores corresponding to single-file diffusion. Here, limited penetration of the reactant into the pores and the associated low reaction yield is impacted by strong spatial correlations induced by both reaction (non-equilibrium correlations) and also by intermolecular interactions (thermodynamic correlations). We develop a generalized hydrodynamic formulation to effectively describe inhibited transport accounting for the effect of these correlations, and incorporate this description of transport into appropriate reaction-diffusion equations. These equations accurately describe both shorter-range concentration oscillations near the pore opening and the longer-range mesoscale variation of concentration profiles in the pore (and thus also describe reaction yield). Success of the analytic theory is validated by comparison with a precise kinetic Monte Carlo simulation of an appropriate molecular-level stochastic reaction-diffusion model. This work elucidates unconventional chemical kinetics in interacting confined systems.

The asymptotic behavior of a stochastic vaccination model with backward bifurcation

In this paper, we aim to explore the effect of environmental fluctuation on a deterministic vaccination model that exhibits backward bifurcation. First, we prove that the stochastic model admits a unique and global positive solution. Next, by investigating the asymptotic behavior of the stochastic model around the disease free equilibrium, we find that due to the influence of environmental fluctuation the disease dies out even though it may still persist for the deterministic model with backward bifurcation. We also estimate the probability distribution for the disease extinction time. Then, we derive sufficient conditions for the solution of the stochastic model to fluctuate around the endemic equilibrium, and thus we prove the ergodicity of the stochastic model. Our theoretical results are verified by computer simulations and numerical comparison results of the stochastic model and deterministic version is also given. Finally, numerical simulations of a stochastic reaction diffusion model with multiplicative noise are presented to illustrate the combined effects of spatial movement of individuals and environmental noise on the spread of disease.

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.

Hybrid stochastic simulation of reaction-diffusion systems with slow and fast dynamics.

In this paper, we present a novel hybrid method to simulate discrete stochastic reaction-diffusion models arising in biochemical signaling pathways. We study moderately stiff systems, for which we can partition each reaction or diffusion channel into either a slow or fast subset, based on its propensity. Numerical approaches missing this distinction are often limited with respect to computational run time or approximation quality. We design an approximate scheme that remedies these pitfalls by using a new blending strategy of the well-established inhomogeneous stochastic simulation algorithm and the tau-leaping simulation method. The advantages of our hybrid simulation algorithm are demonstrated on three benchmarking systems, with special focus on approximation accuracy and efficiency.

Small-scale properties of a stochastic cubic-autocatalytic reaction-diffusion model.

We investigate the small-scale properties of a stochastic cubic-autocatalytic reaction-diffusion (CARD) model using renormalization techniques. We renormalize noise-induced ultraviolet divergences and obtain β functions for the decay rate and coupling at one loop. Assuming colored (power-law) noise, our results show that the behavior of both decay rate and coupling with scale depends crucially on the noise exponent. Interpreting the CARD model as a proxy for a (very simple) living system, our results suggest that power-law correlations in environmental fluctuations can both decrease or increase the growth of structures at smaller scales.

Hybrid approaches for multiple-species stochastic reaction–diffusion models

Reaction–diffusion models are used to describe systems in fields as diverse as physics, chemistry, ecology and biology. The fundamental quantities in such models are individual entities such as atoms and molecules, bacteria, cells or animals, which move and/or react in a stochastic manner. If the number of entities is large, accounting for each individual is inefficient, and often partial differential equation (PDE) models are used in which the stochastic behaviour of individuals is replaced by a description of the averaged, or mean behaviour of the system. In some situations the number of individuals is large in certain regions and small in others. In such cases, a stochastic model may be inefficient in one region, and a PDE model inaccurate in another. To overcome this problem, we develop a scheme which couples a stochastic reaction–diffusion system in one part of the domain with its mean field analogue, i.e. a discretised PDE model, in the other part of the domain. The interface in between the two domains occupies exactly one lattice site and is chosen such that the mean field description is still accurate there. In this way errors due to the flux between the domains are small. Our scheme can account for multiple dynamic interfaces separating multiple stochastic and deterministic domains, and the coupling between the domains conserves the total number of particles. The method preserves stochastic features such as extinction not observable in the mean field description, and is significantly faster to simulate on a computer than the pure stochastic model.

Modeling the growth of an altered layer in mineral weathering

A stochastic reaction-diffusion model on a lattice is introduced to describe the growth kinetics of an altered layer in the weathering of a mineral. Particles R represent H2O that permanently fills the outer surface and diffuse on M (mineral) and A (altered) sites with coefficients DM and DA, respectively. The transformation M+R→A occurs with rate r, representing the irreversible formation of the altered material in a region of molecular size, viz. the lattice site of size a. These assumptions agree with predictions of the interfacial dissolution–reprecipitation mechanism, although the model does not describe the chemistry of dissolution reactions or precipitation processes. Scaling concepts are used to distinguish kinetic regimes and their crossovers, and are supported by simulation results. In the short time reactive regime, the thickness of the altered layer increases linearly in time and filling of that layer by particles R is high. In the long time diffusive regime, the altered layer thickness grows as DAt1/2. Modeling of single crystals require very small values of DM, which produces atomically narrow interfaces between the altered material and the mineral and absence of R in the latter, in agreement with recent experimental results. If r<DM/a2, an interface of thickness lAM∼DM/r1/2 is found, with subsurface fluid penetration and increased growth velocity of the altered layer by a factor lAM/a, but no fluid in the bulk mineral. Estimates of the order of magnitude of transformation rates and of diffusion coefficients are obtained by application of the model to some recently studied systems: calcite dissolution, labradorite weathering, and silicate glass weathering. Effects of dissolution of the altered layer are analyzed. Significant differences between the model and leached layer theories are discussed.

Stochastic Turing patterns: analysis of compartment-based approaches.

Turing patterns can be observed in reaction-diffusion systems where chemical species have different diffusion constants. In recent years, several studies investigated the effects of noise on Turing patterns and showed that the parameter regimes, for which stochastic Turing patterns are observed, can be larger than the parameter regimes predicted by deterministic models, which are written in terms of partial differential equations (PDEs) for species concentrations. A common stochastic reaction-diffusion approach is written in terms of compartment-based (lattice-based) models, where the domain of interest is divided into artificial compartments and the number of molecules in each compartment is simulated. In this paper, the dependence of stochastic Turing patterns on the compartment size is investigated. It has previously been shown (for relatively simpler systems) that a modeler should not choose compartment sizes which are too small or too large, and that the optimal compartment size depends on the diffusion constant. Taking these results into account, we propose and study a compartment-based model of Turing patterns where each chemical species is described using a different set of compartments. It is shown that the parameter regions where spatial patterns form are different from the regions obtained by classical deterministic PDE-based models, but they are also different from the results obtained for the stochastic reaction-diffusion models which use a single set of compartments for all chemical species. In particular, it is argued that some previously reported results on the effect of noise on Turing patterns in biological systems need to be reinterpreted.

Molecular dependence of hippocampal bidirectional plasticity

Hippocampal plasticity inducing protocols vary in their activation of signal transduction pathways and resulting dependence on signaling molecules. Many long-term potentiation (LTP) inducing paradigms require both memory kinases: calcium-calmodulin dependent protein kinase II (CaMKII) and protein kinase A (PKA) signaling pathways. Both CaMKII and PKA are activated by calcium bound calmodulin, CaMKII directly, and PKA indirectly via adenylyl cyclase production of cAMP. The adenylyl cyclase is also activated by stimulation of beta-adrenergic receptor (βAR). Curiously, not all forms of cAMP dependent LTP require PKA. Specifically, LTP induced with isoproterenol (ISO) bath application followed by 100 Hz stimulation is PKA independent [1,2], whereas both ISO followed by 3 min of 5 Hz [2] and four trains of 100 Hz stimulation do require PKA [3]. Furthermore, recent experiments suggest a novel βAR signaling pathway, in which PKA phosphorylation prevents activation of G protein subtype s (Gs), and switches the βAR to the Gi subtype of G [4], leading to recruitment of another memory kinase: extracellular signal-regulated kinase (ERK). To investigate how temporal pattern of synaptic activation determines which signaling pathways are activated, we employed a multi-compartmental stochastic reaction-diffusion model of calcium and cAMP activated signaling pathways. We investigated two model variants, one which took into account only βARcoupling to Gs, and the other which included the novel switching pathway. We performed simulations of four trains of 100 Hz, and isoproterenol application followed by either 1 train of 100 Hz or 3 min of 5 Hz, all of which produce long lasting LTP. We also simulated 3 stimulation paradigms which do not produce LTP, including 1 Hz stimulation which produces LTD. Our results using the first model could not explain and predict the direction of synaptic plasticity based only on CaMKII, PKAc and phosphatases activity, especially for paradigms that do not induce plasticity, such as bath application of isoproteronol. Experimentally, isoproteronol does not cause LTP despite βAR activation resulting in cAMP synthesis [2], suggesting that βAR undergo PKA dependent desensitization, that changes their coupling from Gs (stimulating adenylyl cyclases and cAMP synthesis) to Gi. Gi coupling both inhibits adenylyl cyclase, limiting the elevation in cAMP in response to bath application, and provides free Giβγ, which recruits molecules leading to ERK [5]. Including PKA dependent switching from Gs to Gi resulted in a model which could explain which molecules are activated in 5 LTP induction paradigms, 1 LTD induction paradigms, and 2 induction paradigms that do not produce plasticity. The model also allowed for proposing a measure of plasticity direction. The measure took into account activation of CaMKII, PKA, Giβγ (representing ERK), Epac and influence of phosphatase activity.

Editorial: Special Issue on Stochastic Modelling of Reaction–Diffusion Processes in Biology

Editorial: Special Issue on Stochastic Modelling of Reaction–Diffusion Processes in Biology

Analysis of the Two-Regime Method on Square Meshes

The two-regime method (TRM) has been recently developed for optimizing stochastic reaction-diffusion simulations [M. Flegg, J. Chapman, and R. Erban, J. Roy. Soc. Interface, 9 (2012), pp. 859--868]. It is a multiscale (hybrid) algorithm which uses stochastic reaction-diffusion models with different levels of detail in different parts of the computational domain. The coupling condition on the interface between different modeling regimes of the TRM was previously derived for one-dimensional models. In this paper, the TRM is generalized to higher dimensional reaction-diffusion systems. Coupling Brownian dynamics models with compartment-based models on regular (square) two-dimensional lattices is studied in detail. In this case, the interface between different modeling regimes contains either flat parts or right-angle corners. Both cases are studied in the paper. For flat interfaces, it is shown that the one-dimensional theory can be used along the line perpendicular to the TRM interface. In the direction tan...

A First-Passage Kinetic Monte Carlo method for reaction–drift–diffusion processes

Stochastic reaction–diffusion models are now a popular tool for studying physical systems in which both the explicit diffusion of molecules and noise in the chemical reaction process play important roles. The Smoluchowski diffusion-limited reaction model (SDLR) is one of several that have been used to study biological systems. Exact realizations of the underlying stochastic processes described by the SDLR model can be generated by the recently proposed First-Passage Kinetic Monte Carlo (FPKMC) method. This exactness relies on sampling analytical solutions to one and two-body diffusion equations in simplified protective domains.In this work we extend the FPKMC to allow for drift arising from fixed, background potentials. As the corresponding Fokker–Planck equations that describe the motion of each molecule can no longer be solved analytically, we develop a hybrid method that discretizes the protective domains. The discretization is chosen so that the drift–diffusion of each molecule within its protective domain is approximated by a continuous-time random walk on a lattice. New lattices are defined dynamically as the protective domains are updated, hence we will refer to our method as Dynamic Lattice FPKMC or DL-FPKMC. We focus primarily on the one-dimensional case in this manuscript, and demonstrate the numerical convergence and accuracy of our method in this case for both smooth and discontinuous potentials. We also present applications of our method, which illustrate the impact of drift on reaction kinetics.