Generalized binomial Tau-leap method for biochemical kinetics incorporating both delay and intrinsic noise

The delay stochastic simulation algorithm (DSSA) by Barrio et al. [Plos Comput. Biol.2, 117–E (2006)] was developed to simulate delayed processes in cell biology in the presence of intrinsic noise, that is, when there are small-to-moderate numbers of certain key molecules present in a chemical reaction system. These delayed processes can faithfully represent complex interactions and mechanisms that imply a number of spatiotemporal processes often not explicitly modeled such as transcription and translation, basic in the modeling of cell signaling pathways. However, for systems with widely varying reaction rate constants or large numbers of molecules, the simulation time steps of both the stochastic simulation algorithm (SSA) and the DSSA can become very small causing considerable computational overheads. In order to overcome the limit of small step sizes, various τ-leap strategies have been suggested for improving computational performance of the SSA. In this paper, we present a binomial τ- DSSA method that extends the τ-leap idea to the delay setting and avoids drawing insufficient numbers of reactions, a common shortcoming of existing binomial τ-leap methods that becomes evident when dealing with complex chemical interactions. The resulting inaccuracies are most evident in the delayed case, even when considering reaction products as potential reactants within the same time step in which they are produced. Moreover, we extend the framework to account for multicellular systems with different degrees of intercellular communication. We apply these ideas to two important genetic regulatory models, namely, the hes1 gene, implicated as a molecular clock, and a Her1/Her 7 model for coupled oscillating cells.

Look before you leap : a confidence-based method for selecting species criticality while avoiding negative populations in τ-leaping

The stochastic simulation algorithm was introduced by Gillespie and in a different form by Kurtz. There have been many attempts at accelerating the algorithm without deviating from the behavior of the simulated system. The crux of the explicit τ-leaping procedure is the use of Poisson random variables to approximate the number of occurrences of each type of reaction event during a carefully selected time period, τ. This method is acceptable providing the leap condition, that no propensity function changes “significantly” during any time-step, is met. Using this method there is a possibility that species numbers can, artificially, become negative. Several recent papers have demonstrated methods that avoid this situation. One such method classifies, as critical, those reactions in danger of sending species populations negative. At most, one of these critical reactions is allowed to occur in the next time-step. We argue that the criticality of a reactant species and its dependent reaction channels should be related to the probability of the species number becoming negative. This way only reactions that, if fired, produce a high probability of driving a reactant population negative are labeled critical. The number of firings of more reaction channels can be approximated using Poisson random variables thus speeding up the simulation while maintaining the accuracy. In implementing this revised method of criticality selection we make use of the probability distribution from which the random variable describing the change in species number is drawn. We give several numerical examples to demonstrate the effectiveness of our new method.

Modeling intrinsic noise and delays in chemical kinetics of coupled autoregulated oscillating cells

Delays are an important feature in temporal models of genetic regulation due to slow biochemical processes, such as transcription and translation. In this paper, we show how to model intrinsic noise effects in a delayed setting by either using a delay stochastic simulation algorithm (DSSA) or, for larger and more complex systems, a generalized Binomial τ-leap method (Bτ-DSSA). As a particular application, we apply these ideas to modeling somite segmentation in zebra fish across a number of cells in which two linked oscillatory genes (her1 and her7) are synchronized via Notch signaling between the cells.

Stochastic chemical kinetics and the quasi-steady-state assumption : application to the stochastic simulation algorithm and chemical master equation

Recently the application of the quasi-steady-state approximation (QSSA) to the stochastic simulation algorithm (SSA) was suggested for the purpose of speeding up stochastic simulations of chemical systems that involve both relatively fast and slow chemical reactions [Rao and Arkin, J. Chem. Phys. 118, 4999 (2003)] and further work has led to the nested and slow-scale SSA. Improved numerical efficiency is obtained by respecting the vastly different time scales characterizing the system and then by advancing only the slow reactions exactly, based on a suitable approximation to the fast reactions. We considerably extend these works by applying the QSSA to numerical methods for the direct solution of the chemical master equation (CME) and, in particular, to the finite state projection algorithm [Munsky and Khammash, J. Chem. Phys. 124, 044104 (2006)], in conjunction with Krylov methods. In addition, we point out some important connections to the literature on the (deterministic) total QSSA (tQSSA) and place the stochastic analogue of the QSSA within the more general framework of aggregation of Markov processes. We demonstrate the new methods on four examples: Michaelis–Menten enzyme kinetics, double phosphorylation, the Goldbeter–Koshland switch, and the mitogen activated protein kinase cascade. Overall, we report dramatic improvements by applying the tQSSA to the CME solver.

