Function Of Reaction Coordinate Research Articles

Additive eigenvectors as optimal reaction coordinates, conditioned trajectories, and time-reversible description of stochastic processes.

A fundamental way to analyze complex multidimensional stochastic dynamics is to describe it as diffusion on a free energy landscape-free energy as a function of reaction coordinates (RCs). For such a description to be quantitatively accurate, the RC should be chosen in an optimal way. The committor function is a primary example of an optimal RC for the description of equilibrium reaction dynamics between two states. Here, additive eigenvectors (addevs) are considered as optimal RCs to address the limitations of the committor. An addev master equation for a Markov chain is derived. A stationary solution of the equation describes a sub-ensemble of trajectories conditioned on having the same optimal RC for the forward and time-reversed dynamics in the sub-ensemble. A collection of such sub-ensembles of trajectories, called stochastic eigenmodes, can be used to describe/approximate the stochastic dynamics. A non-stationary solution describes the evolution of the probability distribution. However, in contrast to the standard master equation, it provides a time-reversible description of stochastic dynamics. It can be integrated forward and backward in time. The developed framework is illustrated on two model systems-unidirectional random walk and diffusion.

Computational Investigation on the p53-MDM2 Interaction Using the Potential of Mean Force Study.

Murine double minute 2 (MDM2) proteins are found to be overproduced by many human tumors in order to inhibit the functioning of p53 molecules, a tumor suppressor protein. Thus, reactivating p53 functioning in cancer cells by disrupting p53–MDM2 interactions may offer a significant approach in cancer treatment. However, the structural characterization of the p53–MDM2 complex at the atomistic level and the mechanism of binding/unbinding of the p53–MDM2 complex still remain unclear. Therefore, we demonstrate here the probable binding (unbinding) pathway of transactivation domain 1 of p53 during the formation (dissociation) of the p53–MDM2 complex in terms of free energy as a function of reaction coordinate from the potential of mean force (PMF) study using two different force fields: ff99SB and ff99SB-ILDN. From the PMF plot, we noticed the PMF to have a minimum value at a p53–MDM2 separation of 12 Å, with a dissociation energy of 30 kcal mol–1. We also analyzed the conformational dynamics and stability of p53 as a function of its distance of separation from MDM2. The secondary structure content (helix and turns) in p53 was found to vary with its distance of separation from MDM2. The p53–MDM2 complex structure with lowest potential energy was isolated from the ensemble at the reaction coordinate corresponding to the minimum PMF value and subjected to molecular dynamics simulation to identify the interface surface area, interacting residues at the interface, and the stability of the complex. The simulation results highlight the importance of hydrogen bonds and the salt bridge between Lys94 of MDM2 and Glu17 of p53 in the stability of the p53–MDM2 complex. We also carried out the binding free energy calculations and the per residue energy decomposition analyses of the interface residues of the p53–MDM2 complex. We found that the binding affinity between MDM2 and p53 is indeed high [ΔGbind = −7.29 kcal mol–1 from molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) and ΔGbind = −53.29 kcal mol–1 from molecular mechanics/generalized borne surface area]. The total binding energy obtained using the MM/PBSA method was noticed to be closer to the experimental values (−6.4 to −9.0 kcal mol–1). The p53–MDM2 complex binding profile was observed to follow the same trend even in the duplicate simulation run and also in the simulation carried out with different force fields. We found that Lys51, Leu54, Tyr100, and Tyr104 from MDM2 and the residues Phe19, Trp23, and Leu26 from p53 provide the highest energy contributions for the p53–MDM2 interaction. Our findings highlight the prominent structural and binding characteristics of the p53–MDM2 complex that may be useful in designing potential inhibitors to disrupt the p53–MDM2 interactions.

Protein Folding Free Energy Landscape along the Committor - the Optimal Folding Coordinate.

