Low-temperature Molecular Dynamics Research Articles

Communication: Towards first principles theory of relaxation in supercooled liquids formulated in terms of cooperative motion.

A general theory of the long time, low temperature dynamics of glass-forming fluids remains elusive despite the almost 20 years since the famous pronouncement by the Nobel Laureate P. W. Anderson, "The deepest and most interesting unsolved problem in solid state theory is probably the theory of the nature of glass and the glass transition" [Science 267, 1615 (1995)]. While recent work indicates that Adam-Gibbs theory (AGT) provides a framework for computing the structural relaxation time of supercooled fluids and for analyzing the properties of the cooperatively rearranging dynamical strings observed in low temperature molecular dynamics simulations, the heuristic nature of AGT has impeded general acceptance due to the lack of a first principles derivation [G. Adam and J. H. Gibbs, J. Chem. Phys. 43, 139 (1965)]. This deficiency is rectified here by a statistical mechanical derivation of AGT that uses transition state theory and the assumption that the transition state is composed of elementary excitations of a string-like form. The strings are assumed to form in equilibrium with the mobile particles in the fluid. Hence, transition state theory requires the strings to be in mutual equilibrium and thus to have the size distribution of a self-assembling system, in accord with the simulations and analyses of Douglas and co-workers. The average relaxation rate is computed as a grand canonical ensemble average over all string sizes, and use of the previously determined relation between configurational entropy and the average cluster size in several model equilibrium self-associating systems produces the AGT expression in a manner enabling further extensions and more fundamental tests of the assumptions.

Low-temperature molecular dynamics simulations of horse heart cytochrome c and comparison with inelastic neutron scattering data

Molecular dynamics (MD) simulation combined with inelastic neutron scattering can provide information about the thermal dynamics of proteins, especially the low-frequency vibrational modes responsible for large movement of some parts of protein molecules. We performed several 30-ns MD simulations of cytochrome c (Cyt c) in a water box for temperatures ranging from 110 to 300 K and compared the results with those from experimental inelastic neutron scattering. The low-frequency vibrational modes were obtained via dynamic structure factors, S(Q, ω), obtained both from inelastic neutron scattering experiments and calculated from MD simulations for Cyt c in the same range of temperatures. The well known thermal transition in structural movements of Cyt c is clearly seen in MD simulations; it is, however, confined to unstructured fragments of loops Ω1 and Ω2; movement of structured loop Ω3 and both helical ends of the protein is resistant to thermal disturbance. Calculated and experimental S(Q, ω) plots are in qualitative agreement for low temperatures whereas above 200 K a boson peak vanishes from the calculated plots. This may be a result of loss of crystal structure by the protein-water system compared with the protein crystal.

Molecular dynamics in cytochrome c oxidase Mössbauer spectra deconvolution

In this work low temperature molecular dynamics simulations of cytochrome c oxidase are used to predict an experimentally observable, namely Mössbauer spectra width. Predicted lineshapes are used to model Lorentzian doublets, with which published cytochrome c oxidase Mössbauer spectra were simulated. Molecular dynamics imposed constraints to spectral lineshapes permit to obtain useful information, like the presence of multiple chemical species in the binuclear center of cytochrome c oxidase. Moreover, a benchmark of quality for molecular dynamic simulations can be obtained. Despite the overwhelming importance of dynamics in electron–proton transfer systems, limited work has been devoted to unravel how much realistic are molecular dynamics simulations results. In this work, molecular dynamics based predictions are found to be in good agreement with published experimental spectra, showing that we can confidently rely on actual simulations. Molecular dynamics based deconvolution of Mössbauer spectra will lead to a renewed interest for application of this approach in bioenergetics.

Identification of Determinants of Ligand Binding Affinity and Selectivity in the Prostaglandin D2 Receptor CRTH2

The chemoattractant receptor-homologous molecule expressed on Th2 cells (CRTH2) is a G protein-coupled receptor that mediates the pro-inflammatory effects of prostaglandin D(2) (PGD(2)) generated in allergic inflammation. The CRTH2 receptor shares greatest sequence similarity with chemoattractant receptors compared with prostanoid receptors. To investigate the structural determinants of CRTH2 ligand binding, we performed site-directed mutagenesis of putative mCRTH2 ligand-binding residues, and we evaluated mutant receptor ligand binding and functional properties. Substitution of alanine at each of three residues in the transmembrane (TM) helical domains (His-106, TM III; Lys-209, TM V; and Glu-268, TM VI) and one in extracellular loop II (Arg-178) decreased PGD(2) binding affinity, suggesting that these residues play a role in binding PGD(2). In contrast, the H106A and E268A mutants bound indomethacin, a nonsteroidal anti-inflammatory drug, with an affinity similar to the wild-type receptor. HEK293 cells expressing the H106A, K209A, and E268A mutants displayed reduced inhibition of intracellular cAMP and chemotaxis in response to PGD(2), whereas the H106A and E268A mutants had functional responses to indomethacin similar to the wild-type receptor. Binding of PGE(2) by the E268A mutant was enhanced compared with the wild-type receptor, suggesting that Glu-268 plays a role in determining prostanoid ligand selectivity. Replacement of Tyr-261 with phenylalanine did not affect PGD(2) binding but decreased the binding affinity for indomethacin. These results provided the first details of the ligand binding pocket of an eicosanoid-binding chemoattractant receptor.

