Search

Stochastic simulation for in silico studies of large biochemical networks requires a great amount of computational time. We recently proposed a new exact simulation algorithm, called the rejection-based stochastic simulation algorithm (RSSA) [Thanh et al., J. Chem. Phys. 141(13), 134116 (2014)], to improve simulation performance by postponing and collapsing as much as possible the propensity updates. In this paper, we analyze the performance of this algorithm in detail, and improve it for simulating large-scale biochemical reaction networks. We also present a new algorithm, called simultaneous RSSA (SRSSA), which generates many independent trajectories simultaneously for the analysis of the biochemical behavior. SRSSA improves simulation performance by utilizing a single data structure across simulations to select reaction firings and forming trajectories. The memory requirement for building and storing the data structure is thus independent of the number of trajectories. The updating of the data structure when needed is performed collectively in a single operation across the simulations. The trajectories generated by SRSSA are exact and independent of each other by exploiting the rejection-based mechanism. We test our new improvement on real biological systems with a wide range of reaction networks to demonstrate its applicability and efficiency.

In the preceding paper [I. Fukuda, J. Chem. Phys. 139, 174107 (2013)], the zero-multipole (ZM) summation method was proposed for efficiently evaluating the electrostatic Coulombic interactions of a classical point charge system. The summation takes a simple pairwise form, but prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large energetic noises and significant artifacts. The purpose of this paper is to judge the ability of the ZM method by investigating the accuracy, parameter dependencies, and stability in applications to liquid systems. To conduct this, first, the energy-functional error was divided into three terms and each term was analyzed by a theoretical error-bound estimation. This estimation gave us a clear basis of the discussions on the numerical investigations. It also gave a new viewpoint between the excess energy error and the damping effect by the damping parameter. Second, with the aid of these analyses, the ZM method was evaluated based on molecular dynamics (MD) simulations of two fundamental liquid systems, a molten sodium-chlorine ion system and a pure water molecule system. In the ion system, the energy accuracy, compared with the Ewald summation, was better for a larger value of multipole moment l currently induced until l ≲ 3 on average. This accuracy improvement with increasing l is due to the enhancement of the excess-energy accuracy. However, this improvement is wholly effective in the total accuracy if the theoretical moment l is smaller than or equal to a system intrinsic moment L. The simulation results thus indicate L ∼ 3 in this system, and we observed less accuracy in l = 4. We demonstrated the origins of parameter dependencies appearing in the crossing behavior and the oscillations of the energy error curves. With raising the moment l we observed, smaller values of the damping parameter provided more accurate results and smoother behaviors with respect to cutoff length were obtained. These features can be explained, on the basis of the theoretical error analyses, such that the excess energy accuracy is improved with increasing l and that the total accuracy improvement within l ⩽ L is facilitated by a small damping parameter. Although the accuracy was fundamentally similar to the ion system, the bulk water system exhibited distinguishable quantitative behaviors. A smaller damping parameter was effective in all the practical cutoff distance, and this fact can be interpreted by the reduction of the excess subset. A lower moment was advantageous in the energy accuracy, where l = 1 was slightly superior to l = 2 in this system. However, the method with l = 2 (viz., the zero-quadrupole sum) gave accurate results for the radial distribution function. We confirmed the stability in the numerical integration for MD simulations employing the ZM scheme. This result is supported by the sufficient smoothness of the energy function. Along with the smoothness, the pairwise feature and the allowance of the atom-based cutoff mode on the energy formula lead to the exact zero total-force, ensuring the total-momentum conservations for typical MD equations of motion.

Combining simulation and model theories, this paper analyses the impact of pair dynamics on the intermolecular nuclear Overhauser effect (NOE) in liquids. For the first time, we give a distance resolved NOE. When applied to the ionic liquid 1-ethyl-3-methyl-imidazolium tetrafluoroborate the NOE turns out to be of long-range nature. This behaviour translates to the experimentally measured cross- and longitudinal relaxation rates. We were able to calculate the heteronuclear NOE from simulation data, despite the high computational effort. Model theories are computationally less demanding and cover the complete frequency range of the respective spectral density function, they are usually based on a very simple pair distribution function and the solution of the diffusion equation. In order to model the simulated data sufficiently, these simplifications in structure and dynamics have to be generalised considerably.

