An adaptive stepsize method for the chemical Langevin equation

Mathematical and computational modeling are key tools in analyzing important biological processes in cells and living organisms. In particular, stochastic models are essential to accurately describe the cellular dynamics, when the assumption of the thermodynamic limit can no longer be applied. However, stochastic models are computationally much more challenging than the traditional deterministic models. Moreover, many biochemical systems arising in applications have multiple time-scales, which lead to mathematical stiffness. In this paper we investigate the numerical solution of a stochastic continuous model of well-stirred biochemical systems, the chemical Langevin equation. The chemical Langevin equation is a stochastic differential equation with multiplicative, non-commutative noise. We propose an adaptive stepsize algorithm for approximating the solution of models of biochemical systems in the Langevin regime, with small noise, based on estimates of the local error. The underlying numerical method is the Milstein scheme. The proposed adaptive method is tested on several examples arising in applications and it is shown to have improved efficiency and accuracy compared to the existing fixed stepsize schemes.     

A multiple time step algorithm compatible with a large number of distance classes and an arbitrary distance dependence of the time step size for the fast evaluation of nonbonded interactions in molecular simulations

A new algorithm is introduced to perform the multiple time step integration of the equations of motion for a molecular system, based on the splitting of the nonbonded interactions into a series of distance classes. The interactions between particle pairs in successive classes are updated at a progressively decreasing frequency. Unlike previous multiple time-stepping schemes relying on distance classes, the present algorithm sorts interacting particle pairs by their next update times rather than by their update frequencies. For this reason, the proposed scheme is extremely flexible with respect to the number of classes that can be employed (up to hundred or more) and the distance dependence of the relative time step size (arbitrary integer function of the distance). It can also easily be adapted to classes defined based on a criterion other than the interparticle distance (e.g., interaction magnitude). Different variants of the algorithm are tested in terms of accuracy and efficiency for simulations of a pure water system (6167 molecules) under truncated-octahedral periodic boundary conditions, and compared to the twin-range method standardly used with GROMOS96 (short- and long-range cutoff distances of 0.8 and 1.4 nm, pair list and intermediate-range interactions updated every five steps). In particular, multiple time-stepping schemes with an accuracy comparable to that of the twin-range method can be designed, that permit to increase the effective (long-range) cutoff distance from 1.4 to 3.0 nm with a performance loss of only about a factor 2. This result is quite encouraging, considering the benefits of doubling the cutoff radius in the context of (bio-)molecular simulations.     

Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.     

A computational algorithm for the Green's function method of sensitivity analysis in chemical kinetics

The recent interest in numerical modeling of chemical kinetics has generated the need for proper analysis of the system sensitivities in such models. This paper describes the logic for a program developed by the authors to implement the Green's function method of sensitivity analysis in complex kinetic schemes. The relevant equations and numerical details of the algorithm are outlined, two flow charts are provided, and some special programming considerations are discussed in some detail. Computer storage and computational time considerations are also treated. Finally, applications of sensitivity information to understanding complex kinetic system behavior and analyzing experimental results are suggested.     

Locally implicit solution of a reaction-diffusion system with stiff kinetics

The Tyson-Fife reaction-diffusion equations are solved numerically using a locally implicit approach. Since the variables evolve at very different time scales, the resulting system of equations is stiff. The reaction term is responsible for the stiffness and the time step is increased by using an implicit method. The diffusion operator is evaluated explicitly and the system of implicit nonlinear equations is decoupled. The method is particularly useful for parameter values in which the equations are very stiff, such as the values obtained directly from the experimental reaction rate constants. Previous efforts modified the parameters on the equations to avoid stiffness. The equations then become a simplified model of excitable media and, for those cases, the locally implicit method gives a faster although less accurate solution. Nevertheless, since the modified equations no longer represent a particular chemical system an accurate solution is not as important. The algorithm is applied to observe the transition from simple motion to compound motion of a spiral tip.     

Variable step molecular dynamics: An exploratory technique for peptides with fixed geometry

