The XPK Package: A Comparison between the Extended Phenomenological Kinetic (XPK) Method and the Conventional Kinetic Monte Carlo (KMC) Method

Recently, we proposed the extended phenomenological kinetics (XPK) method, which overcomes the notorious timescale separation difficulty between fast diffusion and slow chemical reactions in conventional kinetic Monte Carlo (KMC) simulations. In the present work, we make a comprehensive comparison, based on the newly developed XPK package, between the XPK method and the conventional KMC method using a model hydrogenation reaction system. Two potential energy surfaces with different lateral interactions have been designed to illustrate the advantages of the XPK method in computational costs, parallel efficiency and the convergence behaviors to steady states. The XPK method is shown to be efficient and accurate, holding the great promise for theoretical modelling in heterogeneous catalysis, in particular, when the role of the lateral interactions among adsorbates is crucial.     

Molecular simulation of multistate peptide dynamics: a comparison between microsecond timescale sampling and multiple shorter trajectories

Molecular dynamics simulations of the RN24 peptide, which includes a diverse set of structurally heterogeneous states, are carried out in explicit solvent. Two approaches are employed and compared directly under identical simulation conditions. Specifically, we examine sampling by two individual long trajectories (microsecond timescale) and many shorter (MS) uncoupled trajectories. Statistical analysis of the structural properties indicates a qualitative agreement between these approaches. Microsecond timescale sampling gives large uncertainties on most structural metrics, while the shorter timescale of MS simulations results in slight structural memory for beta-structure starting states. Additionally, MS sampling detects numerous transitions on a relatively short timescale that are not observed in microsecond sampling, while long simulations allow for detection of a few transitions on significantly longer timescales. A correlation between the complex free energy landscape and the kinetics of the equilibrium is highlighted by principal component analysis on both simulation sets. This report highlights the increased precision of the MS approach when studying the kinetics of complex conformational change, while revealing the complementary insight and qualitative agreement offered by far fewer individual simulations on significantly longer timescales.     

Army ants algorithm for rare event sampling of delocalized nonadiabatic transitions by trajectory surface hopping and the estimation of sampling errors by the bootstrap method

The most widely used algorithm for Monte Carlo sampling of electronic transitions in trajectory surface hopping (TSH) calculations is the so-called anteater algorithm, which is inefficient for sampling low-probability nonadiabatic events. We present a new sampling scheme (called the army ants algorithm) for carrying out TSH calculations that is applicable to systems with any strength of coupling. The army ants algorithm is a form of rare event sampling whose efficiency is controlled by an input parameter. By choosing a suitable value of the input parameter the army ants algorithm can be reduced to the anteater algorithm (which is efficient for strongly coupled cases), and by optimizing the parameter the army ants algorithm may be efficiently applied to systems with low-probability events. To demonstrate the efficiency of the army ants algorithm, we performed atom-diatom scattering calculations on a model system involving weakly coupled electronic states. Fully converged quantum mechanical calculations were performed, and the probabilities for nonadiabatic reaction and nonreactive deexcitation (quenching) were found to be on the order of 10(-8). For such low-probability events the anteater sampling scheme requires a large number of trajectories ( approximately 10(10)) to obtain good statistics and converged semiclassical results. In contrast by using the new army ants algorithm converged results were obtained by running 10(5) trajectories. Furthermore, the results were found to be in excellent agreement with the quantum mechanical results. Sampling errors were estimated using the bootstrap method, which is validated for use with the army ants algorithm.     

Evolution of conformational changes in the dynamics of small biological molecules: a hybrid MD/RRK approach

The dynamics of long timescale evolution of conformational changes in small biological molecules is described by a hybrid molecular dynamics/RRK algorithm. The approach employs classical trajectories for transitions between adjacent structures separated by a low barrier, and the classical statistical RRK approximation when the barrier involved is high. In determining the long-time dynamics from an initial structure to a final structure of interest, an algorithm is introduced for determining the most efficient pathways (sequence of the intermediate conformers). This method uses the Dijkstra algorithm for finding optimal paths on networks. Three applications of the method using an AMBER force field are presented: a detailed study of conformational transitions in a blocked valine dipeptide; a multiple reaction path study of the blocked valine tripeptide; and the evolution in time from the beta hairpin to alpha helix structure of a blocked alanine hexapeptide. Advantages and limitations of the method are discussed in light of the results.     

Prospects of transition interface sampling simulations for the theoretical study of zeolite synthesis

