Normal mode analysis using the driven molecular dynamics method. II. An application to biological macromolecules

Rare Events Sampling Methods for Quantum and Classical Ab Initio Molecular Dynamics

We provide an approach to sample rare events during classical ab initio molecular dynamics and quantum wavepacket dynamics. For classical AIMD, a set of fictitious degrees of freedom are introduced that may harmonically interact with the electronic and nuclear degrees of freedom to steer the dynamics in a conservative fashion toward energetically forbidden regions. A similar approach when introduced for quantum wavepacket dynamics has the effect of biasing the trajectory of the wavepacket centroid toward the regions of the potential surface that are difficult to sample. The approach is demonstrated for a phenol-amine system, which is a prototypical problem for condensed phase-proton transfer, and for model potentials undergoing wavepacket dynamics. In all cases, the approach yields trajectories that conserve energy while sampling rare events.

Permutationally invariant polynomial representation of polarizability tensor surfaces for linear regression analysis

A linearly parameterized functional form for a Cartesian representation of molecular dipole polarizability tensor surfaces (PTS) is described. The proposed expression for the PTS is a linearization of the recently reported power series ansatz of the original Applequist model, which by construction is non-linear in parameter space. This new approach possesses (i) a unique solution to the least-squares fitting problem; (ii) a low level of the computational complexity of the resulting linear regression procedure, comparable to those of the potential energy and dipole moment surfaces; and (iii) a competitive level of accuracy compared to the non-linear PTS model. Calculations of CH₄ PTS, with polarizabilities fitted to 9000 training set points with the energies up to 14,000 cm^-1 show an impressive level of accuracy of the linear PTS model obtained with ~1600 parameters: ~1% versus 0.3% RMSE for the non-linear vs. linear model on a test set of 1000 configurations.

Some opinions on MD-based vibrational spectroscopy of gas phase molecules and their assembly: An overview of what has been achieved and where to go

We hereby review molecular dynamics simulations for anharmonic gas phase spectroscopy and provide some of our opinions of where the field is heading. With these new directions, the theoretical IR/Raman spectroscopy of large (bio)-molecular systems will be more easily achievable over longer time-scale MD trajectories for an increase in accuracy of the MD-IR and MD-Raman calculated spectra. With the new directions presented here, the high throughput 'decoding' of experimental IR/Raman spectra into 3D-structures should thus be possible, hence advancing e.g. the field of MS-IR for structural characterization by spectroscopy. We also review the assignment of vibrational spectra in terms of anharmonic molecular modes from the MD trajectories, and especially introduce our recent developments based on Graph Theory algorithms. Graph Theory algorithmic is also introduced in this review for the identification of the molecular 3D-structures sampled over MD trajectories.

Analysis of the Proton Transfer Bands in the Infrared Spectra of Linear N2H+···OC and N2D+···OC Complexes Using Electric Field-Driven Classical Trajectories

In this work, we describe ab initio calculations and assignment of infrared (IR) spectra of hydrogen-bonded ion-molecular complexes that involve a fluxional proton: the linear N₂H⁺···OC and N₂D⁺···OC complexes. Given the challenges of describing fluxional proton dynamics and especially its IR activity, we use electric field-driven classical trajectories, i.e., the driven molecular dynamics (DMD) method that was developed by us in recent years and for similar applications, in conjunction with high-level electronic structure theory. Namely, we present a modified and a numerically efficient implementation of DMD specifically for direct (or "on the fly") calculations, which we carry out at the MP2-F12/AVDZ level of theory for the potential energy surface (PES) and MP2/AVDZ for the dipole moment surfaces (DMSs). Detailed analysis of the PES, DMS, and the time-dependence of the first derivative of the DMS, referred to as the driving force, for the highly fluxional vibrations involving H⁺/D⁺ revealed that the strongly non-harmonic PES and non-linear DMS yield remarkably complex vibrational spectra. Interestingly, the classical trajectories reveal a doublet in the proton transfer part of the spectrum with the two peaks at 1800 and 1980 cm^-1. We find that their shared intensity is due to a Fermi-like resonance interaction, within the classical limit, of the H⁺ parallel stretch fundamental and an H⁺ perpendicular bending overtone. This doublet is also observed in the deuterated species at 1360 and 1460 cm^-1.