A new algorithm is presented for performing molecular dynamics simulations of peptides with fixed geometry, with the aim of simulating conformational changes and of exploring conformational space. The principle of the method is to expand the potential energy as a Taylor's series in the coordinates around the current point, retaining the force and its first two derivatives, and obtain a series solution of the resulting differential equations using a method due to Lyapunov. By choosing the time step so that the second term in the series is small compared to the first, the true solution can in principle be approximated to any desired degree of accuracy. The algorithm has been used to solve numerically Lagrange's equations of motion for N-acetyl alanine amide and N-acetyl methionide amide, regarded as fixed at their C-termini, under the influence of the ECEPP/2 potential energy function, and time steps of 15–30 fsec have been achieved with little variation in the total energy. Possible directions for future development are discussed.     

Implementation of a symplectic multiple-time-step molecular dynamics algorithm, based on the united-residue mesoscopic potential energy function

A symplectic multiple-time-step (MTS) algorithm has been developed for the united-residue (UNRES) force field. In this algorithm, the slow-varying forces (which contain most of the long-range interactions and are, therefore, expensive to compute) are integrated with a larger time step, termed the basic time step, and the fast-varying forces are integrated with a shorter time step, which is an integral fraction of the basic time step. Based on the split operator formalism, the equations of motion were derived. Separation of the fast- and slow-varying forces leads to stable molecular dynamics with longer time steps. The algorithms were tested with the Ala(10) polypeptide chain and two versions of the UNRES force field: the current one in which the energy components accounting for the energetics of side-chain rotamers (U(rot)) can lead to numerically unstable forces and a modified one in which the the present U(rot) was replaced by a numerically stable expression which, at present, is parametrized only for polyalanine chains. With the modified UNRES potential, stable trajectories were obtained even when extending the basic time step to 15 fs and, with the original UNRES potentials, the basic time step is 1 fs. An adaptive multiple-time-step (A-MTS) algorithm is proposed to handle instabilities in the forces; in this method, the number of substeps in the basic time step varies depending on the change of the magnitude of the acceleration. With this algorithm, the basic time step is 1 fs but the number of substeps and, consequently, the computational cost are reduced with respect to the MTS algorithm. The use of the UNRES mesoscopic energy function and the algorithms derived in this work enables one to increase the simulation time period by several orders of magnitude compared to conventional atomic-resolution molecular dynamics approaches and, consequently, such an approach appears applicable to simulating protein-folding pathways, protein functional dynamics in a real molecular environment, and dynamical molecular recognition processes.     

Modeling continuous aqueous two-phase systems for control purposes

A mathematical model was proposed to allow the analysis of steady-state and transient behaviors of single-stage continuous aqueous two-phase systems. Since the complete system of simultaneous equations contains more equations than unknown variables, a program based on the method of least squares was developed to solve the problem. The methodology was tested using a system composed of thaumatin, sodium chloride, and a contaminant protein. A poly(ethylene-glycol)/phosphate salt/water system was selected to isolate the thaumatin. For the steady-state and transient operations, a constrained optimization procedure--from the Matlab Optimization Toolbox (MathWork, Inc.)--was implemented after recasting the system of equations as a minimization problem. Euler's method was used in the transient case to discretize the differential equations. The steady-state concentrations agreed with published data. An input-output model based on a 4% step change decrease in the inlet stream flow rate showed that output variables such as concentrations of sodium chloride and phosphate salt settled to their final values in different time periods. The proposed analysis may be helpful in the dynamic control of large-scale commercial extractor units using advanced control schemes.     

Numerical Modeling of Non-Steady-State Ion Transfer in Electrochemical Systems with Allowance for Migration