The transition interface sampling (TIS) technique allows large free energy barriers to be overcome within reasonable simulation time, which is impossible for straightforward molecular dynamics. Still, the method does not impose an artificial driving force, but it surmounts the timescale problem by an importance sampling of true dynamical pathways. Recently, it was shown that the efficiency of TIS when calculating reaction rates is less sensitive to the choice of reaction coordinate than those of the standard free energy based techniques. This could be an important advantage in complex systems for which a good reaction coordinate is usually very difficult to find. We explain the principles of this method and discuss some of the promising applications related to zeolite formation.     

Extending the Sensitivity of CEST NMR Spectroscopy to Micro-to-Millisecond Dynamics in Nucleic Acids Using High-Power Radio-Frequency Fields

Biomolecules undergo motions on the micro-to-millisecond timescale to adopt low-populated transient states that play important roles in folding, recognition, and catalysis. NMR techniques, such as Carr–Purcell–Meiboom–Gill (CPMG), chemical exchange saturation transfer (CEST), and R_1ρ are the most commonly used methods for characterizing such transitions at atomic resolution under solution conditions. CPMG and CEST are most effective at characterizing motions on the millisecond timescale. While some implementations of the R_1ρ experiment are more broadly sensitive to motions on the micro-to-millisecond timescale, they entail the use of selective irradiation schemes and inefficient 1D data acquisition methods. Herein, we show that high-power radio-frequency fields can be used in CEST experiments to extend the sensitivity to faster motions on the micro-to-millisecond timescale. Given the ease of implementing high-power fields in CEST, this should make it easier to characterize micro-to-millisecond dynamics in biomolecules.     

An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks

Stochastic chemical kinetics more accurately describes the dynamics of "small" chemical systems, such as biological cells. Many real systems contain dynamical stiffness, which causes the exact stochastic simulation algorithm or other kinetic Monte Carlo methods to spend the majority of their time executing frequently occurring reaction events. Previous methods have successfully applied a type of probabilistic steady-state approximation by deriving an evolution equation, such as the chemical master equation, for the relaxed fast dynamics and using the solution of that equation to determine the slow dynamics. However, because the solution of the chemical master equation is limited to small, carefully selected, or linear reaction networks, an alternate equation-free method would be highly useful. We present a probabilistic steady-state approximation that separates the time scales of an arbitrary reaction network, detects the convergence of a marginal distribution to a quasi-steady-state, directly samples the underlying distribution, and uses those samples to accurately predict the state of the system, including the effects of the slow dynamics, at future times. The numerical method produces an accurate solution of both the fast and slow reaction dynamics while, for stiff systems, reducing the computational time by orders of magnitude. The developed theory makes no approximations on the shape or form of the underlying steady-state distribution and only assumes that it is ergodic. We demonstrate the accuracy and efficiency of the method using multiple interesting examples, including a highly nonlinear protein-protein interaction network. The developed theory may be applied to any type of kinetic Monte Carlo simulation to more efficiently simulate dynamically stiff systems, including existing exact, approximate, or hybrid stochastic simulation techniques.     

Time-dependent master equation simulation of complex elementary reactions in combustion: application to the reaction of 1CH2 with C2H2 from 300-2000 K

Computational simulations of the title reaction are presented, covering a temperature range from 300 to 2000 K. At lower temperatures we find that initial formation of the cyclopropene complex by addition of methylene to acetylene is irreversible, as is the stabilisation process via collisional energy transfer. Product branching between propargyl and the stable isomers is predicted at 300 K as a function of pressure for the first time. At intermediate temperatures (1200 K), complex temporal evolution involving multiple steady states begins to emerge. At high temperatures (2000 K) the timescale for subsequent unimolecular decay of thermalized intermediates begins to impinge on the timescale for reaction of methylene, such that the rate of formation of propargyl product does not admit a simple analysis in terms of a single time-independent rate constant until the methylene supply becomes depleted. Likewise, at the elevated temperatures the thermalized intermediates cannot be regarded as irreversible product channels. Our solution algorithm involves spectral propagation of a symmetrized version of the discretized master equation matrix, and is implemented in a high precision environment which makes hitherto unachievable low-temperature modelling a reality.     

Quantum dynamics of the femtosecond photoisomerization of retinal in bacteriorhodopsin

