K-leap method for accelerating stochastic simulation of coupled chemical reactions

Leap methods are very promising for accelerating stochastic simulation of a well stirred chemically reacting system, while providing acceptable simulation accuracy. In Gillespie's tau-leap method [D. Gillespie, J. Phys. Chem. 115, 1716 (2001)], the number of firings of each reaction channel during a leap is a Poisson random variable, whose sample values are unbounded. This may cause large changes in the populations of certain molecular species during a leap, thereby violating the leap condition. In this paper, we develop an alternative leap method called the K-leap method, in which we constrain the total number of reactions occurring during a leap to be a number K calculated from the leap condition. As the number of firings of each reaction channel during a leap is upper bounded by a properly chosen number, our K-leap method can better satisfy the leap condition, thereby improving simulation accuracy. Since the exact stochastic simulation algorithm (SSA) is a special case of our K-leap method when K=1, our K-leap method can naturally change from the exact SSA to an approximate leap method during simulation, whenever the leap condition allows to do so.     

Accelerated stochastic simulation algorithm for coupled chemical reactions with delays

Some biochemical processes do not occur instantaneously but have considerably delays associated with them. In the existed methods which solve these chemically reacting systems with delays, averaging over a great deal of simulations is needed for reliable statistical characters. Here we present an accelerating approach, called the "Delay Final All Possible Steps" (DFAPS) approach, which does not alter the course of stochastic simulation, but reduces the required running times. Numerical simulation results indicate that the proposed method can be applied to a wide range of chemically reacting systems with delays and obtain a significant improvement on efficiency and accuracy over the existed methods.     

Abstract Next Subvolume Method: a logical process-based approach for spatial stochastic simulation of chemical reactions

The spatial stochastic simulation of biochemical systems requires significant calculation efforts. Parallel discrete-event simulation is a promising approach to accelerate the execution of simulation runs. However, achievable speedup depends on the parallelism inherent in the model. One of our goals is to explore this degree of parallelism in the Next Subvolume Method type simulations. Therefore we introduce the Abstract Next Subvolume Method, in which we decouple the model representation from the sequential simulation algorithms, and prove that state trajectories generated by its executions statistically accord with those generated by the Next Subvolume Method. The experimental performance analysis shows that optimistic synchronization algorithms, together with careful controls over the speculative execution, are necessary to achieve considerable speedup and scalability in parallel spatial stochastic simulation of chemical reactions. Our proposed method facilitates a flexible incorporation of different synchronization algorithms, and can be used to select the proper synchronization algorithm to achieve the efficient parallel simulation of chemical reactions.     

Discrete-time stochastic modeling and simulation of biochemical networks

Since inherent randomness in chemically reacting systems is evident, stochastic modeling and simulation are exceedingly important for investigating complex biological networks. Within the most common stochastic approach a network is modeled by a continuous-time Markov chain governed by the chemical master equation. We show how the continuous-time Markov chain can be converted to a stochastically identical discrete-time Markov chain and obtain a discrete-time version of the chemical master equation. Simulating the discrete-time Markov chain is equivalent to the Gillespie algorithm but requires less effort in that it eliminates the generation of exponential random variables. Thus, exactness as possessed by the Gillespie algorithm is preserved while the simulation can be performed more efficiently.     

Internal stochastic resonance in two coupled chemical oscillators

The dynamics of two coupled chemical oscillators was investigated numerically when the first subsystem was subjected to external parametric noise. The signal-to-noise ratio (SNR) of the response of each subsystem to external noise shows internal stochastic resonance (SR). In addition, the SNR also shows resonance behavior with the variation of coupling strength.     

Deterministic modeling and stochastic simulation of poly-alkoxylation reactions

The main scope of this work is to show the feasibility and the advantage of using a stochastic approach to describe the poly-alkoxylation kinetics of different substrates. For this purpose, the reactions of ethylene and propylene oxides with respectively ethylene glycol, 1-octanol, and 2-octanol were considered. Two kinetic models were used for interpreting all the kinetic runs available in the literature, one deterministic and another one stochastic, for a useful comparison between the two different approaches. As the adopted reaction mechanism, rate laws, and related kinetic parameters were the same for both the kinetic models, the obtained results for what concerns the substrate consumption, and the oligomers distribution profiles were the same in both cases. In the case of the stochastic kinetic approach, the calculations must be made on a small volume of the reacting mixture containing a sufficiently high number of molecules that is suitable for the statistical analysis but as small as possible for reducing the calculation time. The calculations made have allowed individuating this optimal volume. This study is propaedeutic to the application of a stochastic kinetic approach to the study of ethylene-propylene oxides copolymerization that cannot be faced with a deterministic model for the extremely long or impracticable calculation time due to the great number of material balance differential equations that must be integrated.     