Stochastic analysis of the VEGF receptor response curve

Recently, an analysis of the response curve of the vascular endothelial growth factor (VEGF) receptor and its application to cancer therapy was described in [T. Alarcón, and K. Page, J. R. Soc. Lond. Interface 4, 283–304 (2007)]. The analysis is significantly extended here by demonstrating that an alternative computational strategy, namely the Krylov FSP algorithm for the direct solution of the chemical master equation, is feasible for the study of the receptor model. The new method allows us to further investigate the hypothesis of symmetry in the stochastic fluctuations of the response. Also, by augmenting the original model with a single reversible reaction we formulate a plausible mechanism capable of realizing a bimodal response, which is reported experimentally but which is not exhibited by the original model. The significance of these findings for mechanisms of tumour resistance to antiangiogenic therapy is discussed.

Accelerated leap methods for simulating discrete stochastic chemical kinetics

Biologists are increasingly conscious of the critical role that noise plays in cellular functions such as genetic regulation, often in connection with fluctuations in small numbers of key regulatory molecules. This has inspired the development of models that capture this fundamentally discrete and stochastic nature of cellular biology - most notably the Gillespie stochastic simulation algorithm (SSA). The SSA simulates a temporally homogeneous, discrete-state, continuous-time Markov process, and of course the corresponding probabilities and numbers of each molecular species must all remain positive. While accurately serving this purpose, the SSA can be computationally inefficient due to very small time stepping so faster approximations such as the Poisson and Binomial τ-leap methods have been suggested. This work places these leap methods in the context of numerical methods for the solution of stochastic differential equations (SDEs) driven by Poisson noise. This allows analogues of Euler-Maruyuma, Milstein and even higher order methods to be developed through the Itô-Taylor expansions as well as similar derivative-free Runge-Kutta approaches. Numerical results demonstrate that these novel methods compare favourably with existing techniques for simulating biochemical reactions by more accurately capturing crucial properties such as the mean and variance than existing methods.

Investigation of phenol degradation : True reaction kinetics on fixed film titanium dioxide photocatalyst

This study of photocatalytic oxidation of phenol over titanium dioxide films presents a method for the evaluation of true reaction kinetics. A flat plate reactor was designed for the specific purpose of investigating the influence of various reaction parameters, specifically photocatalytic film thickness, solution flow rate (1–8 l min−1), phenol concentration (20, 40 and 80 ppm), and irradiation intensity (70.6, 57.9, 37.1and 20.4 W m−2), in order to further understand their impact on the reaction kinetics. Special attention was given to the mass transfer phenomena and the influence of film thickness. The kinetics of phenol degradation were investigated with different irradiation levels and initial pollutant concentration. Photocatalytic degradation experiments were performed to evaluate the influence of mass transfer on the reaction and, in addition, the benzoic acid method was applied for the evaluation of mass transfer coefficient. For this study the reactor was modelled as a batch-recycle reactor. A system of equations that accounts for irradiation, mass transfer and reaction rate was developed to describe the photocatalytic process, to fit the experimental data and to obtain kinetic parameters. The rate of phenol photocatalytic oxidation was described by a Langmuir–Hinshelwood type law that included competitive adsorption and degradation of phenol and its by-products. The by-products were modelled through their additive effect on the solution total organic carbon.

Influence of physical and chemical parameters on the treatment of heavy metals in polluted stormwater using zeolite : A review

Zeolite-based technology can provide a cost effective solution for stormwater treatment for the removal of toxic heavy metals under increasing demand of safe water from alternative sources. This paper reviews the currently available knowledge relating to the effect of properties of zeolites such as pore size, surface area and Si:Al ratio and the physico-chemical conditions of the system such as pH, temperature, initial metal concentration and zeolite concentration on heavy metal removal performance. The primary aims are, to consolidate available knowledge and identify knowledge gaps. It was established that an in-depth understanding of operational issues such as, diffusion of metal ions into the zeolite pore structure, pore clogging, zeolite surface coverage by particulates in stormwater as well as the effect of pH on stormwater quality in the presence of zeolites is essential for developing a zeolite-based technology for the treatment of polluted stormwater. The optimum zeolite concentration to treat typical volumes of stormwater and initial heavy metal concentrations in stormwater should also be considered as operational issues in this regard. Additionally, leaching of aluminium and sodium ions from the zeolite structure to solution were identified as key issues requiring further research in the effort to develop cost effective solutions for the removal of heavy metals from stormwater.

A novel method for the synthesis of monodisperse gold-coated silica nanoparticles

