An adaptive tau-leaping method for stochastic simulations of reaction-diffusion systems

Computational and Systems Biology are recently experiencing a rapid development. Mathematical modeling is a key tool in analyzing critical biological processes in cells and living organisms. In particular, stochastic models are essential to accurately describe the biochemical processes in the cell. However, stochastic models are computationally much more challenging than the traditional deterministic models. Also, the sheer scale of biological processes makes efficiency of the simulation a key issue. In this thesis we study the numerical solution of a continuous stochastic model of biochemical systems, the Chemical Langevin Equation. We propose an adaptive step-size method for the Euler-Maruyama scheme applied to small noise problems. The adaptive technique is p-mean convergent and computes simultaneously all the trajectories, using the same time-step on all trajectories. Our adaptive algorithm is tested on several examples arising in applications and shown to perform much better than the fixed step-size schemes.

Stochastic modeling and simulation of biochemical systems are topics of high interest in Computational Biology. Stochastic mathematical models are critical in accurately capturing the variability observed experimentally in cellular processes, in particular when some species have low molecular numbers. Many, realistic biochemical networks exhibit stiffness, due to the presence of multiple time-scales. For such networks explicit simulation methods are computationally quite intensive. In this thesis, we introduce an improved implicit tau-leaping strategy for the simulation of stochastic biochemical kinetic models. Numerical tests on various biochemical systems of interest in applications show the efficiency of our method.

Stochastic modeling and simulation are essential tools for studying cellular processes. The dynamics of spatially heterogeneous biochemical systems with species in low amounts is governed by a discrete, stochastic model, the Reaction-Diffusion Master Equation (RDME). The Inhomogeneous Stochastic Simulation Algorithm (ISSA) is an exact numerical method for the RDME, but is prohibitively slow as it simulates every chemical reaction and diffusion event. To overcome this difficulty, an approximate strategy, the tau-leaping scheme, was developed that steps over multiple reactions and diffusion events. Mathematical models of biochemical systems are often prone to stiffness, thus computationally challenging. In this thesis, we propose an adaptive time-stepping scheme for the tau-leaping method for the RDME. This strategy is compared to the ISSA, for several models of interest. The numerical results show that the proposed adaptive technique significantly speeds-up the simulation, while maintaining excellent accuracy of the numerical solution.

Fine-Grained Parallel and Distributed Spatial Stochastic Simulation of Biological Reactions

This study examines the performance of integration methods for hybrid simulation of large and complex structural systems in the context of structural collapse due to seismic excitations. The target application is not necessarily for real-time testing, but rather for models that involve large-scale physical sub-structures and highly nonlinear numerical models. Four case studies are presented and discussed. In the first case study, the accuracy of integration schemes including two widely used methods, namely, modified version of the implicit Newmark with fixed-number of iteration (iterative) and the operator-splitting (non-iterative) is examined through pure numerical simulations. The second case study presents the results of 10 hybrid simulations repeated with the two aforementioned integration methods considering various time steps and fixed-number of iterations for the iterative integration method. The physical sub-structure in these tests consists of a single-degree-of-freedom (SDOF) cantilever column with replaceable steel coupons that provides repeatable highlynonlinear behavior including fracture-type strength and stiffness degradations. In case study three, the implicit Newmark with fixed-number of iterations is applied for hybrid simulations of a 1:2 scale steel moment frame that includes a relatively complex nonlinear numerical substructure. Lastly, a more complex numerical substructure is considered by constructing a nonlinear computational model of a moment frame coupled to a hybrid model of a 1:2 scale steel gravity frame. The last two case studies are conducted on the same porotype structure and the selection of time steps and fixed number of iterations are closely examined in pre-test simulations. The generated unbalance forces is used as an index to track the equilibrium error and predict the accuracy and stability of the simulations.

The random batch Ewald (RBE) is an efficient and accurate method for molecular dynamics (MD) simulations of physical systems at the nano/microscale. The method shows great potential to solve the computational bottleneck of long-range interactions, motivating a necessity to accelerate short-range components of the nonbonded interactions for a further speedup of MD simulations. In this work, we present an improved RBE method for the nonbonding interactions by introducing the random batch idea to constructing neighbor lists for the treatment of both the short-range part of the Ewald splitting and the Lennard-Jones potential. The efficiency that the novel neighbor list algorithm owes to the stochastic minibatch strategy can significantly reduce the total number of neighbors. We obtain the error estimate and convergence by theoretical analysis and implement the improved RBE method in the LAMMPS package. Benchmark simulations are performed to demonstrate the accuracy and stability of the algorithm. Numerical tests on computer performance by conducting large-scaled MD simulations for systems including up to 0.1 billion water molecules are run on massive clusters with up to 50000 CPU cores, demonstrating the attractive features such as the high parallel scalability and memory-saving of the method in comparison to the existing methods.

