Application of torsion angle molecular dynamics for efficient sampling of protein conformations

We investigate the application of torsion angle molecular dynamics (TAMD) to augment conformational sampling of peptides and proteins. Interesting conformational changes in proteins mainly involve torsional degrees of freedom. Carrying out molecular dynamics in torsion space does not only explicitly sample the most relevant degrees of freedom, but also allows larger integration time steps with elimination of the bond and angle degrees of freedom. However, the covalent geometry needs to be fixed during internal coordinate dynamics, which can introduce severe distortions to the underlying potential surface in the extensively parameterized modern Cartesian-based protein force fields. A "projection" approach (Katritch et al. J Comput Chem 2003, 24, 254-265) is extended to construct an accurate internal coordinate force field (ICFF) from a source Cartesian force field. Torsion crossterm corrections constructed from local molecular fragments, together with softened van der Waals and electrostatic interactions, are used to recover the potential surface and incorporate implicit bond and angle flexibility. MD simulations of dipeptide models demonstrate that full flexibility in both the backbone phi/psi and side chain chi1 angles are virtually restored. The efficacy of TAMD in enhancing conformational sampling is then further examined by folding simulations of small peptides and refinement experiments of protein NMR structures. The results show that an increase of several fold in conformational sampling efficiency can be reliably achieved. The current study also reveals some complicated intrinsic properties of internal coordinate dynamics, beyond energy conservation, that can limit the maximum size of the integration time step and thus the achievable gain in sampling efficiency.     

Comparison between self-guided Langevin dynamics and molecular dynamics simulations for structure refinement of protein loop conformations

This article presents a comparative analysis of two replica‐exchange simulation methods for the structure refinement of protein loop conformations, starting from low‐resolution predictions. The methods are self‐guided Langevin dynamics (SGLD) and molecular dynamics (MD) with a Nosé–Hoover thermostat. We investigated a small dataset of 8‐ and 12‐residue loops, with the shorter loops placed initially from a coarse‐grained lattice model and the longer loops from an enumeration assembly method (the Loopy program). The CHARMM22 + CMAP force field with a generalized Born implicit solvent model (molecular‐surface parameterized GBSW2) was used to explore conformational space. We also assessed two empirical scoring methods to detect nativelike conformations from decoys: the all‐atom distance‐scaled ideal‐gas reference state (DFIRE‐AA) statistical potential and the Rosetta energy function. Among the eight‐residue loop targets, SGLD out performed MD in all cases, with a median of 0.48 Å reduction in global root‐mean‐square deviation (RMSD) of the loop backbone coordinates from the native structure. Among the more challenging 12‐residue loop targets, SGLD improved the prediction accuracy over MD by a median of 1.31 Å, representing a substantial improvement. The overall median RMSD for SGLD simulations of 12‐residue loops was 0.91 Å, yielding refinement of a median 2.70 Å from initial loop placement. Results from DFIRE‐AA and the Rosetta model applied to rescoring conformations failed to improve the overall detection calculated from the CHARMM force field. We illustrate the advantage of SGLD over the MD simulation model by presenting potential‐energy landscapes for several loop predictions. Our results demonstrate that SGLD significantly outperforms traditional MD in the generation and populating of nativelike loop conformations and that the CHARMM force field performs comparably to other empirical force fields in identifying these conformations from the resulting ensembles. Published 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Predicting hydration free energies using all-atom molecular dynamics simulations and multiple starting conformations

Molecular dynamics simulations in explicit solvent were applied to predict the hydration free energies for 23 small organic molecules in blind SAMPL2 test. We found good agreement with experimental results, with an RMS error of 2.82 kcal/mol over the whole set and 1.86 kcal/mol over all the molecules except several hydroxyl-rich compounds where we find evidence for a systematic error in the force field. We tested two different solvent models, TIP3P and TIP4P-Ew, and obtained very similar hydration free energies for these two models; the RMS difference was 0.64 kcal/mol. We found that preferred conformation of the carboxylic acids in water differs from that in vacuum. Surprisingly, this conformational change is not adequately sampled on simulation timescales, so we apply an umbrella sampling technique to include free energies associated with the conformational change. Overall, the results of this test reveal that the force field parameters for some groups of molecules (such as hydroxyl-rich compounds) still need to be improved, but for most compounds, accuracy was consistent with that seen in our previous tests.     

Conformational analysis of bradykinin by annealed molecular dynamics and comparison to NMR-derived conformations