The membrane protein bacteriorhodopsin contains all-trans-retinal in a binding site lined by amino acid side groups and water molecules that guide the photodynamics of retinal. Upon absorption of light, retinal undergoes a subpicosecond all-trans-->13-cis phototransformation involving torsion around a double bond. The main reaction product triggers later events in the protein that induce pumping of a proton through bacteriorhodopsin. Quantum-chemical calculations suggest that three coupled electronic states, the ground state and two closely lying excited states, are involved in the motion along the torsional reaction coordinate phi. The evolution of the protein-retinal system on these three electronic surfaces has been modelled using the multiple spawning method for non-adiabatic dynamics. We find that, although most of the population transfer occurs on a timescale of 300 fs, some population transfer occurs on a longer timescale, occasionally extending well beyond 1 ps.     

Atoms of multistationarity in chemical reaction networks

Chemical reaction systems are dynamical systems that arise in chemical engineering and systems biology. In this work, we consider the question of whether the minimal (in a precise sense) multistationary chemical reaction networks, which we propose to call ‘atoms of multistationarity,’ characterize the entire set of multistationary networks. Our main result states that the answer to this question is ‘yes’ in the context of fully open continuous-flow stirred-tank reactors (CFSTRs), which are networks in which all chemical species take part in the inflow and outflow. In order to prove this result, we show that if a subnetwork admits multiple steady states, then these steady states can be lifted to a larger network, provided that the two networks share the same stoichiometric subspace. We also prove an analogous result when a smaller network is obtained from a larger network by ‘removing species.’ Our results provide the mathematical foundation for a technique used by Siegal- Gaskins et al. of establishing bistability by way of ‘network ancestry.’ Additionally, our work provides sufficient conditions for establishing multistationarity by way of atoms and moreover reduces the problem of classifying multistationary CFSTRs to that of cataloging atoms of multistationarity. As an application, we enumerate and classify all 386 bimolecular and reversible two-reaction networks. Of these, exactly 35 admit multiple positive steady states. Moreover, each admits a unique minimal multistationary subnetwork, and these subnetworks form a poset (with respect to the relation of ‘removing species’) which has 11 minimal elements (the atoms of multistationarity).     

Activated Self-Resolution and Error-Correction in Catalytic Reaction Networks**

Understanding the emergence of function in complex reaction networks is a primary goal of systems chemistry and origin-of-life studies. Especially challenging is to create systems that simultaneously exhibit several emergent functions that can be independently tuned. In this work, a multifunctional complex reaction network of nucleophilic small molecule catalysts for the Morita-Baylis-Hillman (MBH) reaction is demonstrated. The dynamic system exhibited triggered self-resolution, preferentially amplifying a specific catalyst/product set out of a many potential alternatives. By utilizing selective reversibility of the products of the reaction set, systemic thermodynamically driven error-correction could also be introduced. To achieve this, a dynamic covalent MBH reaction based on adducts with internal H-transfer capabilities was developed. By careful tuning of the substituents, rate accelerations of retro-MBH reactions of up to four orders of magnitude could be obtained. This study thus demonstrates how efficient self-sorting of catalytic systems can be achieved through an interplay of several complex emergent functionalities.     

Replica exchange with guided annealing for accelerated sampling of disordered protein conformations

We critically examine a recently proposed convective replica exchange (cRE) method for enhanced sampling of protein conformation based on theoretical and numerical analysis. The results demonstrate that cRE and related replica exchange with guided annealing (RE‐GA) schemes lead to unbalanced exchange attempt probabilities and break detailed balance whenever the system undergoes slow conformational transitions (relative to the temperature diffusion timescale). Nonetheless, numerical simulations suggest that approximate canonical ensembles can be generated for systems with small conformational transition barriers. This suggests that RE‐GA maybe suitable for simulating intrinsically disordered proteins, an important class of newly recognized functional proteins. The efficacy of RE‐GA is demonstrated by calculating the conformational ensembles of intrinsically disordered kinase inducible domain protein. The results show that RE‐GA helps the protein to escape nonspecific compact states more efficiently and provides several fold speedups in generating converged and largely correct ensembles compared to the standard temperature RE. © 2014 Wiley Periodicals, Inc.     

Quasiequilibrium approximation of fast reaction kinetics in stochastic biochemical systems