The constrained density functional theory-configuration interaction (CDFT-CI) method has previously been used to calculate ground-state energies and barrier heights, and to describe electronic excited states, in particular conical intersections. However, the method has been limited to evaluating the electronic energy at just a single nuclear configuration, with the gradient of the energy being available only via finite difference. In this paper, we present analytic gradients of the CDFT-CI energy with respect to nuclear coordinates, which gives the potential for accurate geometry optimization and molecular dynamics on both the ground and excited electronic states, a realm which is currently quite challenging for electronic structure theory. We report the performance of CDFT-CI geometry optimization for representative reaction transition states as well as molecules in an excited state. The overall accuracy of CDFT-CI for computing barrier heights is essentially unchanged whether the energies are evaluated at geometries obtained from quadratic configuration-interaction singles and doubles (QCISD) or CDFT-CI, indicating that CDFT-CI produces very good reaction transition states. These results open up tantalizing possibilities for future work on excited states.

In this paper, we propose two methods to carry out the quasi-steady state approximation in stochastic models of enzyme catalytic regulation, based on WKB asymptotics of the chemical master equation or of the corresponding partial differential equation for the generating function. The first of the methods we propose involves the development of multiscale generalisation of a WKB approximation of the solution of the master equation, where the separation of time scales is made explicit which allows us to apply the quasi-steady state approximation in a straightforward manner. To the lowest order, the multi-scale WKB method provides a quasi-steady state, Gaussian approximation of the probability distribution. The second method is based on the Hamilton-Jacobi representation of the stochastic process where, as predicted by large deviation theory, the solution of the partial differential equation for the corresponding characteristic function is given in terms of an effective action functional. The optimal transition paths between two states are then given by those paths that maximise the effective action. Such paths are the solutions of the Hamilton equations for the Hamiltonian associated to the effective action functional. The quasi-steady state approximation is applied to the Hamilton equations thus providing an approximation to the optimal transition paths and the transition time between two states. Using this approximation we predict that, unlike the mean-field quasi-steady approximation result, the rate of enzyme catalysis depends explicitly on the initial number of enzyme molecules. The accuracy and validity of our approximated results as well as that of our predictions regarding the behaviour of the stochastic enzyme catalytic models are verified by direct simulation of the stochastic model using Gillespie stochastic simulation algorithm.

Stochastic modeling is essential for an accurate description of the biochemical network dynamics at the level of a single cell. Biochemically reacting systems often evolve on multiple time-scales, thus their stochastic mathematical models manifest stiffness. Stochastic models which, in addition, are stiff and computationally very challenging, therefore the need for developing effective and accurate numerical methods for approximating their solution. An important stochastic model of well-stirred biochemical systems is the chemical Langevin Equation. The chemical Langevin equation is a system of stochastic differential equation with multidimensional non-commutative noise. This model is valid in the regime of large molecular populations, far from the thermodynamic limit. In this paper, we propose a variable time-stepping strategy for the numerical solution of a general chemical Langevin equation, which applies for any level of randomness in the system. Our variable stepsize method allows arbitrary values of the time-step. Numerical results on several models arising in applications show significant improvement in accuracy and efficiency of the proposed adaptive scheme over the existing methods, the strategies based on halving/doubling of the stepsize and the fixed step-size ones.

Tau leaping methods enable efficient simulation of discrete stochastic chemical systems. Stiff stochastic systems are particularly challenging since implicit methods, which are good for stiffness, result in noninteger states. The occurrence of negative states is also a common problem in tau leaping. In this paper, we introduce the implicit Minkowski-Weyl tau (IMW-τ) methods. Two updating schemes of the IMW-τ methods are presented: implicit Minkowski-Weyl sequential (IMW-S) and implicit Minkowski-Weyl parallel (IMW-P). The main desirable feature of these methods is that they are designed for stiff stochastic systems with molecular copy numbers ranging from small to large and that they produce integer states without rounding. This is accomplished by the use of a split step where the first part is implicit and computes the mean update while the second part is explicit and generates a random update with the mean computed in the first part. We illustrate the IMW-S and IMW-P methods by some numerical examples, and compare them with existing tau methods. For most cases, the IMW-S and IMW-P methods perform favorably.

