Search

Abstract Background Experiments in silico using stochastic reaction-diffusion models have emerged as an important tool in molecular systems biology. Designing computational software for such applications poses several challenges. Firstly, realistic lattice-based modeling for biological applications requires a consistent way of handling complex geometries, including curved inner- and outer boundaries. Secondly, spatiotemporal stochastic simulations are computationally expensive due to the fast time scales of individual reaction- and diffusion events when compared to the biological phenomena of actual interest. We therefore argue that simulation software needs to be both computationally efficient, employing sophisticated algorithms, yet in the same time flexible in order to meet present and future needs of increasingly complex biological modeling. Results We have developed URDME, a flexible software framework for general stochastic reaction-transport modeling and simulation. URDME uses Unstructured triangular and tetrahedral meshes to resolve general geometries, and relies on the Reaction-Diffusion Master Equation formalism to model the processes under study. An interface to a mature geometry and mesh handling external software (Comsol Multiphysics) provides for a stable and interactive environment for model construction. The core simulation routines are logically separated from the model building interface and written in a low-level language for computational efficiency. The connection to the geometry handling software is realized via a Matlab interface which facilitates script computing, data management, and post-processing. For practitioners, the software therefore behaves much as an interactive Matlab toolbox. At the same time, it is possible to modify and extend URDME with newly developed simulation routines. Since the overall design effectively hides the complexity of managing the geometry and meshes, this means that newly developed methods may be tested in a realistic setting already at an early stage of development. Conclusions In this paper we demonstrate, in a series of examples with high relevance to the molecular systems biology community, that the proposed software framework is a useful tool for both practitioners and developers of spatial stochastic simulation algorithms. Through the combined efforts of algorithm development and improved modeling accuracy, increasingly complex biological models become feasible to study through computational methods. URDME is freely available at http://www.urdme.org.

We present a new approach for efficient process synchronization in parallel discrete event simulation on multicore computers. We aim specifically at simulation of spatially extended stochastic system models where time intervals between successive inter-process events are highly variable and without lower bounds: This includes models governed by the mesoscopic Reaction-Diffusion Master Equation (RDME). A central part of our approach is a mechanism for optimism control, in which each process disseminates accurate information about timestamps of its future outgoing interprocess events to its neighbours. This information gives each process a precise basis for deciding when to pause local processing to reduce the risk of expensive rollbacks caused by future "delayed" incoming events. We apply our approach to a natural parallelization of the Next Subvolume Method (NSM) for simulating systems obeying RDME. Since this natural parallelization does not expose accurate timestamps of future interprocess events, we restructure it to expose such information, resulting in a simulation algorithm called Refined Parallel NSM (Refined PNSM). We have implemented Refined PNSM in a parallel simulator for spatial extended Markovian processes. On 32 cores, it achieves an efficiency ranging between 43—95% for large models, and on average 37% for small models, compared to an efficient sequential simulation without any code for parallelization. It is shown that the gain of restructuring the naive parallelization into Refined PNSM more than outweighs its overhead. We also show that our resulting simulator is superior in performance to existing simulators on multicores for comparable models. [ABSTRACT FROM AUTHOR]

Abstract We consider a multiscale strategy addressing the disparate scales in the Landau–Lifschitz equations in micromagnetism. At the microscopic scale, the dynamics of magnetic moments are driven by a high frequency field. On the macroscopic scale we are interested in simulating the dynamics of the magnetisation without fully resolving the microscopic scales. The method follows the framework of heterogeneous multiscale methods and it has two main ingredients: a micro- and a macroscale model. The microscopic model is assumed to be known exactly whereas the macromodel is incomplete as it lacks effective quantities. The two models use different temporal and spatial scales and effective parameter values for the macromodel are computed on the fly, allowing for improved efficiency over traditional one-scale schemes. For the analysis, we consider a single spin under a high frequency field and show that effective quantities can be obtained accurately with step-sizes much larger than the size of the microscopic scales required to resolve the microscopic features. Numerical results both for a single magnetic particle as well as a chain of interacting magnetic particles are given to validate the theory. [ABSTRACT FROM AUTHOR]