Prediction of the terahertz absorption features with a straightforward molecular dynamics method

In this paper, we provide a straightforward method to predict the terahertz absorption spectrum based on a fixed charge model with classic molecular dynamics calculations. The absorption features in the frequency range between 1 and 3.4 THz of stearic acid B-form and between 1 and 2.7 THz of C-form were successfully calculated. Most of the absorption peaks from the simulation correspond well with those from the measurements. By calculating the spatial and time-dependent energy accumulation in the molecular system, the core idea of our calculation method is further validated. Compared with the ab initio calculations, our method provides a computationally inexpensive way to accurately predict the locations of absorption features. With regard to the traditional molecular dynamic simulations, our method is able to extract the spatial distribution of the energy accumulation as well as the local motions in the molecular system.

Gas-Phase Infrared Spectroscopy of Neutral Peptides: Insights from the Far-IR and THz Domain

Gas-phase, double resonance IR spectroscopy has proven to be an excellent approach to obtain structural information on peptides ranging from single amino acids to large peptides and peptide clusters. In this review, we discuss the state-of-the-art of infrared action spectroscopy of peptides in the far-IR and THz regime. An introduction to the field of far-IR spectroscopy is given, thereby highlighting the opportunities that are provided for gas-phase research on neutral peptides. Current experimental methods, including spectroscopic schemes, have been reviewed. Structural information from the experimental far-IR spectra can be obtained with the help of suitable theoretical approaches such as dynamical DFT techniques and the recently developed Graph Theory. The aim of this review is to underline how the synergy between far-IR spectroscopy and theory can provide an unprecedented picture of the structure of neutral biomolecules in the gas phase. The far-IR signatures of the discussed studies are summarized in a far-IR map, in order to gain insight into the origin of the far-IR localized and delocalized motions present in peptides and where they can be found in the electromagnetic spectrum.

Conformational assignment of gas phase peptides and their H-bonded complexes using far-IR/THz: IR-UV ion dip experiment, DFT-MD spectroscopy, and graph theory for mode assignment

Graph theory based vibrational modes as new entities for vibrational THz spectroscopy.

Deconstructing Prominent Bands in the Terahertz Spectra of H7O3+ and H9O4+: Intermolecular Modes in Eigen Clusters

We report experimental vibrational action spectra (210-2200 cm^-1) and calculated IR spectra, using recent ab initio potential energy and dipole moment surfaces, of H₇O₃⁺ and H₉O₄⁺. We focus on prominent far-IR bands, which postharmonic analyses show, arise from two types of intermolecular motions, i.e., hydrogen bond stretching and terminal water wagging modes, that are similar in both clusters. The good agreement between experiment and theory further validates the accuracy of the potential and dipole moment surfaces, which was used in a recent theoretical study that concluded that infrared photodissociation spectra of the cold clusters correspond to the Eigen isomer. The comparison between theory and experiment adds further confirmation of the need of postharmonic analysis for these clusters.

Strategy for Modeling the Infrared Spectra of Ion-Containing Water Drops