The so-called "ultrafast" nuclear magnetic resonance (NMR) methods enable the collection of multidimensional spectra within a single scan. These experiments operate by replacing traditional t(1) time increments, with a series of combined radiofrequency-irradiation/magnetic-field-gradient manipulations that spatially encode the effects of the indirect-domain spin interactions. Barring the presence of sizable displacements, the spatial patterns thus imparted can be read out following a mixing period with the aid of oscillating acquisition gradients, leading to a train of t(2)-modulated echoes carrying in their positions and phases the indirect- and the direct-domain spin interactions. Both the initial spatial encoding as well as the subsequent spatial decoding procedures underlying ultrafast NMR were designed under the assumption that spins remain static within the sample during their execution. Most often this is not the case, and motion-related effects can be expected to affect the outcome of these experiments. The present paper focuses on analyzing the effects of diffusion in ultrafast two-dimensional (2D) NMR. Toward this end both analytical and numerical formalisms are derived, capable of dealing with the nonuniform spin manipulations, macroscopic sample sizes, and microscopic displacements involved in this kind of sequences. After experimentally validating the correctness of these formalisms these were used to analyze the effects of diffusion for a variety of cases, including ultrafast experiments on both rapidly and slowly diffusing molecules. A series of prototypical schemes were considered including discrete and continuous encoding modes, constant- and real-time manipulations, homo- and heteronuclear acquisitions, and single versus multiple quantum modalities. The effects of molecular diffusion were also compared against typical relaxation-driven losses as they happen in these various prototypical situations; from all these situations, general guidelines for choosing the optimal ultrafast 2D NMR scheme for a particular sample and condition could be deduced.

The tau-leap method first developed by Gillespie [D. T. Gillespie, J. Chem. Phys. 115, 1716 (2001)] can significantly speed up stochastic simulation of certain chemically reacting systems with acceptable losses in accuracy. Recently, several improved tau-leap methods, including the binomial, multinomial, and modified tau-leap methods, have been developed. However, in all these tau-leap methods, the mean of the number of times, K(m), that the mth reaction channel fires during a leap is not equal to the true mean. Therefore, all existing tau-leap methods produce biased simulation results, which limit the simulation accuracy and speed. In this paper, we analyze the mean of K(m) based on the chemical master equation. Using this analytical result, we develop unbiased Poisson and binomial tau-leap methods. Moreover, we analyze the variance of K(m), and then develop an unbiased Poisson/Gaussian/binomial tau-leap method to correct the errors in both the mean and variance of K(m). Simulation results demonstrate that our unbiased tau-leap method can significantly improve simulation accuracy without sacrificing speed.

We present an extension of the molecular mechanics-valence bond (MMVB) hybrid method to study ground and excited states of planar conjugated hydrocarbon cations. Currently, accurate excited state calculations on these systems are limited to expensive ab initio studies of smaller systems: up to 15 active electrons in 16 pi orbitals with complete active space self-consistent field (CASSCF) theory using high symmetry. The new MMVB extension provides a faster, cheaper treatment to investigate larger cation systems with more than 24 active orbitals. Extension requires both new matrix elements and new parameters: In this paper we present both, for the limited planar case. The scheme is tested for the planar radical cations of benzene, naphthalene, anthracene, and phenanthrene. Calculated MMVB relative energies are in good agreement with CASSCF results for equilibrium geometries on the ground and first excited states, and conical intersections.

The thermodynamics and kinetics of a many-body system can be described in terms of a potential energy landscape in multidimensional configuration space. The partition function of such a landscape can be written in terms of a density of states, which can be computed using a variety of Monte Carlo techniques. In this paper, a new self-consistent Monte Carlo method for computing density of states is described that uses importance sampling and a multiplicative update factor to achieve rapid convergence. The technique is then applied to compute the equilibrium quench probability of the various inherent structures (minima) in the landscape. The quench probability depends on both the potential energy of the inherent structure and the volume of its corresponding basin in configuration space. Finally, the methodology is extended to the isothermal-isobaric ensemble in order to compute inherent structure quench probabilities in an enthalpy landscape.

This paper reports the first computational estimation of the comb polymers' third virial coefficients. The number of the chains in the comb polymers range from 5 to 11. An algorithm that counts the contributing terms of the third virial coefficients in an accelerated manner is presented along with its efficiency dependence on the polymers' size. In addition, the second virial coefficients are estimated for the comb polymers and compared to previously reported results.