Finite-size effect at both high and low temperatures in molecular dynamics calculations ofthe self-diffusion coefficient and viscosity of liquid silica

Molecular dynamics simulations are performed at 6543 and 3310 K to investigate how theviscosity and self-diffusion coefficient scale with system size in the liquid BKS silica system.We find that at high temperature the finite-size effect on shear viscosity is negligible.However, the size effect on the diffusion coefficient still exists and scales linearly with1/N1/3,where N is the total number of particles in the system. At low temperature, thesize effect on the viscosity becomes stronger than that on diffusion, andthe logarithm of the viscosity and the diffusion coefficient scale linearly with1/N. These results are consistent with previous theoretical developments, and demonstrate thatthe finite-size effect should be considered in both high- and low-temperature moleculardynamics simulations of liquid silica.

NMR solution structure of a potent cyclic nonapeptide inhibitor of ICAM-1-mediated leukocyte adhesion produced by homologous amino acid substitution.

We have previously described a disulfide-linked cyclic nonapeptide (inhibitory peptide-01, IP01), with the sequence CLLRMRSIC, which binds to intercellular adhesion molecule-1 (ICAM-1), and blocks binding to its counter-structure, the integrin alphaLbeta2 (leukocyte functional antigen-1, LFA-1) (Sillerud et al., J. Peptide Res. 62, 2003: 97). We now report the optimization of this peptide by means of single homologous amino acid substitutions to yield a new peptide (IP02-K6; CLLRMKSAC) which shows an approximately sixfold improvement in inhibitory activity of multivalent leukocyte binding (inhibition constant for 50% inhibition, IC50 = 90 microm) compared with IP01 (IC50 = 580 microm). This improvement in activity gives IP02-K6 potent in vivo activity in a murine model of ischemia reperfusion injury (Merchant et al., Am. J. Physiol. Heart Circ. 284, 2003: H1260). In order to determine the structural features relevant to ICAM-1-binding, we have determined the structure of IP02-K6 using proton nuclear magnetic resonance (NMR) spectroscopy and restrained molecular modeling. In our previously reported study of solution models of IP01, we observed three interconverting conformations during low-temperature molecular dynamics simulation. In the present study, we find a single conformation of IP02-K6 similar to one of the previously found conformations of IP01 (family C). In particular, an R4-S7 beta-turn is present in similar proportions in both conformation C of IP01 and in IP02-K6; this motif is important in binding to ICAM-1 because this turn enables the IP02-K6 backbone to drape over proline-36 on ICAM-1. The NMR-derived solution model of IP02-K6 was found to dock at the alphaLbeta2-binding site on ICAM-1 with no changes in peptide backbone conformation. This docking model displaced five of the 15 alphaLbeta2 residues at the ICAM-1-binding site and provided a rationale for understanding the quantitative relationship between IP02-K6 structure and biologic activity.

Real-time NMR Kinetic Studies Provide Global and Residue-specific Information on the Non-cooperative Unfolding of the β-Trefoil Protein, Interleukin-1β

The interleukin-1β (IL-1β) structural motif is a β-trefoil super fold created by six two-stranded β-hairpins. Turns are thus particularly important in creating the topology and the arrangement of β-strands in this structural motif. In contrast to the signals observed in optical studies, real-time NMR kinetic investigations of the denaturant-induced unfolding of interleukin-1β provide direct, global, and residue-specific information on the structural nature of the unfolding reaction. Heterogeneity in the individual amino acid residue kinetics reveals a rugged unfolding landscape. The relative kinetic stability of native-like turns supports low temperature molecular dynamics predictions of turn-controlled unfolding.

The pathway to reorientation in ammonium fluoride

Two alternative transition states for ammonium ion rotation in crystalline NH 4F have been characterised using ab initio calculations. In one all the cations rotate in phase, while in the other the rotating ions are isolated from each other by non-rotated ions. These two transition states are characterised by the cation displacement and changes in the bond lengths and bond angles within the ammonium ion. The pathway for passing over the barrier was investigated by low temperature molecular dynamics runs from each transition state. In both cases the order of events in passing from the ground state to the transition state is first cation displacement towards a neighbouring fluoride ion, secondly, lengthening of the NH bond pointing to this ion, and thirdly, rotation about this bond.