Numerical methods that are used for modeling non-steady-state ion transfer in electrochemical systems and account for the diffusion, migration, convection, and homogeneous chemical reactions are analyzed. It is shown that the violation of the electroneutrality condition (ENC) in the process of numerical solution is due to the difference equations being inconsistent with the initial differential equations. Difference schemes for numerical calculation of transfer processes, which make it possible to split a set of coupled equations, are designed and conditions for their stability are determined. The explicit difference scheme is self-consistent, i.e. it ensures that ENC is rigorously met. In the implicit difference scheme, ENC is probably violated when splitting the set of equations. To restore electroneutrality of the medium, it is proposed to use a physically substantiated analytical relation for the space charge relaxation under the action of a strong electric field.     

Finite difference method of simulation of non-steady-state ion transfer in electrochemical systems with allowance for migration

Finite difference methods of the second order of accuracy are elaborated for numerical calculation of non-steady-state ion transfer, which is caused by diffusion, migration, and convection in the unidimensional electrochemical systems. The methods of decoupling a set of coupled continuity equations of the electrolyte species are proposed, which ensures that the discrete equations are consistent with the initial differential equations and the electroneutrality condition is rigorously met. The methods of approximation of the boundary conditions of the second order temporal and spatial accuracy and the method of decoupling the transfer equations in the boundary nodes are elaborated. The explicit, fully implicit, and semi-implicit finite difference schemes are elaborated. For semi-implicit schemes, two versions of difference equation closure are proposed, which assure the unambiguity of determination of the distribution of electrical potential. Comparison analysis of the accuracy of elaborated finite difference methods of calculation of non-steady-state ion transfer is performed.     

An automated method for graph-based chemical space exploration and transition state finding

Algorithms that automatically explore the chemical space have been limited to chemical systems with a low number of atoms due to expensive involved quantum calculations and the large amount of possible reaction pathways. The method described here presents a novel solution to the problem of chemical exploration by generating reaction networks with heuristics based on chemical theory. First, a second version of the reaction network is determined through molecular graph transformations acting upon functional groups of the reacting. Only transformations that break two chemical bonds and form two new ones are considered, leading to a significant performance enhancement compared to previously presented algorithm. Second, energy barriers for this reaction network are estimated through quantum chemical calculations by a growing string method, which can also identify non-octet species missed during the previous step and further define the reaction network. The proposed algorithm has been successfully applied to five different chemical reactions, in all cases identifying the most important reaction pathways.     

A family of parametric schemes of arbitrary even order for solving nonlinear models: CMMSE2016

Many problems related to gas dynamics, heat transfer or chemical reactions are modeled by means of partial differential equations that usually are solved by using approximation techniques. When they are transformed in nonlinear systems of equations via a discretization process, this system is big-sized and high-order iterative methods are specially useful. In this paper, we construct a new family of parametric iterative methods with arbitrary even order, based on the extension of Ostrowski’ and Chun’s methods for solving nonlinear systems. Some elements of the proposed class are known methods meanwhile others are new schemes with good properties. Some numerical tests confirm the theoretical results and allow us to compare the numerical results obtained by applying new methods and known ones on academical examples. In addition, we apply one of our methods for approximating the solution of a heat conduction problem described by a parabolic partial differential equation.     

Reduction and solution of the chemical master equation using time scale separation and finite state projection

The dynamics of chemical reaction networks often takes place on widely differing time scales--from the order of nanoseconds to the order of several days. This is particularly true for gene regulatory networks, which are modeled by chemical kinetics. Multiple time scales in mathematical models often lead to serious computational difficulties, such as numerical stiffness in the case of differential equations or excessively redundant Monte Carlo simulations in the case of stochastic processes. We present a model reduction method for study of stochastic chemical kinetic systems that takes advantage of multiple time scales. The method applies to finite projections of the chemical master equation and allows for effective time scale separation of the system dynamics. We implement this method in a novel numerical algorithm that exploits the time scale separation to achieve model order reductions while enabling error checking and control. We illustrate the efficiency of our method in several examples motivated by recent developments in gene regulatory networks.     

Design and parameterization of a stochastic cellular automaton describing a chemical reaction