Molecular ground-state energies and two-electron reduced density matrices (2-RDMs) have recently been computed without the many-electron wave function by constraining the 2-RDM to satisfy a complete set of three-positivity conditions for N representability [D. A. Mazziotti, Phys. Rev. A 74, 032501 (2006)]. Energies at both equilibrium and nonequilibrium geometries are obtained within 0.3% of the correlation energy. In this paper the authors extend this work to examine the accuracy of molecular properties, including multipole moments and components of the ground-state energy, relative to full configuration interaction (FCI). Comparisons are also made with 2-RDM methods with two-positivity conditions and two-positivity plus the generalized T1T2 conditions as well as several approximate wave function methods. Using the 2-RDM method with three-positivity conditions, the authors obtain dipole, quadrupole, and octupole moments for BeH2, BH, H2O, CO, and NH3 at equilibrium geometries that are within 0.04% of their FCI values. In addition, for the potential energy surface of N2, the 2-RDM method with three-positivity yields not only accurate total ground-state energies but also accurate expectation values of the kinetic energy operator, the electron-nuclei potential, and electron-electron repulsion.

Spatially distributed problems are often approximately modeled in terms of partial differential equations (PDEs) for appropriate coarse-grained quantities (e.g., concentrations). The derivation of accurate such PDEs starting from finer scale, atomistic models, and using suitable averaging is often a challenging task; approximate PDEs are typically obtained through mathematical closure procedures (e.g., mean field approximations). In this paper, we show how such approximate macroscopic PDEs can be exploited in constructing preconditioners to accelerate stochastic computations for spatially distributed particle-based process models. We illustrate how such preconditioning can improve the convergence of equation-free coarse-grained methods based on coarse timesteppers. Our model problem is a stochastic reaction-diffusion model capable of exhibiting Turing instabilities.

The equilibrium free energy difference between two long-lived molecular species or "conformational states" of a protein (or any other molecule) can in principle be estimated by measuring the work needed to shuttle the system between them, independent of the irreversibility of the process. This is the meaning of the Jarzynski equality (JE), which we test in this paper by performing simulations that unfold a protein by pulling two atoms apart. Pulling is performed fast relative to the relaxation time of the molecule and is thus far from equilibrium. Choosing a simple protein model for which we can independently compute its equilibrium properties, we show that the free energy can be exactly and effectively estimated from nonequilibrium simulations. To do so, one must carefully and correctly determine the ensemble of states that are pulled, which is more important the farther from equilibrium one performs simulations; this highlights a potential problem in using the JE to extract the free energy from forced unfolding experiments. The results presented here also demonstrate that the free energy difference between the native and denatured states of a protein measured in solution is not always equal to the free energy profile that can be estimated from forced unfolding simulations (or experiments) using the JE.

Special structures in a chemical reaction network can give rise to bistability, oscillations, and chaos. It has been shown recently [A. Sensse and M. Eiswirth, J. Chem. Phys. 122, 044516 (2005)] that the introduction of an additional species in a supplementary feedback loop to a minimal autocatalytic oscillator gives rise to chaotic dynamics in a certain range of parameters, independent of the particular realization of the additional loop. This provides a possibility to decide if chaos may occur just by analyzing the network structure of an existing model. Here, we apply this concept to analyze the complex dynamics in several essential subsystems of the peroxidase-oxidase reaction system. The aim of the present paper is to determine the nature of the occurring chaos and its location in the parameter space by numerical bifurcation analysis and simulations.

When a minimum on the potential energy surface is surrounded by multiple saddle points with similar energy barriers, the transition pathways with greater prefactors are more important than those that have similar energy barriers but smaller prefactors. In this paper, we present a theoretical formulation for the prefactors, computing the probabilities for transition paths from a minimum to its surrounding saddle points. We apply this formulation to a system of 2 degrees of freedom and a system of 14 degrees of freedom. The first is Brownian motion in a two-dimensional potential whose global anharmonicities play a dominant role in determining the transition rates. The second is a Lennard-Jones (LJ) cluster of seven particles in two dimensions. Low lying transition states of the LJ cluster, which can be reached directly from a minimum without passing through another minimum, are identified without any presumption of their characteristics nor of the product states they lead to. The probabilities are computed for paths going from an equilibrium ensemble of states near a given minimum to the surrounding transition states. These probabilities are directly related to the prefactors in the rate formula. This determination of the rate prefactors includes all anharmonicities, near or far from transition states, which are pertinent in the very sophisticated energy landscape of LJ clusters and in many other complex systems.