Recent advances in simulation and experiment have led to dramatic increases in the quantity and complexity of produced data, which makes the development of automated analysis tools very important. A powerful approach to analyze dynamics contained in such data sets is to describe/approximate it by diffusion on a free energy landscape - free energy as a function of reaction coordinates (RC). For the description to be quantitatively accurate, RCs should be chosen in an optimal way. Recent theoretical results show that such an optimal RC exists; however, determining it for practical systems is a very difficult unsolved problem. Here we describe a solution to this problem. We describe an adaptive nonparametric approach to accurately determine the optimal RC (the committor) for an equilibrium trajectory of a realistic system. In contrast to alternative approaches, which require a functional form with many parameters to approximate an RC and thus extensive expertise with the system, the suggested approach is nonparametric and can approximate any RC with high accuracy without system specific information. To avoid overfitting for a realistically sampled system, the approach performs RC optimization in an adaptive manner by focusing optimization on less optimized spatiotemporal regions of the RC. The power of the approach is illustrated on a long equilibrium atomistic folding simulation of HP35 protein. We have determined the optimal folding RC - the committor, which was confirmed by passing a stringent committor validation test. It allowed us to determine a first quantitatively accurate protein folding free energy landscape. We have confirmed the recent theoretical results that diffusion on such a free energy profile can be used to compute exactly the equilibrium flux, the mean first passage times, and the mean transition path times between any two points on the profile. We have shown that the mean squared displacement along the optimal RC grows linear with time as for simple diffusion. The free energy profile allowed us to obtain a direct rigorous estimate of the pre-exponential factor for the folding dynamics.

Multiscale molecular dynamics simulations of rotary motor proteins.

Protein functions require specific structures frequently coupled with conformational changes. The scale of the structural dynamics of proteins spans from the atomic to the molecular level. Theoretically, all-atom molecular dynamics (MD) simulation is a powerful tool to investigate protein dynamics because the MD simulation is capable of capturing conformational changes obeying the intrinsically structural features. However, to study long-timescale dynamics, efficient sampling techniques and coarse-grained (CG) approaches coupled with all-atom MD simulations, termed multiscale MD simulations, are required to overcome the timescale limitation in all-atom MD simulations. Here, we review two examples of rotary motor proteins examined using free energy landscape (FEL) analysis and CG-MD simulations. In the FEL analysis, FEL is calculated as a function of reaction coordinates, and the long-timescale dynamics corresponding to conformational changes is described as transitions on the FEL surface. Another approach is the utilization of the CG model, in which the CG parameters are tuned using the fluctuation matching methodology with all-atom MD simulations. The long-timespan dynamics is then elucidated straightforwardly by using CG-MD simulations.

Interrogation of bimetallic particle oxidation in three dimensions at the nanoscale.

An understanding of bimetallic alloy oxidation is key to the design of hollow-structured binary oxides and the optimization of their catalytic performance. However, one roadblock encountered in studying these binary oxide systems is the difficulty in describing the heterogeneities that occur in both structure and chemistry as a function of reaction coordinate. This is due to the complexity of the three-dimensional mosaic patterns that occur in these heterogeneous binary systems. By combining real-time imaging and chemical-sensitive electron tomography, we show that it is possible to characterize these systems with simultaneous nanoscale and chemical detail. We find that there is oxidation-induced chemical segregation occurring on both external and internal surfaces. Additionally, there is another layer of complexity that occurs during the oxidation, namely that the morphology of the initial oxide surface can change the oxidation modality. This work characterizes the pathways that can control the morphology in binary oxide materials.

Influence of the conserved active site residues of histidyl tRNA synthetase on the mechanism of aminoacylation reaction

The relation between the conservation of active site residues and the molecular mechanism of aminoacylation reaction is an unexplored problem. In the present paper, the influences of the conserved active site residues on the reaction mechanism as well as the electrostatic potential near the reaction center are analyzed for Histidyl tRNA synthetase from Escherichia coli, Thermus thermophilus and Staphylococcus aureus. While the primary structures show both convergence as well as divergence, the secondary level structures of the active sites of the three species show considerable conservation in the respective structural organizations. The conserved active site residues near the reaction center, which have a major role in the reaction mechanism and catalysis, retain their specific position and orientation relative to the substrate in the three species. In order to understand the influence of different conserved and nonconserved residues near the reaction center, two different models are considered. First, a large model of active site with the substrates, Mg 2+ ions and water is constructed in which the first shell residues (including both conserved as well as nonconserved) near the reaction center are studied. From the large model, a smaller model is constructed for reaction path modeling individually for three species. Validation of the smaller model is carried out by comparing the energy surfaces of large and small models as a function of reaction coordinates. Further, the electrostatic potential near the reaction center for the large and small model are compared. The transition state structures of the activation step of aminoacylation reaction for E. coli, T. thermophilus and S. aureus are calculated using the combined ab-initio/semi-empirical calculation. The similarity of the energy profiles as a function of the relevant reaction coordinate and the orientation of the catalytic residue, Arg259, indicate that the reaction mechanisms are identical which are guided by the strikingly similar structural pattern formed by conserved residues for three species. The energy surfaces have close resemblance in three species and present a clear perspective that how the reaction proceeds with the aid of different conserved residues. The study of electrostatic potential confirms this view. The present study provides an understanding of the relationship between the conservation of residues and the efficient reaction mechanism of aminoacylation reaction.