Hydrated ions are ubiquitous in environmental and biological media. Understanding the perturbation exerted by an ion on the water hydrogen bond network is possible in the nanodrop regime by recording vibrational spectra in the O-H bond stretching region. This has been achieved experimentally in recent years by forming gaseous ions containing tens to hundreds of water molecules and recording their infrared photodissociation spectra. In this paper, we demonstrate the capabilities of a modeling strategy based on an extension of the AMOEBA polarizable force field to implement water atomic charge fluctuations along with those of intramolecular structure along the dynamics. This supplementary flexibility of nonbonded interactions improves the description of the hydrogen bond network and, therefore, the spectroscopic response. Finite temperature IR spectra are obtained from molecular dynamics simulations by computing the Fourier transform of the dipole moment autocorrelation function. Simulations of 1-2 ns are required for extensive sampling in order to reproduce the experimental spectra. Furthermore, bands are assigned with the driven molecular dynamics approach. This method package is shown to compare successfully with experimental spectra for 11 ions in water drops containing 36-100 water molecules. In particular, band frequency shifts of the free O-H stretching modes at the cluster surface are well reproduced as a function of both ion charge and drop size.

Ab Initio Studies of Structural and Vibrational Properties of Protonated Water Cluster H7O3(+) and Its Deuterium Isotopologues: An Application of Driven Molecular Dynamics

In this work, we present infrared (IR) spectra of H7O3(+) and its deuterium isotopomers calculated by direct molecular dynamics (MD) simulations at the B3LYP/6-31+G** computational level. The calculated spectra obtained at 100, 300, and 500 K were compared to available experimental observations, and spectral features were assigned using normal-mode analysis (NMA) and driven molecular dynamics (DMD). Spectral peaks at 2410 and 2540 cm(-1) were assigned to asymmetric and symmetric stretches of the bridging hydrogen (BH) using NMA. The weak spectral features at 2166 and 2275 cm(-1) were assigned to a combination band of BH asymmetric stretch, H2O in phase wagging, OO stretch, and H3O(+) rocking vibrations by DMD simulations. Our observation of BH stretch vibrations as low as 2166 cm(-1) is in good agreement with the assignment of the low-resolution spectrum obtained by Schwarz at 2200-2300 cm(-1) [Schwarz, H. A. J. Chem. Phys. 1977, 67, 5525-5534] and vibrational predissociation spectrum by Lee et al. ∼2300 cm(-1) [Okumura, M.; Yeh, L. I.; Myers, J. D.; Lee, Y. T. J. Chem. Phys. 1990, 94, 3416-3427].

Dynamic Allostery of the Catabolite Activator Protein Revealed by Interatomic Forces

The Catabolite Activator Protein (CAP) is a showcase example for entropic allostery. For full activation and DNA binding, the homodimeric protein requires the binding of two cyclic AMP (cAMP) molecules in an anti-cooperative manner, the source of which appears to be largely of entropic nature according to previous experimental studies. We here study at atomic detail the allosteric regulation of CAP with Molecular dynamics (MD) simulations. We recover the experimentally observed entropic penalty for the second cAMP binding event with our recently developed force covariance entropy estimator and reveal allosteric communication pathways with Force Distribution Analyses (FDA). Our observations show that CAP binding results in characteristic changes in the interaction pathways connecting the two cAMP allosteric binding sites with each other, as well as with the DNA binding domains. We identified crucial relays in the mostly symmetric allosteric activation network, and suggest point mutants to test this mechanism. Our study suggests inter-residue forces, as opposed to coordinates, as a highly sensitive measure for structural adaptations that, even though minute, can very effectively propagate allosteric signals.

The Catabolite Activator Protein (CAP) is a well-studied example for how cellular catabolite levels are integrated into the gene regulation. Its affinity for a specific stretch of DNA can be switched on by the binding of two nucleotide molecules termed cAMP to its two protomers. Even though the nucleotides occupy structurally identical binding pockets, the second cAMP binding occurs at an affinity orders of magnitude lower than the first cAMP binding. The question arises how, in the absence of structural changes, the first binding can affect the second. An answer from experiments has been that the communication is largely of entropic nature, i.e. the second cAMP binding would lead to a pronounced reduction in atomic fluctuations of the protein without affecting the atomic mean positions. We here revisited this question by performing Molecular Dynamics simulations. By measuring correlations of forces, a newly derived method outperforming the more common coordinate-based approach, we could recover the previously determined entropic penalty. In addition, however, we observed unobtrusive structural changes of side-chain interactions leading to the occlusion of the second binding pocket that add a critical ‘enthalpic’ component hitherto overlooked. Our study provides a mechanistic view onto the intriguing anti-cooperativity of CAP.