Conformational analysis of bradykinin (BK), a nonapeptide of the sequence RPPGFSPFR, was accomplished using annealed molecular dynamics (AMD) at 1000 K in BIOGRAF 2.2. One hundred anneal cycles produced 100 conformations over approximately 2000 ps. These conformations were compared to structures derived by nuclear magnetic resonance (NMR) methods for similar shape and energy. Energy minimization of relevant conformations using both BIOGRAF 2.2 and AMBER 3.0a revealed that the AMD-determined conformations are in the same energy range as the NMR-determined structures. Also, the shape of the relevant conformations appeared similar, suggesting that AMD is a good tool for the conformational analysis of small peptide ligands. © 1993 John Wiley & Sons, Inc.     

Error analysis of molecular dynamics and fractal time approximants from a combinatorial perspective

Trotter's theorem forms the theoretical basis of most modern molecular dynamics. In essence this theorem states that a time displacement operator (a Lie operator) constructed by exponentiating a sum of noncommuting operators can be approximated by a product of single operators provided the time interval is "very small." In theory "very small" implies infinitesimally small (at which point the approximate product becomes exact), while in practical analysis a finite time interval is divided into several small subintervals or steps. It follows, therefore, that the larger the number of steps the better the approximation to the exact time displacement operator. The question therefore arises: How many steps are sufficient? For bounded operators, standard theorems are available to provide the answer. In this paper we show that a very simple combinatorial formula can be derived which allows the computation of the global differences (as a function of the number of steps) between the Taylor coefficients of the exact time displacement operator and an approximate one constructed by using a finite number of steps. The formula holds for both bounded and nonbounded operators and shows, quantitatively, what is qualitatively expected-that the error decreases with increasing number of steps. Furthermore, the formula applies irrespective of the complexity of the system, boundary conditions, or the thermodynamic ensemble employed for averaging the initial conditions. The analysis yields explicit expressions for the Taylor coefficients which are then used to compute the errors. In the case of the algorithmically based practical numerical simulations in which fixed, albeit small, steps are repeatedly applied, the rise in the number of steps does not reduce the size of the steps but increases the total time of interest. The combinatorial formula shows that, here, the errors diverge. Furthermore, this work can be used to supplement other efforts such as the use of shadow Hamiltonians where the truncation of the series expansion of the latter will produce errors in the higher order propagator moments.     

Cluster analysis of molecular conformations

We describe a method for locating clusters of geometrically similar conformers in ensembles of chemical conformations. We first calculate the pairwise interconformational distance matrix in either torsional or Cartesian space and then use an agglomerative, single-link clustering method to define a hierarchy of clusterings in the same space. Especially good clusterings are distinguished by high values of the separation ratio: the ratio of the shortest intercluster distance to the characteristic threshold distance defining the clustering. We also discuss other statistics. The method has been embodied in a program called XCluster, which can display the distance matrix, the hierarchy of clusterings, and the clustering statistics in a variety of formats. XCluster can also write out the clustered conformations for subsequent or simultaneous viewing with a molecular visualization program. We demonstrate the sorts of insight that this approach affords with examples obtained from conformational search and molecular dynamics procedures. © 1994 by John Wiley & Sons, Inc.     

Spreading dynamics of chain-like monolayers: a molecular dynamics study

Using large-scale molecular dynamics simulations, we have shown previously that the spreading dynamics of sessile drops on solid surfaces can be described in detail using the molecular-kinetic theory of dynamic wetting. Here we present our first steps in extending this approach to investigate the spreading dynamics of Langmuir-Blodgett monolayers. We make use of a monolayer model originally developed by Karaborni and Toxvaerd, but somewhat simplified to facilitate large-scale simulations. Our preliminary results are in good agreement with recent experimental observations and also support a molecular-kinetic interpretation in which the driving force for spreading is the lateral pressure in the monolayer. Away from equilibrium, initial spreading rates are constant and logarithmically dependent on pressure. However, near equilibrium, spreading is pseudo-diffusive and follows the square root of time. In both regimes the controlling factor is the equilibrium frequency of molecular displacements within the monolayer.     

Polarization effects on peptide conformations at water–membrane interface by molecular dynamics simulation