Energy Landscapes Associated with the Self-Aggregation of an Alanine-Based Oligopeptide (AAKA)4

The alanine-based 16-mer peptide (AAKA)(4) has recently been shown to aggregate into a hydrogel at centimolar concentrations and high ionic strength ( Measey, T. J.; Schweitzer-Stenner, R. J. Am. Chem. Soc. 2006 , 128 , 13324. ). This is a surprising result since the current understanding of peptide self-aggregation would lead one to expect that the positively charged lysine residues inhibit the formation of beta-sheet structures. The present study was aimed at shedding some light on the mechanism governing the initial phase of the self-aggregation process. To this end we measured CD and FTIR spectra of the (AAKA)(n) at millimolar concentration and low ionic strength and found that the peptide forms soluble aggregates rather than a hydrogel under these conditions. To analyze the initial phase of this aggregation process, we carried out simulations on the dimerization and trimerization of (AAKA)(4) using replica exchange molecular dynamics (MD) methods based on an implicit solvent model. The results indicate that the peptide aggregates into stable antiparallel beta-sheet structures under conditions comparable to experiments. Furthermore, the existence of trimers is a very sensitive function of the peptide concentration and temperature; no significant amount of trimers is formed at low concentration and the trimer population drops sharply as temperature increases above about 320 K. These results on oligomers are consistent with the trend observed in our experiments (reported elsewhere). We have also obtained the free energy landscapes of the trimers as a function of reaction coordinates, which reveal a few basins of minimal free energy. Associated with them, six major types of stable conformers were identified and the structures of transition states between them were also determined. We suggest a possible transition mechanism between these conformers based on the landscape and structural analysis of the aggregates.

Frontier orbital and density functional study of the variation of the hard–soft behavior of monoborane (BH3) and boron trifluoride (BF3) as a function of angles of reorganization from planar (D3h) to pyramidal (C3v) shape

AbstractIn chemical response the BH3 and BF3 molecules undergo the physical process of planar (D3h) to pyramidal (C3v) reorganization in shape as the condition precedent to the event of chemical reaction under the requirement of symmetry. A frontier orbital and density functional study of the variation of the stability of electronic structures and chemical reactivity of associated with the physical process of D3h to C3v geometry reorganization has been performed. The theoretical parameters viz. eigenvalues of HOMO and LUMO, the HOMO and LUMO energy gap, the global hardness and global softness, the chemical potential, the condensed Fukui function, and local softness of B atom, the reaction site, have been computed over a wide range of ∠XBX angles. The nature of variation in the intrinsic chemical reactivity, global and local, of the molecules associated with their geometry reorganization during the chemical event of charge transfer interaction involving their frontier molecular orbitals has been quantitatively explored. The hardness profiles as a function of reaction coordinates are consistent with the principle of maximum hardness (PMH). Results demonstrate that the hardness and softness are not a static and invariable property of molecules but a dynamic and variable function of molecular structure. The hardness parameters and the HOMO–LUMO gap of the molecules are so modified with the distortion of molecular geometry that, after a certain stage of molecular deformation, the profiles of such parameters of the molecules intersect and cross each other, signifying that the relative order of the intrinsic hardness of their equilibrium geometry is reversed. The intrinsically hard molecule BF3 becomes softer than the intrinsically soft molecule BH3 as a consequence of structural distortion. The increase in chemical reactivity computed in terms of density functional parameters are transparent and justified in terms of the profiles of the eigenvalues of the frontier orbitals. The profiles of chemical potential reveal the inherent difference in the tendency of backdonation from two molecules. The computed values of Fukui functions and local softness parameters of the B atom site demonstrate that the concept of local softness can be exploited for a theoretical analysis and understanding of the characteristic chemical events of the molecules under consideration. The profiles of the Fukui functions and local softness parameters of the two molecules seem to reflect and reveal their intrinsic difference in the tendency of receiving donation in the LUMO (electrophilicity) and that of backdonation from the HOMO (nucleophilicity) and the inherent difference of overall reactivity of the two molecules by a simultaneous operation of two opposing processes of charge transfer. © 2003 Wiley Periodicals, Inc. Int J Quantum Chem, 2003

Canonical Variational Theory for Enzyme Kinetics with the Protein Mean Force and Multidimensional Quantum Mechanical Tunneling Dynamics. Theory and Application to Liver Alcohol Dehydrogenase