Eliminating fast reactions in stochastic simulations of biochemical networks: a bistable genetic switch

In many stochastic simulations of biochemical reaction networks, it is desirable to "coarse grain" the reaction set, removing fast reactions while retaining the correct system dynamics. Various coarse-graining methods have been proposed, but it remains unclear which methods are reliable and which reactions can safely be eliminated. We address these issues for a model gene regulatory network that is particularly sensitive to dynamical fluctuations: a bistable genetic switch. We remove protein-DNA and/or protein-protein association-dissociation reactions from the reaction set using various coarse-graining strategies. We determine the effects on the steady-state probability distribution function and on the rate of fluctuation-driven switch flipping transitions. We find that protein-protein interactions may be safely eliminated from the reaction set, but protein-DNA interactions may not. We also find that it is important to use the chemical master equation rather than macroscopic rate equations to compute effective propensity functions for the coarse-grained reactions.     

State-dependent doubly weighted stochastic simulation algorithm for automatic characterization of stochastic biochemical rare events

In recent years there has been substantial growth in the development of algorithms for characterizing rare events in stochastic biochemical systems. Two such algorithms, the state-dependent weighted stochastic simulation algorithm (swSSA) and the doubly weighted SSA (dwSSA) are extensions of the weighted SSA (wSSA) by H. Kuwahara and I. Mura [J. Chem. Phys. 129, 165101 (2008)]. The swSSA substantially reduces estimator variance by implementing system state-dependent importance sampling (IS) parameters, but lacks an automatic parameter identification strategy. In contrast, the dwSSA provides for the automatic determination of state-independent IS parameters, thus it is inefficient for systems whose states vary widely in time. We present a novel modification of the dwSSA--the state-dependent doubly weighted SSA (sdwSSA)--that combines the strengths of the swSSA and the dwSSA without inheriting their weaknesses. The sdwSSA automatically computes state-dependent IS parameters via the multilevel cross-entropy method. We apply the method to three examples: a reversible isomerization process, a yeast polarization model, and a lac operon model. Our results demonstrate that the sdwSSA offers substantial improvements over previous methods in terms of both accuracy and efficiency.     

Electroanalytical simulations. Orthogonal collocation simulation of fast second-order chemical reactions coupled to an electron transfer with a heterogeneous equivalent formulation

The combination of orthogonal collocation and the heterogeneous equivalent technique is extended to simulate cyclic voltammograms of fast second-order follow-up reactions coupled to an electron transfer at an electrode surface. The $\text{ED}$ (reversible electron transfer with irreversible follow-up dimerization) and $\text{EC}$₂ (reversible electron transfer with irreversible second-order follow-up reaction) models are considered. The non-linear boundary equations are solved numerically. No linear approximation of the concentration profiles is required. The use of non-linear space coordinate transformations is described. Peak potential and peak current function results are compared with literature values and agreement is found. The transition between the second-order $\text{EC}$₂ and the corresponding first-order mechanisms is discussed.     

Thermodynamics of coupled reactions by the amplitudes of chemical relaxation

The amplitudes of the relaxation curves, as obtained by the Temperature-jump method have been used to measure simultaneously equilibrium constant and enthalpy for the reaction of complex formation of Ni²⁺ ion by 2,6-dihydroxobenzoic acid in the presence of a buffer. The experiments have been performed by changing the concentration of metal ion at constant ligand concentration andpH as in a complexometric titration. The points of such ‘dynamic titrations’ have been analysed by means of the concept of ‘normal reactions’ which enabled us to transform a set of coupled individual steps into a set of kinetically independent reactions. The potentialities of the dynamic titrations are discussed.     

Stiffness detection and reduction in discrete stochastic simulation of biochemical systems

Typical multiscale biochemical models contain fast-scale and slow-scale reactions, where "fast" reactions fire much more frequently than "slow" ones. This feature often causes stiffness in discrete stochastic simulation methods such as Gillespie's algorithm and the Tau-Leaping method leading to inefficient simulation. This paper proposes a new strategy to automatically detect stiffness and identify species that cause stiffness for the Tau-Leaping method, as well as two stiffness reduction methods. Numerical results on a stiff decaying dimerization model and a heat shock protein regulation model demonstrate the efficiency and accuracy of the proposed methods for multiscale biochemical systems.     

Domain catalyzed chemical reactions: a molecular dynamics simulation of isomerization kinetics

The model for domain catalyzed isomerization kinetics in condensed fluids is applied for a diluted mixture of a chiral solute with a consolute temperature. The solution is quench to phase separation at temperatures below the consolute temperature. The droplet coalescence enhances the isomerization kinetics due to the substantial excess pressure inside the small droplets given by the Laplace equation. The domain catalyzed isomerization kinetics breaks the symmetry, and the droplets end with only one dominating species. We argue that D-glyceraldehyde which is only moderately solvable in water and which has played a crucial role in the evolution is a candidate for the stereo specific ordering in bio-organic matter.     