Neuronal models based on the Hodgkin-Huxley equation form a fundamental framework in the field of computational neuroscience. While the neuronal state is often modelled deterministically, experimental recordings show stochastic fluctuations, presumably driven by molecular noise from the underlying microphysical conditions. In turn, the firing of individual neurons gives rise to an electric field in extracellular space, also thought to affect the firing pattern of nearby neurons. We develop a multiscale model which combines a stochastic ion channel gating process taking place on the neuronal membrane, together with the propagation of an action potential along the neuronal structure. We also devise a numerical method relying on a split-step strategy which effectively couples these two processes and we experimentally test the feasibility of this approach. We finally also explain how the approach can be extended with Maxwell's equations to allow the potential to be propagated in extracellular space. [ABSTRACT FROM AUTHOR]

At sites of hypoxia, neutrophils migrate and accumulate at angiogenic vasculature, where they closely interact with endothelial sprouts and macrophages. Angiogenesis, the growth of new blood vessels, is a complex process requiring the orchestration of numerous different cell types, growth factors, and chemokines. Some of the recently acknowledged actors in this process are immune cells. They accumulate at hypoxic sites, but the kinetics, dynamics, and regulation of that trafficking are unknown. In this study, we used intravital and live cell imaging to understand how neutrophils and macrophages migrate and behave at angiogenic sites. We developed two reproducible models of angiogenesis: one by transplanting isolated and hypoxic pancreatic islets into the cremaster muscles of mice, and another by in vitro coculturing of mouse aortic rings with neutrophils. In vivo imaging of the hypoxic site revealed recruitment of neutrophils and macrophages, which occurred in parallel, with depletion of one subset not affecting the accumulation of the other. We found, by cell tracking and statistical analyses, that neutrophils migrated in a directional manner to "angiogenic hotspots" around the islet where endothelial sprouting occurs, which was confirmed in the in vitro model of angiogenesis and is dependent on CXCL12 signaling. Intimate interactions between neutrophils and neovessels were prevalent, and neutrophil depletion greatly hampered vessel growth. Macrophages were less motile and attained supportive positions around the neovessels. Here, we present two novel in vivo and in vitro imaging models to study leukocyte behavior and actions during angiogenesis. These models unveiled that neutrophil migration at a hypoxic site was guided by signals emanating from sprouting endothelium where these immune cells gathered at "angiogenic hotspots" at which vascular growth occurred. [ABSTRACT FROM AUTHOR]

We present a computational modeling framework for data-driven simulations and analysis of infectious disease spread in large populations. For the purpose of efficient simulations, we devise a parallel solution algorithm targeting multi-socket shared-memory architectures. The model integrates infectious dynamics as continuous-time Markov chains and available data such as animal movements or aging are incorporated as externally defined events. To bring out parallelism and accelerate the computations, we decompose the spatial domain and optimize cross-boundary communication using dependency-aware task scheduling. Using registered livestock data at a high spatiotemporal resolution, we demonstrate that our approach not only is resilient to varying model configurations but also scales on all physical cores at realistic workloads. Finally, we show that these very features enable the solution of inverse problems on national scales. [ABSTRACT FROM AUTHOR]

Subdifision has been proposed as an explanation of various kinetic phenomena inside living cells. In order to fascilitate large-scale computational studies of subdiffusive chemical processes, we extend a recently suggested mesoscopic model of subdiffusion into an accurate and consistent reaction-subdiffusion computational framework. Two different possible models of chemical reaction are revealed and some basic dynamic properties are derived. In certain cases those mesoscopic models have a direct interpretation at the macroscopic level as fractional partial differential equations in a bounded time interval. Through analysis and numerical experiments we estimate the macroscopic effects of reactions under subdiffusive mixing. The models display properties also observed in experiments: for a short time interval the behavior of the diffusion and the reaction is ordinary, in an intermediate interval the behavior is anomalous, and at long times the behavior is ordinary again. [ABSTRACT FROM AUTHOR]