We present a theoretical framework for the calculation of rate constants of enzyme-catalyzed reactions that combines variational optimization of the dynamical bottleneck for overbarrier reactive events and multidimensional quantum mechanical tunneling dynamics for through-barrier reactive events, both in the presence of the protein environment. The theory features a two-zone, three-stage procedure called ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) with the transmission coefficient based on the equilibrium secondary-zone (ESZ) approximation for including the effects of the protein on a catalytic reaction center, called the primary zone. The dynamics is calculated by canonical variational theory with optimized multidimensional tunneling contributions, and the formalism allows for Boltzmann averaging over an ensemble of reactant and transition state conformations. In the first stage of the calculations, we assume that the generalized transition states can be well described by a single progress coordinate expressed in primary-zone internal coordinates; in subsequent steps, the transmission coefficient is averaged over a set of primary-zone reaction paths that depend on the protein configuration, and each reaction path has its own reaction coordinate and optimized tunneling path. We also present a simpler approximation to the transmission coefficient that is called the static secondary-zone (SSZ) approximation. We illustrate both versions of this method by carrying out calculations of the reaction rate constants and kinetic isotope effects for oxidation of benzyl alcoholate to benzaldehyde by horse liver alcohol dehydrogenase. The potential energy surface is modeled by a combined generalized hybrid orbital/quantum mechanical/molecular mechanical/semiempirical valence bond (GHO-QM/MM/SEVB) method. The multidimensional tunneling calculations are microcanonically optimized by employing both the small-curvature tunneling approximation and version 4 of the large-curvature tunneling approximation. We find that the variation of the protein mean force as a function of reaction coordinate is quantitatively significant, but it does not change the qualitative conclusions for the present reaction. We obtain good agreement with experiment for both kinetic isotope effects and Swain−Schaad exponents.

Quantum chemical study of the umbrella inversion of the ammonia molecule

The physical event of the umbrella inversion of ammonia has been studied in detail by application of the formalisms of frontier orbital theory, the density functional theory, the localized molecular orbital method, and the energy partitioning analysis. An intuitive structure for the transition state and dynamics of the physical process of structural reorganization prior to inversion have been suggested. The computation starts with the CNDO/2 equilibrium geometry, and thereafter the cycle proceeds through all the conformations of ammonia obtained by varying the ∠HNH angle in steps of 2° from its equilibrium value up to the transition state. The geometry of each conformation is optimized with respect to the length of the N–H bond. The glimpses of the charge density reorganization during the movement of the molecule from equilibrium conformation toward the transition state is computed in terms of dipole moment and the quantum mechanical hybridizations of bond pair and lone pair of N atom through the localized molecular orbitals (LMOs) of all the conformations. Results demonstrate that as the geometry of the molecule begins to evolve through the reorganization of structure, the N–H bond length and the dipole moment begin to decrease, and the trend continues up to the transition state. The dipole moment of the molecule at the suggested transition state is zero. The computed nature of quantum mechanical hybridization of bond pair and lone pair of the N atom as a function of reaction coordinates of the ∠HNH angles reveals that the percentage of s character of the lone pair hybrid decreases and that of the bond pair hybrid forming the σ(N–H) bond increases during the process of geometry reorganization from the equilibrium shape to the transition state. The rationale of the zero dipole moment of the transition state for inversion is not straightforward from its point-group symmetry, but is self-evident from its electronic structure drawn in terms of LMOs. The electronic structure of the transition state, which may be drawn in terms of the LMOs, seems to closely reproduce its suggested intuitive structure. The pattern of variation of dipole moment and nature of the changes of the percentage of the s character in the lone pair hybrid creating dipole with the evolution of geometry during the physical process of structural reorganization for the inversion are found to be nicely correlated according to the suggestion of Coulson. The profiles of the increasing strength of the N–H bond and the increasing percentage of s character of the bond pair hybrid of N atom forming this bond as a function of reaction coordinates are also found to be correlated in accordance with the suggestion of Coulson. The profile of global hardness as a function of reaction coordinate seems to demonstrate that the dynamics of the evolution of the molecular structure from equilibrium shape to the transition state following the course of suggested mode of structural reorganization conforms to the principle of maximum hardness (PMH). The profiles of parameters like the energies of highest occupied and lowest unoccupied molecular orbital (HOMO and LUMO), the gap in energy between HOMO and LUMO, the global hardness, the global softness, and chemical potential as a function of reaction coodinates of a continuous structural evolution from equilibrium shape to the transition state mimic the potential energy diagram of the total energy. Both the frontier orbital parameters and the density functional quantities are found to be equally effective and reliable to monitor the process of necessary structural reorganization prior to the inversion of mofecules. An energy partitioning analysis demonstrates that the origin of barrier has no unique single source rather as many as four mutually exclusive but interacting one- and two-center energy terms within the molecule entail the origin and the height of the barrier. From a close analysis of the results, it seems highly probable that the necessary structural reorganization prior to umbrella inversion of ammonia most realistically occurs following the course of normal modes of vibration of the molecule. © 2000 John Wiley & Sons, Inc. Int J Quant Chem 80: 1–26, 2000