The chemical master equation is a continuous time discrete space Markov model of chemical reactions. The chemical master equation is derived mathematically and it is shown that the corresponding initial value problem has a unique solution. Conditions are given under which this solution is a probability distribution. We present finite state and aggregation-disaggregation approximations and provide error bounds for the case of piecewise constant disaggregation. The aggregation-disaggregation approximation allows the solution of the chemical master equation for larger state spaces and is also an important tool for the solution of multidimensional problems. References Adam Arkin, John Ross, and Harley H. McAdams. {Stochastic Kinetic Analysis of Developmental Pathway Bifurcation in Ph age lambda-Infected Escherichia coli Cells}. Genetics, 149(4):1633--1648, 1998. http://www.genetics.org/cgi/content/abstract/149/4/1633 P. Deuflhard, W. Huisinga, T. Jahnke, and M. Wulkow. Adaptive discrete Galerkin methods applied to the chemical master equation. Technical Report ZIB-Report 07-04, Konrad--Zuse Zentrum fuer Informationstechnik Berlin, 2007. http://www.zib.de/Publications/Reports/ZR-07-04.pdf Stefan Engblom. Computing the moments of high dimensional solutions of the master equation. Appl. Math. Comput., 180(2):498--515, 2006. doi:10.1016/j.amc.2005.12.032 Lars Ferm and Per Loetstedt. Numerical method for coupling the macro and meso scales in stochastic chemical kinetics. BIT Numerical Mathematics, 47:735--762, 2007. doi:10.1007/s10543-007-0150-z Daniel T. Gillespie. Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry, 81(25):2340--2361, December 1977. Markus Hegland, Conrad Burden, Lucia Santoso, Shev MacNamara, and Hilary Booth. A solver for the stochastic master equation applied to gene regulatory networks. J. Comput. Appl. Math., 205(2):708--724, 2007. doi:10.1016/j.cam.2006.02.053 Leo Huber and Markus Hegland. Dimension adaptive sparse grids for the chemical master equations. submitted, 2008. B. Munsky and M. Khammash. The finite state projection algorithm for the solution of the chemical master equation. The Journal of Chemical Physics, 124:044104--1--044104--1, 2006. doi:10.1063/1.2145882 Roger B. Sidje. EXPOKIT: Software package for computing matrix exponentials. ACM Transactions on Mathematical Software, 24(1):130--156, March 1998. doi:10.1145/285861.285868 T. Tian and K. Burrage. Bistability and switching in the lysis/lysogeny genetic regulatory network of bacteriophage-$\lambda $. Journal of Theoretical Biology, 227:229--237, 2004. doi:10.1016/j.jtbi.2003.11.003 N. G. van Kampen. Stochastic Processes in Physics and Chemistry. North Holland, Amsterdam, the Netherlands, 1981.

Stochastic modelling and simulations play a major role in Systems Biology because, at molecular level, biological systems exhibit noise coming both from within the cell (intrinsic) and from the environment (extrinsic). Stochastic modelling takes into account the effects of noise over the system dynamics, that can strongly affect the behavior of the system in conditions of relatively low amounts of molecular species. Stochastic simulations provide an effective way to describe the system dynamics, and can be applied on systems where specified chemical species are processed by a set of biochemical reactions, each one characterized by a stochastic constant. In the context of stochastic modelling, Parameter Sweep Applications (PSAs) can be a useful way to explore the huge spaces generated by the combinations of variables and parameters values in order to test their effects on systems dynamics. PSAs are common in the scientific community and are structured as sets of instances, each one characterized by a distinct parametrisation. A PSA that aims to sample such large spaces must involve a large number of instances and hence the problem becomes very time consuming. However, the independence of each instance of a particular PSA makes the distributed computing paradigm a very useful solution for large scale PSAs. In this work we present a grid based version of a multi-volume stochastic simulator, tau-DPP, implemented on the EGEE project platform. The aim of the proposed work is to exploit this platform for testing the reliability of PSAs over the grid, pointing out critical factors, bottlenecks and scalability by providing data about our experience in this kind of biological modelling and simulations. As a case study, we present a number of PSAs for a stochastic model of bacterial chemotaxis composed of 59 reactions and 31 chemical species.

Fixed-time control design for nonlinear uncertain systems via adaptive method

Gillespie's stochastic simulation algorithm (SSA) has been a conventional method for stochastic modeling and simulation of biochemical systems. However, as a population-based algorithm it faces the challenge of combinatorial complexity in many biochemical models where species may present with multiple states. To solve this problem, the rules-based modeling was proposed by Hlavacek's group and the particle-based Network-Free Algorithm (NFA) was used for its simulation. In this paper, we first improve the NFA to efficiently integrate the population-based and particle-based features. Then we propose a population-based SSA scheme for the stochastic simulation of rule-based models. Complexity analysis is presented for the proposed methods. Numerical experiments on two rule-based models demonstrate the power of the proposed methods.