Driven molecular dynamics studies of the shared proton motion in the H5O2+·Ar cluster: the effect of argon tagging and deuteration on vibrational spectra

We report IR spectra of H5O2(+) and H5O2(+)·Ar and their deuterium isotopologues using ab initio molecular dynamics. The trajectories were propagated as microcanonical (NVE) ensembles at energies corresponding to temperatures 50 and 100 K. The potential energy surface is calculated on-the-fly at the MP2/aug-cc-pVDZ level of theory. The calculations show that adding an argon atom to H5O2(+) introduces symmetry breaking in the Zundel core ion, causes blueshift in the shared proton vibration by about 200 cm(-1), and leads to the splitting of the OH stretch vibrations into four bands. Driven molecular dynamics (DMD) method is used to assign the spectrum by coupling the dipole moment to an external electric field oscillating at frequency ω. The broad feature at 1100 cm(-1) in the H5O2(+)·Ar spectrum is ascribed to the large amplitude shared proton vibration coupled with torsion and wag modes. MD MP2 simulations predict the H/D redshift in the shared proton vibration and water bending vibration to be about 280 and 460 cm(-1), respectively, in good agreement with experimental observations.

Elucidating structural basis of benzofuran pyrrolidine pyrazole derivatives for enhancing potency against both the InhA enzyme and intact M. tuberculosis cells: a combined MD simulations and 3D-QSAR study

The structural concept for enhancing both IC₅₀ and MIC₉₀ activities summarized from MD simulations and CoMSIA results.

Vibrational mode assignment of finite temperature infrared spectra using the AMOEBA polarizable force field

The Driven Molecular Dynamics approach has been adapted and associated with the AMOEBA polarizable force field to assign and visualize vibrational modes in infrared spectra obtained by molecular dynamics simulations.

Elucidating the structural basis of diphenyl ether derivatives as highly potent enoyl-ACP reductase inhibitors through molecular dynamics simulations and 3D-QSAR study

Diphenyl ether derivatives are good candidates for anti-tuberculosis agents that display a promising potency for inhibition of InhA, an essential enoyl-acyl carrier protein (ACP) reductase involved in fatty acid biosynthesis pathways in Mycobacterium tuberculosis. In this work, key structural features for the inhibition were identified by 3D-QSAR CoMSIA models, constructed based on available experimental binding properties of diphenyl ether inhibitors, and a set of four representative compounds was subjected to MD simulations of inhibitor-InhA complexes for the calculation of binding free energies. The results show that bulky groups are required for the R1 substituent on the phenyl A ring of the inhibitors to favor a hydrophobic pocket formed by residues Phe149, Met155, Pro156, Ala157, Tyr158, Pro193, Met199, Val203, Leu207, Ile215, and Leu218. Small substituents with a hydrophilic property are required at the R3 and R4 positions of the inhibitor phenyl B rings to form hydrogen bonds with the backbones of Gly96 and Met98, respectively. For the R2 substituent, small substituents with simultaneous hydrophilic or hydrophobic properties are required to favor the interaction with the pyrophosphate moiety of NAD(+) and the methyl side chain of Ala198, respectively. The reported data provide structural guidance for the design of new and potent diphenyl ether-based inhibitors with high inhibitory activities against M. tuberculosis InhA.

Rational design of InhA inhibitors in the class of diphenyl ether derivatives as potential anti-tubercular agents using molecular dynamics simulations