Quantum chemical study of the umbrella inversion of the ammonia molecule

The physical event of the umbrella inversion of ammonia has been studied in detail by application of the formalisms of frontier orbital theory, the density functional theory, the localized molecular orbital method, and the energy partitioning analysis. An intuitive structure for the transition state and dynamics of the physical process of structural reorganization prior to inversion have been suggested. The computation starts with the CNDO/2 equilibrium geometry, and thereafter the cycle proceeds through all the conformations of ammonia obtained by varying the ∠HNH angle in steps of 2° from its equilibrium value up to the transition state. The geometry of each conformation is optimized with respect to the length of the N–H bond. The glimpses of the charge density reorganization during the movement of the molecule from equilibrium conformation toward the transition state is computed in terms of dipole moment and the quantum mechanical hybridizations of bond pair and lone pair of N atom through the localized molecular orbitals (LMOs) of all the conformations. Results demonstrate that as the geometry of the molecule begins to evolve through the reorganization of structure, the N–H bond length and the dipole moment begin to decrease, and the trend continues up to the transition state. The dipole moment of the molecule at the suggested transition state is zero. The computed nature of quantum mechanical hybridization of bond pair and lone pair of the N atom as a function of reaction coordinates of the ∠HNH angles reveals that the percentage of s character of the lone pair hybrid decreases and that of the bond pair hybrid forming the σ(N–H) bond increases during the process of geometry reorganization from the equilibrium shape to the transition state. The rationale of the zero dipole moment of the transition state for inversion is not straightforward from its point-group symmetry, but is self-evident from its electronic structure drawn in terms of LMOs. The electronic structure of the transition state, which may be drawn in terms of the LMOs, seems to closely reproduce its suggested intuitive structure. The pattern of variation of dipole moment and nature of the changes of the percentage of the s character in the lone pair hybrid creating dipole with the evolution of geometry during the physical process of structural reorganization for the inversion are found to be nicely correlated according to the suggestion of Coulson. The profiles of the increasing strength of the N–H bond and the increasing percentage of s character of the bond pair hybrid of N atom forming this bond as a function of reaction coordinates are also found to be correlated in accordance with the suggestion of Coulson. The profile of global hardness as a function of reaction coordinate seems to demonstrate that the dynamics of the evolution of the molecular structure from equilibrium shape to the transition state following the course of suggested mode of structural reorganization conforms to the principle of maximum hardness (PMH). The profiles of parameters like the energies of highest occupied and lowest unoccupied molecular orbital (HOMO and LUMO), the gap in energy between HOMO and LUMO, the global hardness, the global softness, and chemical potential as a function of reaction coodinates of a continuous structural evolution from equilibrium shape to the transition state mimic the potential energy diagram of the total energy. Both the frontier orbital parameters and the density functional quantities are found to be equally effective and reliable to monitor the process of necessary structural reorganization prior to the inversion of mofecules. An energy partitioning analysis demonstrates that the origin of barrier has no unique single source rather as many as four mutually exclusive but interacting one- and two-center energy terms within the molecule entail the origin and the height of the barrier. From a close analysis of the results, it seems highly probable that the necessary structural reorganization prior to umbrella inversion of ammonia most realistically occurs following the course of normal modes of vibration of the molecule. © 2000 John Wiley & Sons, Inc. Int J Quant Chem 80: 1–26, 2000

Generalized born fragment charge model for solvation effects as a function of reaction coordinate

A generalized Born dielectric continuum model is used to calculate the solvation energy of the S N2 reaction Cl − (aq) + CH 3Cl(aq) at three points along the reaction path, using fragment charges and gas-phase energies computed by electronic structure calculations at the MP2/6–31G** level. The results are comparable in accuracy to much more difficult computations involving many-body treatments of the solvent.