The electrostatic image method was applied to investigate the conformation of peptides characterized by different hydrophobicities in a water–membrane interface model. The interface was represented by a surface of discontinuity between two media with different dielectric constants, taking into account the difference between the polarizabilities of the aqueous medium and the hydrocarbon one. The method consists of a substitution of the real problem, which involves the charges and the induced polarization at the surface of discontinuity, by a simpler problem formed with charges and their images. The electric field due to the polarization induced at the surface by charge q was calculated using a hypothetical charge q′ (image of q), symmetrically located on the opposite side of the surface. The value of q′ was determined using the appropriate electrostatic boundary conditions at the surface. By means of this procedure, the effect of the interface can be introduced easily in the usual force field. We included this extension in the computational package that we are developing for molecular dynamics simulations (Thor ). The peptides studied included hydrophilic tetraaspartic acid (Asp–Asp–Asp–Asp), tetralysine (Lys–Lys–Lys–Lys), hydrophobic tetrapeptide (His–Phe–Arg–Trp), an amphiphilic fragment of β-endorphin, and the signal sequence of the E. coli λ-receptor. The simulation results are in agreement with known experimental data regarding the behavior of peptides at the water–membrane interface. An analysis of the conformational dynamics of the signal sequence peptide at the interface was performed over the course of a few nanoseconds. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 971–982, 1999     

H-Bond in methanol: a molecular dynamics study

A set of molecular dynamics (MD) simulations of methanol-d₄ at three temperatures in the liquid range (200, 250 and 300K) have been carried out. The equations of motion of 256 molecules interacting through a potential model due to Haughney et al. [J. Phys. Chem., 91 (1987) 4934] were solved using the velocity version of the Verlet algorithm. This rather large number of molecules was required for studying the behaviour of the system at momentum transfers as low as 0.25 Å⁻¹. It was found that the system experiences long period fluctuations, and therefore very long MD runs (of the order of 100 ps) are necessary in order to obtain accurate statistical averages. Computed static properties are in good agreement with those reported by Haughney et al. and the neutron weighted g(r) and the static structure factor compare favourably with available neutron diffraction data. The study of time-dependent properties through centre-of-mass autocorrelation functions (VACF, F_s(Q,t) and F(Q,t)) and their memory functions reveals features unknown in simple liquids and very similar to those found in liquid water. A close agreement between centre-of-mass single-particle autocorrelation functions and the translational part of QENS data is also observed. The dynamic structure factor for the centres of mass show distinctive side peaks in the same region of the (Q,ω) plane where recently collective excitations have been studied using coherent neutron scattering thus establishing the presence of propagating short wavelength modes. fa]Presented at the International Symposium on Hydrogen Bond Physics held at Il Ciocco, Barga, Italy, 11–14 September 1990.     

Probing polar solvation dynamics in proteins: a molecular dynamics simulation analysis

Measurements of time-resolved Stokes shifts on picosecond to nanosecond time scales have been used to probe the polar solvation dynamics of biological systems. Since it is difficult to decompose the measurements into protein and solvent contributions, computer simulations are useful to aid in understanding the details of the molecular behavior. Here we report the analysis of simulations of the electrostatic interactions of the rest of the protein and the solvent with 11 residues of the immunoglobulin binding domain B1 of protein G. It is shown that the polar solvation dynamics are position-dependent and highly heterogeneous. The contributions due to interactions with the protein and with the solvent are determined. The solvent contributions are found to vary from negligible after a few picoseconds to dominant on a scale of hundreds of picoseconds. The origin for the latter is found to involve coupled hydration and protein conformational dynamics. The resulting microscopic picture demonstrates that a wide range of possibilities have to be considered in the interpretation of time-resolved Stokes shift measurements.     

Ligand-to-metal charge-transfer dynamics in a blue copper protein plastocyanin: a molecular dynamics study

Equilibrium and nonequilibrium dynamics of a blue copper protein plastocyanin in an oxidized state are studied by molecular dynamics (MD) simulation. Potential energy functions of the lowest seven electronic states, including ligand-to-metal charge-transfer (LMCT) and copper d --> d excited states, were taken from our previous work (Ando, K. J. Phys. Chem. B 2004, 108, 3940), which employed ab initio molecular orbital and density functional calculations on the active-site model. The equilibrium MD simulations in the ground state indicate that ligand motions coupled to transition from the ground state to the LMCT state are mostly represented by stretching and bending vibrations of the Cu-S(Cys) distance, Ndelta(His)-Cu-Ndelta(His) angle, and S(Cys)-Cu-[Ndelta(His)]2 trigonal pyramid structure. The nonequilibrium dynamics on the LMCT potential exhibit rapid decays in which surface crossings to the d --> d and the first excited states occur in 70-80 fs. The crossing dynamics mostly correlate with cleavage of the Cu-S(Cys) bond and the associated response in the Ndelta(His)-Cu-Ndelta(His) moiety. The average dynamics of the vertical energy gap coordinates exhibit an overdamped decay with a recurrence oscillation in 500 fs, which shows clear coherence surviving after the ensemble averaging. This oscillation stems mostly from the recoiling motion of the Ndelta(His)-Cu-Ndelta(His) part. The dynamics of the energy gaps after this coherent oscillation are randomized such that the ensemble average yields flat profiles along time, although each single trajectory exhibits fluctuations with amplitudes large enough to reach surface crossings. These indicate that the relaxation from the LMCT state first occurs via ballistic and coherent potential crossings in 70-80 and 500 fs, followed by thermally activated random transitions.     