European Union legislation requires member states to keep national databases of all bovine animals. This allows for disease spread models that includes the time-varying contact network and population demographic. However, performing data-driven simulations with a high degree of detail are computationally challenging. We have developed an efficient and flexible discrete-event simulator SimInf for stochastic disease spread modelling that divides work among multiple processors to accelerate the computations. The model integrates disease dynamics as continuous-time Markov chains and livestock data as events. In this study, all Swedish livestock data (births, movements and slaughter) from July 1st 2005 to December 31st 2013 were included in the simulations. Verotoxigenic Escherichia coli O157:H7 (VTEC O157) are capable of causing serious illness in humans. Cattle are considered to be the main reservoir of the bacteria. A better understanding of the epidemiology in the cattle population is necessary to be able to design and deploy targeted measures to reduce the VTEC O157 prevalence and, subsequently, human exposure. To explore the spread of VTEC O157 in the entire Swedish cattle population during the period under study, a within- and betweenherd disease spread model was used. Real livestock data was incorporated to model demographics of the population. Cattle were moved between herds according to real movement data. The results showed that the spatial pattern in prevalence may be due to regional differences in livestock movements. However, the movements, births and slaughter of cattle could not explain the temporal pattern of VTEC O157 prevalence in cattle, despite their inherently distinct seasonality. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

The classical method of determining the atomic structure of complex molecules by analyzing diffraction patterns is currently undergoing drastic developments. Modern techniques for producing extremely bright and coherent X-ray lasers allow a beam of streaming particles to be intercepted and hit by an ultrashort high-energy X-ray beam. Through machine learning methods the data thus collected can be transformed into a three-dimensional volumetric intensity map of the particle itself. The computational complexity associated with this problem is very high such that clusters of data parallel accelerators are required. We have implemented a distributed and highly efficient algorithm for the inversion of large collections of diffraction patterns targeting clusters of hundreds of GPUs. With the expected enormous amount of diffraction data to be produced in the foreseeable future, this is the required scale to approach real-time processing of data at the beam site. Using both real and synthetic data we look at the scaling properties of the application and discuss the overall computational viability of this exciting and novel imaging technique. [ABSTRACT FROM PUBLISHER]

Background: Experiments in silico using stochastic reaction-diffusion models have emerged as an important tool in molecular systems biology. Designing computational software for such applications poses several challenges. Firstly, realistic lattice-based modeling for biological applications requires a consistent way of handling complex geometries, including curved inner- and outer boundaries. Secondly, spatiotemporal stochastic simulations are computationally expensive due to the fast time scales of individual reaction- and diffusion events when compared to the biological phenomena of actual interest. We therefore argue that simulation software needs to be both computationally efficient, employing sophisticated algorithms, yet in the same time flexible in order to meet present and future needs of increasingly complex biological modeling. Results: We have developed URDME, a flexible software framework for general stochastic reaction-transport modeling and simulation. URDME uses Unstructured triangular and tetrahedral meshes to resolve general geometries, and relies on the Reaction-Diffusion Master Equation formalism to model the processes under study. An interface to a mature geometry and mesh handling external software (Comsol Multiphysics) provides for a stable and interactive environment for model construction. The core simulation routines are logically separated from the model building interface and written in a low-level language for computational efficiency. The connection to the geometry handling software is realized via a Matlab interface which facilitates script computing, data management, and post-processing. For practitioners, the software therefore behaves much as an interactive Matlab toolbox. At the same time, it is possible to modify and extend URDME with newly developed simulation routines. Since the overall design effectively hides the complexity of managing the geometry and meshes, this means that newly developed methods may be tested in a realistic setting already at an early stage of development. Conclusions: In this paper we demonstrate, in a series of examples with high relevance to the molecular systems biology community, that the proposed software framework is a useful tool for both practitioners and developers of spatial stochastic simulation algorithms. Through the combined efforts of algorithm development and improved modeling accuracy, increasingly complex biological models become feasible to study through computational methods. URDME is freely available at http://www.urdme.org. [ABSTRACT FROM AUTHOR]

A version of the time-parallel algorithm parareal is analyzed and applied to stochastic models in chemical kinetics. A fast predictor at the macroscopic scale (evaluated in serial) is available in the form of the usual reaction rate equation. A stochastic simulation algorithm is used to obtain an exact realization of the process at the mesoscopic scale (in parallel). The underlying stochastic description is a jump process driven by the Poisson measure. A convergence result in this arguably difficult setting is established, suggesting that a homogenization of the solution is advantageous. We devise a simple but highly general such technique. Three numerical experiments on models representative to the field of computational systems biology illustrate the method. For nonstiff problems, it is shown that the method is able to quickly converge even when stochastic effects are present. For stiff problems, we are instead able to obtain fast convergence to a homogenized solution. Overall, the method builds an attractive bridge between, on the one hand, macroscopic deterministic scales and, on the other hand, mesoscopic stochastic ones. This construction is clearly possible to apply also to stochastic models within other fields. [ABSTRACT FROM AUTHOR]