We address the problem of eliminating fast reaction kinetics in stochastic biochemical systems by employing a quasiequilibrium approximation. We build on two previous methodologies developed by [Haseltine and Rawlings, J. Chem. Phys. 117, 6959 (2002)] and by [Rao and Arkin, J. Chem. Phys. 118, 4999 (2003)]. By following Haseltine and Rawlings, we use the numbers of occurrences of the underlying reactions to characterize the state of a biochemical system. We consider systems that can be effectively partitioned into two distinct subsystems, one that comprises "slow" reactions and one that comprises "fast" reactions. We show that when the probabilities of occurrence of the slow reactions depend at most linearly on the states of the fast reactions, we can effectively eliminate the fast reactions by modifying the probabilities of occurrence of the slow reactions. This modification requires computation of the mean states of the fast reactions, conditioned on the states of the slow reactions. By assuming that within consecutive occurrences of slow reactions, the fast reactions rapidly reach equilibrium, we show that the conditional state means of the fast reactions satisfy a system of at most quadratic equations, subject to linear inequality constraints. We present three examples which allow analytical calculations that clearly illustrate the mathematical steps underlying the proposed approximation and demonstrate the accuracy and effectiveness of our method.     

Extending the Sensitivity of CEST NMR Spectroscopy to Micro‐to‐Millisecond Dynamics in Nucleic Acids Using High‐Power Radio‐Frequency Fields

Biomolecules undergo motions on the micro‐to‐millisecond timescale to adopt low‐populated transient states that play important roles in folding, recognition, and catalysis. NMR techniques, such as Carr–Purcell–Meiboom–Gill (CPMG), chemical exchange saturation transfer (CEST), and R_1ρ are the most commonly used methods for characterizing such transitions at atomic resolution under solution conditions. CPMG and CEST are most effective at characterizing motions on the millisecond timescale. While some implementations of the R_1ρ experiment are more broadly sensitive to motions on the micro‐to‐millisecond timescale, they entail the use of selective irradiation schemes and inefficient 1D data acquisition methods. Herein, we show that high‐power radio‐frequency fields can be used in CEST experiments to extend the sensitivity to faster motions on the micro‐to‐millisecond timescale. Given the ease of implementing high‐power fields in CEST, this should make it easier to characterize micro‐to‐millisecond dynamics in biomolecules.     

Advances in Enhanced Sampling Molecular Dynamics Simulations for Biomolecules

Molecular dynamics simulation has emerged as a powerful computational tool for studying biomolecules as it can provide atomic insights into the conformational transitions involved in biological functions. However, when applied to complex biological macromolecules, the conformational sampling ability of conventional molecular dynamics is limited by the rugged free energy landscapes, leading to inherent timescale gaps between molecular dynamics simulations and real biological processes. To address this issue, several advanced enhanced sampling methods have been proposed to improve the sampling efficiency in molecular dynamics. In this review, the theoretical basis, practical applications, and recent improvements of both constraint and unconstrained enhanced sampling methods are summarized. Furthermore, the combined utilizations of different enhanced sampling methods that take advantage of both approaches are also briefly discussed.     

Single-molecule enzyme kinetics in the presence of inhibitors

Recent studies in single-molecule enzyme kinetics reveal that the turnover statistics of a single enzyme is governed by the waiting time distribution that decays as mono-exponential at low substrate concentration and multi-exponential at high substrate concentration. The multi-exponentiality arises due to protein conformational fluctuations, which act on the time scale longer than or comparable to the catalytic reaction step, thereby inducing temporal fluctuations in the catalytic rate resulting in dynamic disorder. In this work, we study the turnover statistics of a single enzyme in the presence of inhibitors to show that the multi-exponentiality in the waiting time distribution can arise even when protein conformational fluctuations do not influence the catalytic rate. From the Michaelis-Menten mechanism of inhibited enzymes, we derive exact expressions for the waiting time distribution for competitive, uncompetitive, and mixed inhibitions to quantitatively show that the presence of inhibitors can induce dynamic disorder in all three modes of inhibitions resulting in temporal fluctuations in the reaction rate. In the presence of inhibitors, dynamic disorder arises due to transitions between active and inhibited states of enzymes, which occur on time scale longer than or comparable to the catalytic step. In this limit, the randomness parameter (dimensionless variance) is greater than unity indicating the presence of dynamic disorder in all three modes of inhibitions. In the opposite limit, when the time scale of the catalytic step is longer than the time scale of transitions between active and inhibited enzymatic states, the randomness parameter is unity, implying no dynamic disorder in the reaction pathway.     

Microsecond timescale MD simulations at the transition state of PmHMGR predict remote allosteric residues