The finite-temperature string method proposed by E, et al. [W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66, 052301 (2002)] is a very effective way of identifying transition mechanisms and transition rates between metastable states in systems with complex energy landscapes. In this paper, we discuss the theoretical background and algorithmic details of the finite-temperature string method, as well as the application to the study of isomerization reaction of the alanine dipeptide, both in vacuum and in explicit solvent. We demonstrate that the method allows us to identify directly the isocommittor surfaces, which are approximated by hyperplanes, in the region of configuration space where the most probable transition trajectories are concentrated. These results are verified subsequently by computing directly the committor distribution on the hyperplanes that define the transition state region.

Jarzynski's relation and the fluctuation theorem have established important connections between nonequilibrium statistical mechanics and equilibrium thermodynamics. In particular, an exact relationship between the equilibrium free energy and the nonequilibrium work is useful for computer simulations. In this paper, we exploit the fact that the free energy is a state function, independent of the pathway taken to change the equilibrium ensemble. We show that a generalized expression is advantageous for computer simulations of free energy differences. Several methods based on this idea are proposed. The accuracy and efficiency of the proposed methods are evaluated with a model problem., ((c) 2004 American Institute of Physics.)

It is now widely accepted that there is a class of enzymatic proton transfer reactions, which proceed through quantum tunneling. In a series of papers we have argued that some experimental features of these reactions can be explained by assuming the presence of a "rate-promoting" vibration which brings donor and acceptor closer together, thus leading to rate enhancement. There has never been a study of this effect for classical systems. We used transition path sampling to study the equivalent classical problem and found a complicated dynamical behavior that cannot be captured by transition state theory. Slow promoting vibrations lead to reactive trajectories that overshoot the saddle point, but on the other hand the short period of fast oscillations allows the reactants to stay only briefly in a low-barrier regime. There is a competition between these effects, which results to an intermediate value for the frequency of the rate-promoting vibration that is optimal for enhancing the rate., ((c) 2004 American Institute of Physics)

The semiclassical method is characterized by finite forces and smooth, well-behaved trajectories, but also by multivalued representational functions that are ill behaved at caustics. In contrast, quantum trajectory methods--based on Bohmian mechanics (quantum hydrodynamics)--are characterized by divergent forces and erratic trajectories near nodes, but also well-behaved, single-valued representational functions. In this paper, we unify these two approaches into a single method that captures the best features of both, and in addition, satisfies the correspondence principle. Stationary eigenstates in one degree of freedom are the primary focus, but more general applications are also anticipated., (Copyright 2004 American Institute of Physics)

The determination of reaction paths for enzyme systems remains a great challenge for current computational methods. In this paper we present an efficient method for the determination of minimum energy reaction paths with the ab initio quantum mechanical/molecular mechanical approach. Our method is based on an adaptation of the path optimization procedure by Ayala and Schlegel for small molecules in gas phase, the iterative quantum mechanical/molecular mechanical (QM/MM) optimization method developed earlier in our laboratory and the introduction of a new metric defining the distance between different structures in the configuration space. In this method we represent the reaction path by a discrete set of structures. For each structure we partition the atoms into a core set that usually includes the QM subsystem and an environment set that usually includes the MM subsystem. These two sets are optimized iteratively: the core set is optimized to approximate the reaction path while the environment set is optimized to the corresponding energy minimum. In the optimization of the core set of atoms for the reaction path, we introduce a new metric to define the distances between the points on the reaction path, which excludes the soft degrees of freedom from the environment set and includes extra weights on coordinates describing chemical changes. Because the reaction path is represented by discrete structures and the optimization for each can be performed individually with very limited coupling, our method can be executed in a natural and efficient parallelization, with each processor handling one of the structures. We demonstrate the applicability and efficiency of our method by testing it on two systems previously studied by our group, triosephosphate isomerase and 4-oxalocrotonate tautomerase. In both cases the minimum energy paths for both enzymes agree with the previously reported paths.