A series of diphenyl ether derivatives were developed and showed promising potency for inhibiting InhA, an essential enoyl acyl carrier protein reductase involved in mycolic acid biosynthesis, leading to the lysis of Mycobacterium tuberculosis. To understand the structural basis of diphenyl ether derivatives for designing more potent inhibitors, molecular dynamics (MD) simulations were performed. Based on the obtained results, the dynamic behaviour in terms of flexibility, binding free energy, binding energy decomposition, conformation, and the inhibitor-enzyme interaction of diphenyl ether inhibitors were elucidated. Phe149, Tyr158, Met161, Met199, Val203 and NAD+ are the key residues for binding of diphenyl ether inhibitors in the InhA binding pocket. Our results could provide the structural concept to design new diphenyl ether inhibitors with better enzyme inhibitory activity against M. tuberculosis InhA. The present work facilitates the design of new and potentially more effective anti-tuberculosis agents.

Computational design of thermostabilizing D-amino acid substitutions

Judicious incorporation of D-amino acids in engineered proteins confers many advantages such as preventing degradation by endogenous proteases and promoting novel structures and functions not accessible to homochiral polypeptides. Glycine to D-alanine substitutions at the carboxy termini can stabilize α-helices by reducing conformational entropy. Beyond alanine, we propose additional side chain effects on the degree of stabilization conferred by D-amino acid substitutions. A detailed, molecular understanding of backbone and side chain interactions is important for developing rational, broadly applicable strategies in using D-amino acids to increase protein thermostability. Insight from structural bioinformatics combined with computational protein design can successfully guide the selection of stabilizing D-amino acid mutations. Substituting a key glycine in the Trp-cage miniprotein with D-Gln dramatically stabilizes the fold without altering the protein backbone. Stabilities of individual substitutions can be understood in terms of the balance of intramolecular forces both at the α-helix C-terminus and throughout the protein.

Conformational symmetry and vibrational dynamics of polymers

Polymers are an important class of materials, and their conformation dictates their dynamical, thermodynamical, and hydrodynamical behavior. Several spectroscopic and other techniques have been employed to characterize their conformation. However, little use has been made of group-theoretical techniques except in the classification of symmetry species. In the present review, an attempt has been made to correlate normal modes and their dispersion profiles with the conformation of the polymeric systems. This has been attempted in the case of 2-, 3-, 4-fold and α-helical polymers.

Extracting effective normal modes from equilibrium dynamics at finite temperature

A general method for obtaining effective normal modes of a molecular system from molecular dynamics simulations is presented. The method is based on a localization criterion for the Fourier transformed velocity time-correlation functions of the effective modes. For a given choice of the localization function used, the method becomes equivalent to the principal mode analysis (PMA) based on covariance matrix diagonalization. On the other hand, a proper choice of the localization function leads to a novel method with a strong analogy with the usual normal mode analysis of equilibrium structures, where the Hessian system at the minimum energy structure is replaced by the thermal averaged Hessian, although the Hessian itself is never actually calculated. This method does not introduce any extra numerical cost during the simulation and bears the same simplicity as PMA itself. It can thus be readily applied to ab initio molecular dynamics simulations. Three such examples are provided here. First we recover effective normal modes of an isolated formaldehyde molecule computed at 20 K in very good agreement with the results of a normal mode analysis performed at its equilibrium structure. We then illustrate the applicability of the method for liquid phase studies. The effective normal modes of a water molecule in liquid water and of a uracil molecule in aqueous solution can be extracted from ab initio molecular dynamics simulations of these two systems at 300 K.

Pump-probe molecular dynamics as a tool for studying protein motion and long range coupling

A new method for analyzing the dynamics of proteins is developed and tested. The method, pump-probe molecular dynamics, excites selected atoms or residues with a set of oscillating forces, and the transmission of the impulse to other parts of the protein is probed using Fourier transform of the atomic motions. From this analysis, a coupling profile can be determined which quantifies the degree of interaction between pump and probe residues. Various physical properties of the method such as reciprocity and speed of transmission are examined to establish the soundness of the method. The coupling strength can be used to address questions such as the degree of interaction between different residues at the level of dynamics, and identify propagation of influence of one part of the protein on another via "pathways" through the protein. The method is illustrated by analysis of coupling between different secondary structure elements in the allosteric protein calmodulin, and by analysis of pathways of residue-residue interaction in the PDZ domain protein previously elucidated by genomics and mutational studies.