Understanding the mechanisms of enzymatic catalysis requires a detailed understanding of the complex interplay of structure and dynamics of large systems that is a challenge for both experimental and computational approaches. More importantly, the computational demands of QM/MM simulations mean that the dynamics of the reaction can only be considered on a timescale of nanoseconds even though the conformational changes needed to reach the catalytically active state happen on a much slower timescale. Here we demonstrate an alternative approach that uses transition state force fields (TSFFs) derived by the quantum-guided molecular mechanics (Q2MM) method that provides a consistent treatment of the entire system at the classical molecular mechanics level and allows simulations at the microsecond timescale. Application of this approach to the second hydride transfer transition state of HMG-CoA reductase from Pseudomonas mevalonii (PmHMGR) identified three remote residues, R396, E399 and L407, (15–27 Å away from the active site) that have a remote dynamic effect on enzyme activity. The predictions were subsequently validated experimentally via site-directed mutagenesis. These results show that microsecond timescale MD simulations of transition states are possible and can predict rather than just rationalize remote allosteric residues.

Transition state force fields enable MD simulations at the transition state of HMGCoA reductase that sample the transition state ensemble on the μs timescale to identify remote residues that affect the reaction rate.     

Photochemical Z-->E isomerization of a hemithioindigo/hemistilbene omega-amino acid.

The molecule HTI, which combines hemithioindigo and hemistilbene molecular parts, allows reversible switching between two isomeric states. Photochromic behaviour of the HTI molecule is observed by irradiation with UV/Vis light. The photochemical reaction, a Z/E isomerization around the central double bond connecting the two molecular parts, is investigated by transient absorption and emission spectroscopy. For a special HTI molecule, namely, an omega-amino acid, the Z-->E isomerization process occurs on a timescale of 30 ps. In the course of the reaction fast processes on the 1-10 ps timescale are observed which point to motions of the molecule on the potential-energy surface of the excited state. The combination of transient absorption experiments in the visible spectral range with time-resolved fluorescence and infrared measurements reveal a photochemical pathway with three intermediate states. Together with a theoretical modelling procedure the experiments point to a sequential reaction scheme and give indications of the nature of the involved intermediates.     

Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules

Radiationless transitions between electronic excited states in polyatomic molecules take place through unavoided crossings of the potential energy surfaces with substantial non-adiabatic coupling between the respective adiabatic states. While the extent in time of these couplings are large enough, these transitions can be reasonably well simulated through quantum transitions using trajectory surface hopping-like methods. In addition, complex molecular systems may have multiple "trivial" unavoided crossings between noninteracting states. In these cases, the non-adiabatic couplings are described as sharp peaks strongly localized in time. Therefore, their modeling is commonly subjected to the identification of regions close to the particular instantaneous nuclear configurations for which the energy surfaces actually cross each other. Here, we present a novel procedure to identify and treat these regions of unavoided crossings between non-interacting states using the so-called Min-Cost algorithm. The method differentiates between unavoided crossings between interacting states (simulated by quantum hops), and trivial unavoided crossings between non-interacting states (detected by tracking the states in time with Min-Cost procedure). We discuss its implementation within our recently developed non-adiabatic excited state molecular dynamics framework. Fragments of two- and four-ring linear polyphenylene ethynylene chromophore units at various separations have been used as a representative molecular system to test the algorithm. Our results enable us to distinguish and analyze the main features of these different types of radiationless transitions the molecular system undertakes during internal conversion.     

State-selective optimization of local excited electronic states in extended systems

Standard implementations of time-dependent density-functional theory (TDDFT) for the calculation of excitation energies give access to a number of the lowest-lying electronic excitations of a molecule under study. For extended systems, this can become cumbersome if a particular excited state is sought-after because many electronic transitions may be present. This often means that even for systems of moderate size, a multitude of excited states needs to be calculated to cover a certain energy range. Here, we present an algorithm for the selective determination of predefined excited electronic states in an extended system. A guess transition density in terms of orbital transitions has to be provided for the excitation that shall be optimized. The approach employs root-homing techniques together with iterative subspace diagonalization methods to optimize the electronic transition. We illustrate the advantages of this method for solvated molecules, core-excitations of metal complexes, and adsorbates at cluster surfaces. In particular, we study the local π→π(?) excitation of a pyridine molecule adsorbed at a silver cluster. It is shown that the method works very efficiently even for high-lying excited states. We demonstrate that the assumption of a single, well-defined local excitation is, in general, not justified for extended systems, which can lead to root-switching during optimization. In those cases, the method can give important information about the spectral distribution of the orbital transition employed as a guess.