Anomalous fluctuations in low-temperature molecular dynamics simulations.

We present results of a standard (constant energy) molecular dynamics simulation of Lennard-Jones microcrystals at low temperature. The kinetic energy fluctuations exhibit an anomalous behavior, being much larger than expected in a microcanonical ensemble; this is due to a dynamics, which is only weakly chaotic. Such a dynamics does not warrant the equivalence of time and ensemble averages, unless one extends the simulation over exceedingly long times. A similar phenomenon has been recently found in the simulation of a canonical ensemble through a Nos\'e-Hoover dynamics. \textcopyright{} 1996 The American Physical Society.

Conformational study of Met‐enkephalin in its zwitterionic form

AbstractAn extensive exploration of Met‐enkephalin in its zwitterionic form has been carried out, in order to characterize the different low‐energy conformational domains accessible to this pentapeptide. The study builds on previous studies carried out in our lab for the neutral molecule, which provided the initial geometries from which the conformational space of the charged molecule could be scanned. The initial conformations were subjected to a series of high‐ and lowtemperature molecular dynamics simulations. Snapshots along each trajectory were taken, minimized, and used as starting points in further MD trajectories until no lower‐energy conformers could be characterized. The CHARMm force field was used throughout the study for this purpose. The same search strategy was used in these studies simulating two different environmental conditions, a distance‐dependent dielectric ϵ = r and a high constant dielectric ϵ = 80. In the low dielectric environment, the formation of the salt bridge dominates the structure. In the high dielectric environment, the screening of the electrostatic interactions results in weaker intramolecular interactions. In both cases, the Gly2–Gly3 β‐turn‐type structures are preferred over the Gly3–Phe4 turns, in marked contrast to what is found for the neutral molecule. The lowest‐energy structures from both environmental conditions were reoptimized in the presence of a cluster of explicit water molecules. Reoptimization of the structures considering explicit water structures did not result in significant conformational changes for the structures characterized with the ε = r or ε = 80 environments.

Theoretical determination of conformational paths in citrate synthase.

Two methods are developed for the theoretical determination of a conformational path between two well-documented forms, a closed form and the open form [Remington et al. (1982) J. Mol. Biol. 158, 111-152] of pig heart citrate synthase, a dimeric enzyme of 2 x 437 residues. The first method uses the minimization of the sum of the potential energies at a set of equidistant points, according to Elber and Karplus [(1987) Chem. Phys. Lett. 139, 375-380]. The initialization of the algorithm is modified to account for large-angle rotations of many groups by performing the interpolations in the space of internal polar coordinates of a set of generalized Jacobi vectors earlier introduced by Durup [(1991) J. Phys. Chem. 95, 1817-1829] and by carefully testing all choices of directions of rotation for determining the initialized midpoint between the known forms. The path includes intermediate points, created by successive splittings of each interval into two equal parts, with a partial energy minimization performed after each splitting. The minimization encounters the well-known local-minima problem, which here is handled by low-temperature molecular dynamics annealing. It is shown that the best ratio of potential energy decrease to rms deviation is achieved by running the dynamics at 50 K, as compared to 100 K and above. The main character of the path obtained is the occurrence of strong to-and-fro variations of some dihedral angles at specific stages along the path. The second method, which we name directed dynamics, uses only low-temperature molecular dynamics simulations by starting trajectories from each of the two known forms with initial velocities directed toward the other one. The procedure is iterated by restarting trajectory pairs after the points of closest approach of the preceding pair. The two half-paths thus built eventually meet after 70 iterations. This method provides a second path with strong similarities, as well as some differences, with respect to the path obtained by the first method.

Characterization of low-energy conformational domains for Met-enkephalin

An extensive exploration of the conformational hypersurface of Met-enkephalin has been carried out, in order to characterize different low-energy conformational domains accessible to this pentapeptide. The search strategy used consisted of two steps. First, systematic nested rotations were performed using the ECEPP potential. Ninety-two low-energy structures were found and minimized using the CHARMm potential. High and low-temperature molecular dynamics trajectories were then computed for the lowest energy structures in an interative fashion until no lower energy conformers could be found. The same search strategy was used in these studies simulating three different environments, a distance-dependent dielectric epsilon = r, and two constant dielectrics epsilon = 10 and epsilon = 80. The lowest energy structure found in a distance-dependent dielectric is a Gly-Gly beta-II'-type turn. All other structures found for epsilon = r within 10 kcal/mol of this lowest energy structure are also bends. In the more polar environments, the density of conformational states is significantly larger compared to the apolar media. Moreover, fewer hydrogen bonds are formed in the more polar environments, which increases the flexibility of the peptide and results in less structured conformers. Comparisons are made with previous calculations and experimental results.