Thermodynamic parameters of coupled chemical reactions from temperature jump relaxation amplitudes

Equations describing the temperature jump amplitudes associated with a system of two coupled reactions (TRIS-phenol red) as well as the ternary system (Mg²⁺-iminodiacetic acid-phenol red) are presented. The termodynamic parameters calculated from experimentally measured temperature perturbation amplitudes using a multiparametric curve fitting procedure are found to be in good agreement with those determined from pH- and costant rate thermometric titrations. For phenol red, pK _I =7.55, H _I =3.45 kcal, and for Mg²⁺ iminodiacetic acid, log K _M =2.84, H _M =3.25 kcal, were obtained. It is shown that this method can be used to determine accurate thermodynamic enthalpy changes over a narrow temperature interval of less than 1.0°C from a single experiment requiring about 50 l of sample solution.     

Kinetic Monte Carlo modeling of chemical reactions coupled with heat transfer

In this paper, we describe two types of effective events for describing heat transfer in a kinetic Monte Carlo (KMC) simulation that may involve stochastic chemical reactions. Simulations employing these events are referred to as KMC-TBT and KMC-PHE. In KMC-TBT, heat transfer is modeled as the stochastic transfer of "thermal bits" between adjacent grid points. In KMC-PHE, heat transfer is modeled by integrating the Poisson heat equation for a short time. Either approach is capable of capturing the time dependent system behavior exactly. Both KMC-PHE and KMC-TBT are validated by simulating pure heat transfer in a rod and a square and modeling a heated desorption problem where exact numerical results are available. KMC-PHE is much faster than KMC-TBT and is used to study the endothermic desorption of a lattice gas. Interesting findings from this study are reported.     

Accurate dispensing of volatile reagents on demand for chemical reactions in EWOD chips

Digital microfluidic chips provide a new platform for manipulating chemicals for multi-step chemical synthesis or assays at the microscale. The organic solvents and reagents needed for these applications are often volatile, sensitive to contamination, and wetting, i.e. have contact angles of <90° even on the highly hydrophobic surfaces (e.g., Teflon? or Cytop?) typically used on digital microfluidic chips. Furthermore, often the applications dictate that the processes are performed in a gas environment, not allowing the use of a filler liquid (e.g., oil). These properties pose challenges for delivering controlled volumes of liquid to the chip. An automated, simple, accurate and reliable method of delivering reagents from sealed, off-chip reservoirs is presented here. This platform overcomes the issues of evaporative losses of volatile solvents, cross-contamination, and flooding of the chip by combining a syringe pump, a simple on-chip liquid detector and a robust interface design. The impedance-based liquid detection requires only minimal added hardware to provide a feedback signal to ensure accurate volumes of volatile solvents are introduced to the chip, independent of time delays between dispensing operations. On-demand dispensing of multiple droplets of acetonitrile, a frequently used but difficult to handle solvent due to its wetting properties and volatility, was demonstrated and used to synthesize the positron emission tomography (PET) probe [(18)F]FDG reliably.     

Explicit expression for the relaxation spectrum of a model with infinitely many coupled chemical reactions

The explicit expression for the relaxation spectrum of the polymerization reaction A₁ + A_n ⇌ A_n+1 (n  1, 2, …) where all reactions are supposed to have the same values of the thermodynamic and kinetic parameters, is given. The result deviates significantly from what would be predicted if the reactions were considered to be independent.     

Infrared glass fibers for in-situ sensing, chemical and biochemical reactions

Infrared optical fibres based on chalcogenide glasses have been designed for evanescent wave spectroscopy. The sensitivity of the optical sensor is improved in tapering the sensing zone by chemical etching and the working optical domain of the system has been tested on a chloroform sample. This original remote sensor, based on the analysis of infrared signatures, has been applied to follow the fermentation process in cider fabrication as well as to detect and monitor a bacterial biofilm.     

Transient complex oscillations in a closed chemical system with coupled autocatalysis

In this study, hydroquinone was introduced to the classic Belousov-Zhabotinsky (BZ) reaction to build up coupled autocatalytic feedbacks. Various complex dynamical behaviors including successive period-adding bifurcations, irregular oscillations, and frequency modulations were observed in the coupled reaction system. Not only the complexity of oscillations but also the time period during which complex oscillations persist were found to depend greatly on the initial concentration of hydroquinone, which was expected to manifest the coupling strength in the studied system. Dependence of the observed transient complex oscillations on concentrations of ferroin, sulfuric acid, bromate, and malonic acid was also characterized systematically. Numerical simulations with a modified BZ model via incorporating reactions involving hydroquinone and products of hydroquinone qualitatively reproduced the influence of hydroquinone seen in experiments.