NMR studies of molecular conformations in alpha-cyclodextrin

A new approach for analysis of NMR parameters is proposed. The experimental data set includes scalar couplings, NOEs, and residual dipolar couplings. The method, which aims at construction of the conformational distribution function, is applied to alpha-cyclodextrin in isotropic solution and dissolved in a dilute liquid crystal. An attempt to analyze the experimental data using an average molecular conformation resulted in unacceptable errors. Our approach rests on the maximum entropy method (ME), which gives the flattest possible distribution, consistent with the experimental data. Very good agreement between experimental and calculated NMR parameters was observed. In fact, two conformational states were required in order to obtain a satisfactory agreement between calculated and experimental data. In addition, good agreement with Langevin dynamics computer simulations was obtained.     

Dipole moments and molecular conformations of acetylpyrimidines

By comparing experimentally determined dipole moments with those derived from vector-addition and semiempirical computation using MNDO and AM 1 approximations the conformations of 4- and 5-acetylpyrimidines were found to be close to planar, while that of 2-acetylpyrimidine exhibited a torsion angle of 90°.Translated from Khimiya Geterotsiklicheskikh Soedinenii, No. 10, pp. 1387–1390, October, 1992.     

Calculation of alcohol conformations by molecular mechanics

The geometries and energies of simple alcohols were calculated with a molecular mechanics force field. The force field requires the application of the charge interaction model with charges calculated by the CNDO/2 method, the importance of electrostatic interactions for the equilibrium of rotamers about the C-O bond exceeds that of van der Waals interactions. The calculated rotamer populations are discussed with regard to the value of ¹H NMR coupling constants ³JHCOH and other experimental data.     

Predicting conformations and orientations of guests within a water soluble host: a molecular docking approach

In this study we have examined conformations and orientations of guests within a water-soluble host known by the trivial name Octa Acid (OA). Docking program Vina, which was originally developed for screening drug-like molecules, has been used to identify binding modes and affinities of select guest molecules with OA. Hydrophobic guests were encapsulated into the nonpolar cavity of OA capsule owing to solvophobic interactions. Amphiphilic guests were bound by keeping the nonpolar part within the cavity of OA, while pointing the polar anionic group out of the cavity. All these results obtained from the docking study were in accord with experimental findings. The post-complexation attributes of the guests were regulated by available free space and the specific interactions between guest–OA pair, which led to unusual conformations and orientations. This study showed that scoring function available with Vina, which was derived from protein–ligand data set, could successfully predict post-complexed structural features of guests within OA, thus opening opportunities to modulate physical and chemical behavior of guest molecules.     

Thermophoresis in liquids: a molecular dynamics simulation study

Thermophoresis in liquids is studied by molecular dynamics simulation (MD). A theory is developed that divides the problem in the way consistent with the characteristic scales. MD is then conducted to obtain the solution of each problem, which is to be all combined for macroscopic predictions. It is shown that when the temperature gradient is applied to the nonconducting liquid bath that contains neutral particles, there occurs a pressure gradient tangential to the particle surface at the particle-liquid interface. This may induce the flow in the interfacial region and eventually the particle to move. This applies to the material system that interacts through van der Waals forces and may be a general source of the thermophoresis phenomenon in liquids. The particle velocity is linearly proportional to the temperature gradient. And, in a large part of the given temperature range, the particle motion is in the direction toward the cold end and decreases with respect to the temperature. It is also shown that the particle velocity decreases or even reverses its sign in the lowest limit of the temperature range or with a particle of relatively weak molecular interactions with the liquid. The characteristics of the phenomenon are analyzed in molecular details.