Lattice Defects in Semiconductors: Towards a Finite-Temperature Theory

The progress of Density Functional (DF) calculations of optical, electrical and structural properties of deep defects in semiconductors has been remarkable over the past ten years. Time seems ripe for the extension of such a theory to finite temperature thermodynamical properties (notably free energies of formation and migration). This paper is a step towards that direction. The quasi-harmonic theory of Lattice Dynamics is proposed as a good candidate for such an extension. One important difficulty is the numerical effort required to extract the full dynamical matrix vs. volume from standard DF band calculations for large systems (30-100 atoms/cell), which are on the other hand needed for a thermodynamically consistent theory. To overcome this serious problem we point out some exact relationships existing between differential properties of trajectories and force constants, which may be used to extract accurate dynamical matrices from low-temperature Molecular Dynamics runs. The formalism is presented and successfully tested for a 32-atom argon crystal with and without a vacancy, for which exact results are available. Such a new formalism is proposed as a powerful tool to evaluate dynamical matrices, and thus phonon spectra and thermodynamical properties, of low-symmetry semiconductor systems like defects and amorphous materials, for which a few Carr-Parrinello Molecular Dynamics trajectories are much cheaper than the huge number of self-consistent calculations traditionally required by frozen-phonon or force-constant techniques.

From Electronic Structure To Thermodynamics

ABSTRACTA simple way to extend the remarkable results of Density Functional calculations to finite-temperature properties of materials is the quasi-harmonic theory of Lattice Dynamics. In this framework a thermodynamically consistent theory needs the complete phonon spectrum for a large periodic system (30–100 atoms/cell) at many different volumes, which poses severe practical limitations. In this paper I present the application to a semiconducting system of a method recently proposed by Bachelet and De Lorenzi to overcome these limitations. Based on low-temperature Molecular-Dynamics trajectories (now possible from first principles for semiconducting systems according to the method of Car and Parrinello), the method is shown to provide accurate dynamical matrices for an 8-atom silicon supercell. Such a successful, preliminary test, together with the fact that for larger and/or lower-symmetry systems the computational effort required by the “trajectory approach” is lower than traditional frozen-phonon or force-constant techniques, suggests its use in the determination of dynamical matrices of larger defect or amorphous systems, and thus in the study of their thermodynamics from first principles.

Temperature-dependent molecular dynamics and restrained X-ray refinement simulations of a Z-DNA hexamer

Molecular dynamics simulations of the Z-DNA hexamer 5BrdC-dG-5BrdC-dG-5BrdC-dG ‡ ‡ Abbreviations used: Brd, bromodeoxy; r.m.s., root-mean-square; MD, molecular dynamics. were performed at several temperatures between 100 K and 300 K. Above 250 K, a strong sequence-dependent flexibility in the nucleic acid is observed, with the guanine sugar and the phosphate of GpC sequences much more mobile than the cytosine sugar and phosphate of CpG sequences. At 300 K, the hexamer is in dynamic equilibrium between several Z forms, including the crystallographically determined Z I and Z II forms. The local base-pair geometry, however, is not very variable, except for the roll of the base-pairs. The hexamer molecular dynamics trajectories have been used to test the restrained parameter crystallographic refinement model for nucleic acids. X-ray diffraction intensities corresponding to observed diffraction data were computed. The average structures obtained from the simulations were then refined against the calculated intensities, using a restrained least-squares program developed for nucleic acids in order to analyse the effects of the refinement model on the derived quantities. In general, the temperature dependence of the atomic fluctuations determined directly from the refined Debye-Waller factors is in reasonably good agreement with the results obtained by calculating the atomic fluctuations directly from the Z-DNA molecular dynamics trajectories. The agreement is best for refinement of temperature factors without restraints. At the highest temperature studied (300 K), the effect of the refinement on the most mobile atoms (phosphates) is to significantly reduce the mean-square atomic fluctuations estimated from the refined Debye-Waller factors below the actual values ( < (Δr) 2 > ≈ 0.5 A ̊ 2 ). Analysis of the temperature-dependence of the mean-square atomic fluctuations provides information concerning the conformational potential within which the atoms move. The calculated temperature-dependence and anharmonicity of the Z-DNA helix are compared with the results observed for proteins. The average structures from the simulations were refined against the experimental X-ray intensities. It is found that low-temperature molecular dynamics simulations provide a useful tool for optimizing the refinement of X-ray structures.