Although most of the work concerned with reaction kinetics concentrates on empirical findings, stochastic models, and differential equations, a growing number of researchers is exploring other methods to elucidate reaction kinetics. In this work, the parameterization of an utter discrete spatio-temporal model, more specifically, a cellular automaton (CA), describing the reaction of HCl with CaCO(3) , is suggested. Furthermore, a system of partial differential equations (PDE), deduced from a set of CA rules, is implemented to compare both modeling paradigms. In this article, the experimental setup to acquire time series of data is explained, a stochastic CA-based model and a continuous PDE-based model capable of describing the reaction are proposed, the models are parameterized using the experimental data and, finally, the relationship between a discrete time step of the CA-based model and the physical time is studied. Essentially, the parameterization of both models can be traced back to the quest for a solution of the inverse problem in which a (set of) rule(s), respectively a system of PDE, is deduced starting from the observed data. It is demonstrated that the proposed CA- and PDE-based models are capable of describing the considered chemical reaction with a high accuracy, which is confirmed by a root mean squared error between the simulated and observed data of 0.388 and 0.869 g CO(2) , respectively. Further, it is shown that an exponential or linear relationship can be used to link the physical time to a discrete time step of the CA-based model.     

Improved diffusion Monte Carlo propagators for bosonic systems using Itô calculus

The construction of importance sampled diffusion Monte Carlo (DMC) schemes accurate to second order in the time step is discussed. A central aspect in obtaining efficient second order schemes is the numerical solution of the stochastic differential equation (SDE) associated with the Fokker-Plank equation responsible for the importance sampling procedure. In this work, stochastic predictor-corrector schemes solving the SDE and consistent with It? calculus are used in DMC simulations of helium clusters. These schemes are numerically compared with alternative algorithms obtained by splitting the Fokker-Plank operator, an approach that we analyze using the analytical tools provided by Ito; calculus. The numerical results show that predictor-corrector methods are indeed accurate to second order in the time step and that they present a smaller time step bias and a better efficiency than second order split-operator derived schemes when computing ensemble averages for bosonic systems. The possible extension of the predictor-corrector methods to higher orders is also discussed.     

Parameterization of a reactive force field using a Monte Carlo algorithm

Parameterization of a molecular dynamics force field is essential in realistically modeling the physicochemical processes involved in a molecular system. This step is often challenging when the equations involved in describing the force field are complicated as well as when the parameters are mostly empirical. ReaxFF is one such reactive force field which uses hundreds of parameters to describe the interactions between atoms. The optimization of the parameters in ReaxFF is done such that the properties predicted by ReaxFF matches with a set of quantum chemical or experimental data. Usually, the optimization of the parameters is done by an inefficient single‐parameter parabolic‐search algorithm. In this study, we use a robust metropolis Monte‐Carlo algorithm with simulated annealing to search for the optimum parameters for the ReaxFF force field in a high‐dimensional parameter space. The optimization is done against a set of quantum chemical data for MgSO₄ hydrates. The optimized force field reproduced the chemical structures, the equations of state, and the water binding curves of MgSO₄ hydrates. The transferability test of the ReaxFF force field shows the extend of transferability for a particular molecular system. This study points out that the ReaxFF force field is not indefinitely transferable. © 2013 Wiley Periodicals, Inc.     

Comparison of compuational methods for simulating nuclear Overhauser effects in NMR spectroscopy

Various algorithms for solving the Solomon equations describing nuclear Overhauser effects (nOes) in NMR spectroscopy have been compared. The applicability of the eigenvalue/eigenvector and the numerical integration approaches have been investigated. The eigenvalue/eigenvector approach is not a computationally efficient means of simulating nOe experiments in which a saturating radiofrequency field is applied during the time course. For experiments in which nOes develop in the absence of an RF field, this approach should only be used in simulating a full NOESY spectrum. Integration schemes have been found to be more efficient at simulating nOe experiments in which the nOe evolves in the presence of a saturating field, at simulating a partial set of initial perturbation experiments and at simulating a few rows or columns in a NOESY spectrum. Various integration schemes were applied to a two-spin system for which an analytic solution is available and to a model B-DNA oligonucleotide hexamer. The previously unused Taylor series algorithm was found to be superior to the Euler, midpoint, and fourth-order Runge–Kutta methods with regard to integration accuracy/computation time. An adaptive step size control routine for the Taylor series integration scheme was developed. Integration schemes can be speeded up in a simple fashion by introducing a distance cutoff for the dipolar interaction. Using a cutoff of 8 Å the Taylor series algorithm was able to compute the NOESY spectrum more rapidly than the eigenvalue/eigenvector algorithm for large spin systems at short mixing times. At longer mixing times the eigenvalue/eigenvector approach becomes the more efficient scheme.     