Coarse-grained normal mode analysis in structural biology

The realization that experimentally observed functional motions of proteins can be predicted by coarse-grained normal mode analysis has renewed interest in applications to structural biology. Notable applications include the prediction of biologically relevant motions of proteins and supramolecular structures driven by their structure-encoded collective dynamics; the refinement of low-resolution structures, including those determined by cryo-electron microscopy; and the identification of conserved dynamic patterns and mechanically key regions within protein families. Additionally, hybrid methods that couple atomic simulations with deformations derived from coarse-grained normal mode analysis are able to sample collective motions beyond the range of conventional molecular dynamics simulations. Such applications have provided great insight into the underlying principles linking protein structures to their dynamics and their dynamics to their functions.

Vibrational Analysis of the H5O2+ Infrared Spectrum Using Molecular and Driven Molecular Dynamics

Standard molecular and driven molecular dynamics are used to analyze prominent spectral features in the H5O2+ infrared spectrum. In the driven method, the molecular Hamiltonian is augmented with a time-dependent term, mu x epsilon(0) sin(omegat), where mu is the dipole moment of H5O2+, epsilon0 is the electric field, and omega is the frequency. The magnitude of the electric field determines whether the driving is mild (the harmonic limit) or strong (anharmonic motion and mode coupling). We analyze the spectrum in the wavenumber range from 600 to 1900 cm(-1), where recent experimental measurements are available for H5O2+. On the basis of the simulations, we have assigned the broad feature around 1000 cm(-1) to the proton transfer coupled with the torsion motion. Intense absorption near 1780 cm(-1) is assigned to the H2O monomer bend coupled with proton transfer.

The vibrational predissociation spectra of the H5O2+∙RGn(RG=Ar,Ne) clusters: Correlation of the solvent perturbations in the free OH and shared proton transitions of the Zundel ion

Predissociation spectra of the H(5)O(2) (+)RG(n)(RG = Ar,Ne) cluster ions are reported in energy regions corresponding to both the OH stretching (3350-3850 cm(-1)) and shared proton (850-1950 cm(-1)) vibrations. The two free OH stretching bands displayed by the Ne complex are quite close to the band origins identified earlier in bare H(5)O(2) (+) [L. I. Yeh, M. Okumura, J. D. Myers, J. M. Price, and Y. T. Lee, J. Chem. Phys. 91, 7319 (1989)], indicating that the symmetrical H(5)O(2) (+) "Zundel" ion remains largely intact in H(5)O(2) (+)Ne. The low-energy spectrum of the Ne complex is simpler than that observed previously for H(5)O(2) (+)Ar, and is dominated by two sharp transitions at 928 and 1047 cm(-1), with a weaker feature at 1763 cm(-1). The H(5)O(2) (+)Ar(n),n = 1-5 spectra generally exhibit complex band structures reflecting solvent-induced symmetry breaking of the Zundel core ion. The extent of solvent perturbation is evaluated with electronic structure calculations, which predict that the rare gas atoms should attach to the spectator OH groups of H(5)O(2) (+) rather than to the shared proton. In the asymmetric complexes, the shared proton resides closer to the more heavily solvated water molecule, leading to redshifts in the rare gas atom-solvated OH stretches and to blueshifts in the shared proton vibrations. The experimental spectra are compared with recent full-dimensional vibrational calculations (diffusion Monte Carlo and multimode/vibrational configuration interaction) on H(5)O(2) (+). These results are consistent with assignment of the strong low-energy bands in the H(5)O(2) (+)Ne spectrum to the vibration of the shared proton mostly along the O-O axis, with the 1763 cm(-1) band traced primarily to the out-of-phase, intramolecular bending vibrations of the two water molecules.