Purpose This purpose of this study is to propose and test an adaptive grid scheme for modeling phase change problems. Design/methodology/approach Appealing to literature related to solving hydrogeology problems and classic moving boundary methodologies, a novel adaptive mesh scheme for Stefan melting problems is proposed. This scheme uses a local adaptive mesh. A set of nodes attached to the moving melt front advances along prescribed directions. All other nodes in the domain retain fixed locations. When the distance of a front node to its nearest nonfront neighboring node exceeds a given distance, a new node is inserted in the grid at a location that preserves a regular grid spacing. Findings When applied to Stefan melting test problems, the proposed scheme provides smooth, accurate and conserved results. Practical implications In addition to the solution of phase change problems, the results from the numerical testing in this work verify the conserved nature of the operation of the proposed adaptive scheme. This indicates that the wider application of the adaptive grid scheme to other conservation problems occurring in domains that evolve over time is valid. Social implications The proposed scheme can be use​d to model hydrogeology problems in evolving landscapes; a problem with important consequences for coastal populations. Originality/value The developed method offers a novel and useful approach for numerically solving problems in time dependent domains.

Stochastic simulations have been done for enzyme kinetics reaction with Michaelis-Menten mechanism in low population number. Gillespie and Poisson algorithms have been used for investigation of population number and fluctuation population around their mean values as a function of time. Our result shows that equilibrium time for population dynamics via Poisson algorithm is smaller than Gillespie algorithm. Variations of average population number versus time for all species have the following order: deterministic approach (mean fields) > Gillespie > Poisson. There is asymptotic limit for fluctuation population as a function of time via Poisson algorithm but there is not such trend for fluctuation population via Gillespie algorithm. There is a maximum for fluctuation population for all species for kinetics reaction with Michaelis-Menten mechanism as a function of time via Gillespie algorithm. The stochastic approach has also been used for horse liver alcohol dehydrogenase which catalyses the NAD \( ^{+} \) (nicotinamide heterocyclic ring) oxidation of ethanol to acetaldehyde and three kinds of third order reactions. Probability distribution function and fluctuation population for reactants are calculated as a function of time. Increasing a variety of species for third order reactions leads to decrease of coefficient variation.

Biochemical systems have numerous practical applications, in particular to the study of critical intracellular processes. Frequently, biochemical kinetic models depict cellular processes as systems of chemical reactions. Many biological processes in a cell are inherently stochastic, due to the existence of some low molecular amounts. These stochastic fluctuations may have a great effect on the biochemical system’s behaviour. In such cases, stochastic models are necessary to accurately describe the system’s dynamics. Biochemical systems at the cellular level may entail many species or reactions and their mathematical models may be non-linear and with multiple scales in time. In this work, we provide a numerical technique for simplifying stochastic discrete models of well-stirred biochemical systems, which ensures that the main properties of the original system are preserved. The proposed technique employs sensitivity analysis and requires solving an optimization problem. The numerical tests on several models of practical interest show that our model reduction strategy performs very well.

Stochastic simulation is an essential method for modeling complex geological structures in geosciences. Evaluating the uncertainty of the realizations of stochastic simulations can better describe real phenomena. However, uncertainty evaluation of stochastic simulation methods remains a challenge due to the limited data from geological surveys and the uncertainty in reliability estimation with stochastic simulation models. In addition, understanding the sensitivity of the parameters in stochastic simulation models is invaluable when exploring the parameters with a higher influence on the uncertainty associated with predictions generated from stochastic simulation. To facilitate uncertainty evaluation in stochastic simulation methods, we use the circular treemap as an interactive workflow to explore prediction uncertainty in and the parameter sensitivity of multiple-point geostatistical (MPS) stochastic simulation methods. In this work, we present a novel visualization framework for assessing the uncertainty in MPS stochastic simulation methods and exploring the parameter sensitivity of the MPS methods. We present a new indicator to integrate multiple metrics that characterize geospatial features and visualize these metrics to assist domain experts in making decisions. Parallel coordinates-scatter matrix plots and multi-dimensional scaling (MDS) plots are used to analyze the parametric sensitivity of MPS stochastic simulation methods. The realizations and parameters of two MPS stochastic simulation methods are used to test the applicability of the proposed visualization workflow and the visualization methods. The results demonstrate that our workflow and the visualization methods can assist experts in finding the model with less uncertainty and improve the efficiency of parameter adjustment using different MPS stochastic simulation methods.

This paper provides a comprehensive review of the state-of-the-art research on stochastic simulations. It stimulates the interest of tackling the problem of stochastic simulation using statistical signal processing methods, as well as innovative thinking of stochastic modeling of gene networks from the viewpoint of signal processing