Monodisperse silica nanoparticles were synthesised by the well-known Stober protocol, then dispersed in acetonitrile (ACN) and subsequently added to a bisacetonitrile gold(I) coordination complex ([Au(MeCN)2]?) in ACN. The silica hydroxyl groups were deprotonated in the presence of ACN, generating a formal negative charge on the siloxy groups. This allowed the [Au(MeCN)2]? complex to undergo ligand exchange with the silica nanoparticles and form a surface coordination complex with reduction to metallic gold (Au0) proceeding by an inner sphere mechanism. The residual [Au(MeCN)2]? complex was allowed to react with water, disproportionating into Au0 and Au(III), respectively, with the Au0 adding to the reduced gold already bound on the silica surface. The so-formed metallic gold seed surface was found to be suitable for the conventional reduction of Au(III) to Au0 by ascorbic acid (ASC). This process generated a thin and uniform gold coating on the silica nanoparticles. The silica NPs batches synthesised were in a size range from 45 to 460 nm. Of these silica NP batches, the size range from 400 to 480 nm were used for the gold-coating experiments.

Efficient preconditioning of the method of lines for solving nonlinear two-sided space-fractional diffusion equations

A standard method for the numerical solution of partial differential equations (PDEs) is the method of lines. In this approach the PDE is discretised in space using �finite di�fferences or similar techniques, and the resulting semidiscrete problem in time is integrated using an initial value problem solver. A significant challenge when applying the method of lines to fractional PDEs is that the non-local nature of the fractional derivatives results in a discretised system where each equation involves contributions from many (possibly every) spatial node(s). This has important consequences for the effi�ciency of the numerical solver. First, since the cost of evaluating the discrete equations is high, it is essential to minimise the number of evaluations required to advance the solution in time. Second, since the Jacobian matrix of the system is dense (partially or fully), methods that avoid the need to form and factorise this matrix are preferred. In this paper, we consider a nonlinear two-sided space-fractional di�ffusion equation in one spatial dimension. A key contribution of this paper is to demonstrate how an eff�ective preconditioner is crucial for improving the effi�ciency of the method of lines for solving this equation. In particular, we show how to construct suitable banded approximations to the system Jacobian for preconditioning purposes that permit high orders and large stepsizes to be used in the temporal integration, without requiring dense matrices to be formed. The results of numerical experiments are presented that demonstrate the effectiveness of this approach.

Anomalous subdiffusion equations by the meshless collocation method

This paper aims to develop an implicit meshless collocation technique based on the moving least squares approximation for numerical simulation of the anomalous subdiffusion equation(ASDE). The discrete system of equations is obtained by using the MLS meshless shape functions and the meshless collocation formulation. The stability and convergence of this meshless approach related to the time discretization are investigated theoretically and numerically. The numerical examples with regular and irregular nodal distributions are used to the newly developed meshless formulation. It is concluded that the present meshless formulation is very effective for the modeling of ASDEs.

Alternating direction implicit numerical method for the space and time fractional Bloch-Torrey equation in 3-D

Recently, some authors have considered a new diffusion model–space and time fractional Bloch-Torrey equation (ST-FBTE). Magin et al. (2008) have derived analytical solutions with fractional order dynamics in space (i.e., _ = 1, β an arbitrary real number, 1 < β ≤ 2) and time (i.e., 0 < α < 1, and β = 2), respectively. Yu et al. (2011) have derived an analytical solution and an effective implicit numerical method for solving ST-FBTEs, and also discussed the stability and convergence of the implicit numerical method. However, due to the computational overheads necessary to perform the simulations for nuclear magnetic resonance (NMR) in three dimensions, they present a study based on a two-dimensional example to confirm their theoretical analysis. Alternating direction implicit (ADI) schemes have been proposed for the numerical simulations of classic differential equations. The ADI schemes will reduce a multidimensional problem to a series of independent one-dimensional problems and are thus computationally efficient. In this paper, we consider the numerical solution of a ST-FBTE on a finite domain. The time and space derivatives in the ST-FBTE are replaced by the Caputo and the sequential Riesz fractional derivatives, respectively. A fractional alternating direction implicit scheme (FADIS) for the ST-FBTE in 3-D is proposed. Stability and convergence properties of the FADIS are discussed. Finally, some numerical results for ST-FBTE are given.

Analytical solutions for the two- and three-dimensional time-fractional telegraph equations by the method of separating variables

In this paper, a method of separating variables is effectively implemented for solving a time-fractional telegraph equation (TFTE) in two and three dimensions. We discuss and derive the analytical solution of the TFTE in two and three dimensions with nonhomogeneous Dirichlet boundary condition. This method can be extended to other kinds of the boundary conditions.