Abstract: The master equation of chemical reactions is an accurate stochastic description of general systems in chemistry. For reacting species this is a differential-difference equation in dimensions, exactly soluble for very simple systems only. We propose and analyze a novel solution strategy in the form of a Galerkin spectral method with a favorable choice of basis functions. A spectral approximation theory in the corresponding spaces is developed and the issue of stability is discussed. The convergence properties of the method are demonstrated by the numerical solution of two model problems with known solutions and a third problem for which no solution is known. It is shown that the method is effective and accurate, providing a viable alternative to other solution methods when the dimensionality is not too high. [Copyright &y& Elsevier]

We model stochastic chemical systems with diffusion by a reaction-diffusion master equation. On a macroscopic level, the governing equation is a reaction-diffusion equation for the averages of the chemical species. On a mesoscopic level, the master equation for a well stirred chemical system is combined with a discretized Brownian motion in space to obtain the reaction-diffusion master equation. The space is covered in our method by an unstructured mesh, and the diffusion coefficients on the mesoscale are obtained from a finite element discretization of the Laplace operator on the macroscale. The resulting method is a flexible hybrid algorithm in that the diffusion can be handled either on the meso- or on the macroscale level. The accuracy and the efficiency of the method are illustrated in three numerical examples inspired by molecular biology. [ABSTRACT FROM AUTHOR]

We develop a mesoscopic modeling framework for diffusion in a crowded environment, particularly targeting applications in the modeling of living cells. Through homogenization techniques we effectively coarse grain a detailed microscopic description into a previously developed internal state diffusive framework. The observables in the mesoscopic model correspond to solutions of macroscopic partial differential equations driven by stochastically varying diffusion fields in space and time. Analytical solutions and numerical experiments illustrate the framework. [ABSTRACT FROM AUTHOR]

Modern technology for producing extremely bright and coherent x-ray laser pulses provides the possibility to acquire a large number of diffraction patterns from individual biological nanoparticles, including proteins, viruses, and DNA. These two-dimensional diffraction patterns can be practically reconstructed and retrieved down to a resolution of a few angstroms. In principle, a sufficiently large collection of diffraction patterns will contain the required information for a full three-dimensional reconstruction of the biomolecule. The computational methodology for this reconstruction task is still under development and highly resolved reconstructions have not yet been produced. We analyze the expansion-maximization-compression scheme, the current state of the art approach for this very challenging application, by isolating different sources of resolution-limiting factors. Through numerical experiments on synthetic data we evaluate their respective impact. We reach conclusions of relevance for handling actual experimental data, and we also point out certain improvements to the underlying estimation algorithm. We also introduce a practically applicable computational methodology in the form of bootstrap procedures for assessing reconstruction uncertainty in the real data case. We evaluate the sharpness of this approach and argue that this type of procedure will be critical in the near future when handling the increasing amount of data. [ABSTRACT FROM AUTHOR]

This paper offers a qualitative insight into the convergence of Bayesian parameter inference in a setup which mimics the modeling of the spread of a disease with associated disease measurements. Specifically, we are interested in the Bayesian model's convergence with increasing amounts of data under measurement limitations. Depending on how weakly informative the disease measurements are, we offer a kind of 'best case' as well as a 'worst case' analysis where, in the former case, we assume that the prevalence is directly accessible, while in the latter that only a binary signal corresponding to a prevalence detection threshold is available. Both cases are studied under an assumed so-called linear noise approximation as to the true dynamics. Numerical experiments test the sharpness of our results when confronted with more realistic situations for which analytical results are unavailable.