Developing Itô stochastic differential equation models for neuronal signal transduction pathways

Mathematical modeling and simulation of dynamic biochemical systems are receiving considerable attention due to the increasing availability of experimental knowledge of complex intracellular functions. In addition to deterministic approaches, several stochastic approaches have been developed for simulating the time-series behavior of biochemical systems. The problem with stochastic approaches, however, is the larger computational time compared to deterministic approaches. It is therefore necessary to study alternative ways to incorporate stochasticity and to seek approaches that reduce the computational time needed for simulations, yet preserve the characteristic behavior of the system in question. In this work, we develop a computational framework based on the It? stochastic differential equations for neuronal signal transduction networks. There are several different ways to incorporate stochasticity into deterministic differential equation models and to obtain It? stochastic differential equations. Two of the developed models are found most suitable for stochastic modeling of neuronal signal transduction. The best models give stable responses which means that the variances of the responses with time are not increasing and negative concentrations are avoided. We also make a comparative analysis of different kinds of stochastic approaches, that is the It? stochastic differential equations, the chemical Langevin equation, and the Gillespie stochastic simulation algorithm. Different kinds of stochastic approaches can be used to produce similar responses for the neuronal protein kinase C signal transduction pathway. The fine details of the responses vary slightly, depending on the approach and the parameter values. However, when simulating great numbers of chemical species, the Gillespie algorithm is computationally several orders of magnitude slower than the It? stochastic differential equations and the chemical Langevin equation. Furthermore, the chemical Langevin equation produces negative concentrations. The It? stochastic differential equations developed in this work are shown to overcome the problem of obtaining negative values.     

含时密度矩阵重正化群的理论与应用

密度矩阵重正化群(DMRG)作为计算低维强关联体系强有力的方法为人熟知,在量子化学电子结构计算中得到广泛应用.最近几年,含时密度矩阵重正化群(TD-DMRG)的理论取得较快发展, TD-DMRG逐渐成为复杂体系量子动力学理论模拟的重要新兴方法之一.本文综述了基于矩阵乘积态(MPS)和矩阵乘积算符(MPO)的DMRG基本理论,并重点介绍了若干最常见的TD-DMRG时间演化算法,包括基于演化再压缩(P&C)的算法、基于含时变分原理(TDVP)的算法和时间步瞄准(TST)算法;还对利用TD-DMRG模拟有限温体系的纯化(Purification)算法和最小纠缠典型量子热态(METTS)算法进行了介绍.最后,对近年TD-DMRG在复杂体系量子动力学中的应用进行了总结.     

Correlation regions within a localized molecular orbital approach

A hybrid scheme for the computation of reaction energies in large molecular systems is proposed. The approach is based on localized orbitals and allows for the treatment of different parts of a molecule at different computational levels. The localized orbitals are assigned to regions, and then different local correlation methods, such as local MP2 or local CCSD(T), can be applied to different regions. In contrast to previous hybrid schemes, the molecule does not have to be split into parts and, therefore, it is not necessary to saturate dangling bonds using link atoms. For fixed region sizes, the cost of the high-level calculation becomes independent of the molecular size, and it is demonstrated that O(1) scaling can be achieved. Illustrative applications are presented and the convergence of the results with respect to the size of the regions is investigated for reaction energies, barrier heights, and weakly bound complexes.