Modern Flash X-ray diffraction Imaging (FXI) acquires diffraction signals from single biomolecules at a high repetition rate from X-ray Free Electron Lasers (XFELs), easily obtaining millions of 2D diffraction patterns from a single experiment. Due to the stochastic nature of FXI experiments and the massive volumes of data, retrieving 3D electron densities from raw 2D diffraction patterns is a challenging and time-consuming task. We propose a semi-automatic data analysis pipeline for FXI experiments, which includes four steps: hit-finding and preliminary filtering, pattern classification, 3D Fourier reconstruction, and post-analysis. We also include a recently developed bootstrap methodology in the post-analysis step for uncertainty analysis and quality control. To achieve the best possible resolution, we further suggest using background subtraction, signal windowing, and convex optimization techniques when retrieving the Fourier phases in the post-analysis step. As an application example, we quantified the 3D electron structure of the PR772 virus using the proposed data analysis pipeline. The retrieved structure was above the detector edge resolution and clearly showed the pseudo-icosahedral capsid of the PR772.

Current Flash X-ray single-particle diffraction Imaging (FXI) experiments, which operate on modern X-ray Free Electron Lasers (XFELs), can record millions of interpretable diffraction patterns from individual biomolecules per day. Due to the practical limitations with the FXI technology, those patterns will to a varying degree include scatterings from contaminated samples. Also, the heterogeneity of the sample biomolecules is unavoidable and complicates data processing. Reducing the data volumes and selecting high-quality single-molecule patterns are therefore critical steps in the experimental setup. In this paper, we present two supervised template-based learning methods for classifying FXI patterns. Our Eigen-Image and Log-Likelihood classifier can find the best-matched template for a single-molecule pattern within a few milliseconds. It is also straightforward to parallelize them so as to match the XFEL repetition rate fully, thereby enabling processing at site. The methods perform in a stable way on various kinds of synthetic data. As a practical example we tested our methods on a real mimivirus dataset, obtaining a convincing classification accuracy of 0.9.

The processes taking place inside the living cell are now understood to the point where predictive computational models can be used to gain detailed understanding of important biological phenomena. A key challenge is to extrapolate this detailed knowledge of the individual cell to be able to explain at the population level how cells interact and respond with each other and their environment. In particular, the goal is to understand how organisms develop, maintain and repair functional tissues and organs. In this paper, we propose a novel computational framework for modelling populations of interacting cells. Our framework incorporates mechanistic, constitutive descriptions of biomechanical properties of the cell population, and uses a coarse-graining approach to derive individual rate laws that enable propagation of the population through time. Thanks to its multiscale nature, the resulting simulation algorithm is extremely scalable and highly efficient. As highlighted in our computational examples, the framework is also very flexible and may straightforwardly be coupled with continuous-time descriptions of biochemical signalling within, and between, individual cells., Competing Interests: The authors declare no competing interests.

Nature presents multiple intriguing examples of processes that proceed with high precision and regularity. This remarkable stability is frequently counter to modellers' experience with the inherent stochasticity of chemical reactions in the regime of low-copy numbers. Moreover, the effects of noise and nonlinearities can lead to 'counterintuitive' behaviour, as demonstrated for a basic enzymatic reaction scheme that can display stochastic focusing (SF). Under the assumption of rapid signal fluctuations, SF has been shown to convert a graded response into a threshold mechanism, thus attenuating the detrimental effects of signal noise. However, when the rapid fluctuation assumption is violated, this gain in sensitivity is generally obtained at the cost of very large product variance, and this unpredictable behaviour may be one possible explanation of why, more than a decade after its introduction, SF has still not been observed in real biochemical systems. In this work, we explore the noise properties of a simple enzymatic reaction mechanism with a small and fluctuating number of active enzymes that behaves as a high-gain, noisy amplifier due to SF caused by slow enzyme fluctuations. We then show that the inclusion of a plausible negative feedback mechanism turns the system from a noisy signal detector to a strong homeostatic mechanism by exchanging high gain with strong attenuation in output noise and robustness to parameter variations. Moreover, we observe that the discrepancy between deterministic and stochastic descriptions of stochastically focused systems in the evolution of the means almost completely disappears, despite very low molecule counts and the additional nonlinearity due to feedback. The reaction mechanism considered here can provide a possible resolution to the apparent conflict between intrinsic noise and high precision in critical intracellular processes., (© 2015 The Author(s).)