A kinetic model for the internal motions of proteins: diffusion between multiple harmonic wells.

The dynamics of collective protein motions derived from Molecular Dynamics simulations have been studied for two small model proteins: initiation factor I and the B1 domain of Protein G. First, we compared the structural fluctuations, obtained by local harmonic approximations in different energy minima, with the ones revealed by large scale molecular dynamics (MD) simulations. It was found that a limited set of harmonic wells can be used to approximate the configurational fluctuations of these proteins, although any single harmonic approximation cannot properly describe their dynamics. Subsequently, the kinetics of the main (essential) collective protein motions were characterized. A dual-diffusion behavior was observed in which a fast type of diffusion switches to a much slower type in a typical time of about 1-3 ps. From these results, the large backbone conformational fluctuations of a protein may be considered as "hopping" between multiple harmonic wells on a basically flat free energy surface.     

Locally accessible conformations of proteins: multiple molecular dynamics simulations of crambin.

Multiple molecular dynamics (MD) simulations of crambin with different initial atomic velocities are used to sample conformations in the vicinity of the native structure. Individual trajectories of length up to 5 ns sample only a fraction of the conformational distribution generated by ten independent 120 ps trajectories at 300 K. The backbone atom conformational space distribution is analyzed using principal components analysis (PCA). Four different major conformational regions are found. In general, a trajectory samples only one region and few transitions between the regions are observed. Consequently, the averages of structural and dynamic properties over the ten trajectories differ significantly from those obtained from individual trajectories. The nature of the conformational sampling has important consequences for the utilization of MD simulations for a wide range of problems, such as comparisons with X-ray or NMR data. The overall average structure is significantly closer to the X-ray structure than any of the individual trajectory average structures. The high frequency (less than 10 ps) atomic fluctuations from the ten trajectories tend to be similar, but the lower frequency (100 ps) motions are different. To improve conformational sampling in molecular dynamics simulations of proteins, as in nucleic acids, multiple trajectories with different initial conditions should be used rather than a single long trajectory.     

A molecular dynamics study of the C-terminal fragment of the L7/L12 ribosomal protein

The crystallographic dimer of the C-terminal fragment (CTF) of the L7/L12 ribosomal protein has been subjected to molecular dynamics (MD) simulations. A 90 picosecond (ps) trajectory for the protein dimer, 19 water molecules and two counter ions has been calculated at constant temperature. Effects of intermolecular interactions on the structure and dynamics have been studied. The exact crystallographic symmetry is lost and the atomic fluctuations differ from one monomer to the other. The average MD structure is more stable than the X-ray one, as judged by accessible surface area and energy calculations. Crystal (non-dimeric) interactions have been simulated in another 40 ps trajectory by using harmonic restraints to represent intermolecular hydrogen bonds. The conformational changes with respect ot the X-ray structure are then virtually suppressed.The unrestrained dimer trajectory has been scanned for cooperative motions involving secondary structure elements. The intrinsic collective motions of the monomer are transmitted via intermolecular contacts to the dimer structure.The existence of a stable dimeric form of CTF, resembling the crystallographic one, has been documented. At the cost of fairly small energy expenditure the dimer has considerable conformational flexibility. This flexibility may endow the dimer with some functional potential as an energy transducer.     

Dynamics and relative stabilities of parallel- and antiparallel-stranded DNA duplexes.

The dynamics and stability of four DNA duplexes are studied by means of molecular dynamics simulations. The four molecules studied are combinations of 4, 15 bases long, single-stranded oligomers, F1, F2, F3, and F4. The sequence of these single strand oligomers are chosen such that F1-F2 and F3-F4 form parallel (ps) DNA double helices, whereas F1-F4 and F2-F3 form antiparallel-stranded (aps) DNA double helices. Simulations were done at low (100 K) and room (300 K) temperatures. At low temperatures the dynamics are quasi-harmonic and the analysis of the trajectories gives good estimates of the low frequency vibrational modes and density of states. These are used to estimate the linear (harmonic) contribution of local fluctuations to the configurational entropy of the systems. Estimates of the differences in enthalpy between ps and aps duplexes show that aps double helices are more stable than the corresponding ps duplexes, in agreement with experiments. At higher temperatures, the distribution of the fluctuations around the average structures are multimodal and estimates of the configurational entropy cannot be obtained. The multi-basin, nonlinear character of the dynamics at 300 K is established using a novel method which extracts large amplitude nonlinear motions from the molecular dynamics trajectories. Our analysis shows that both ps DNA exhibit much larger fluctuations than the two aps DNA. The large fluctuations of ps DNA are explained in terms of correlated transitions in the beta, epsilon, and zeta backbone dihedral angles.     

Protein conformational dynamics of crambin in crystal, solution and in the trajectories of molecular dynamics simulations

Atomic displacement parameters — B factors of the eight crambin crystal structures obtained at 0.54–1.5 Å resolution and temperatures of 100–293 K have been analyzed. The comparable contributions to the B factor values are the intramolecular motions which are modeled by the harmonic vibration calculations and derived from the molecular dynamics simulation (MD) as well as rigid body changes in the position of a protein molecule as a whole. In solution for the average NMR structure of crambin the amplitudes of the backbone atomic fluctuations of the most residues of the segments with the regular backbone conformations are close to the amplitudes of the small scale harmonic vibrations. For the same residues the probability of the medium scale fluctuations fixed by the hydrogen exchange method is very low. The restricted conformational mobility of those segments is coupled with the depressed amplitudes of the fluctuation changes of the tertiary structure registered by the residue accessibility changes in an ensemble of NMR structures that forms the average NMR structure of crambin. The amplitudes of temperature fluctuations of backbone atoms and the tertiary structure raise in the segment with the irregular conformations, turn and loops. In the same segments the amplitudes of the calculated harmonic vibrations also increase, but to a lesser extent and especially in the interhelical loop with the most strong and complicated fluctuation changes of the backbone conformation. In solution for the NMR structure in this loop the conformational transitions occur between the conformational substates separated by the energy barriers, but they are not observed even in the long 100 ns trajectories from the MD simulation of crambin. These strong local fluctuation changes of the structure may play a key role in the protein functioning and modern performance improvements in the MD simulation techniques are oriented to increase the probability of protein appearance in the trajectories from the MD simulations.     

Molecular dynamics simulation of Lewis blood groups and related oligosaccharides.

Molecular dynamics simulations without explicit inclusion of solvent molecules have been performed to study the motions of Lewisa and Lewisb blood group oligosaccharides, and two blood group A tetrasaccharides having type I and type II core chains. The blood group H trisaccharide has also been studied and compared with the blood group A type II core chain. The potential energy surface developed by Rasmussen and co-workers was used with the molecular mechanics code CHARMM. The lowest energy minima of the component disaccharide fragments were obtained from conformational energy mapping. The lowest energy minima of these disaccharide fragments were used to build the tri- and tetrasaccharides that were further minimized before the actual heating/equilibration and dynamics simulations. The trajectories of the disaccharide fragments, e.g., Fuc alpha- (1----4)GlcNAc, Gal beta-(1----4)GlcNAc, etc., show transitions among various minima. However, the oligosaccharides were found to be dynamically stable and no transitions to other minimum energy conformations were observed in the time series of the glycosidic dihedral angles even during trajectories as long as 300 ps. The stable conformations of the glycosidic linkages in the oligosaccharides are not necessarily the same as the minimum energy conformation of the corresponding isolated disaccharides. The average fluctuations of the glycosidic angles in the oligosaccharides were well within the range of +/- 15 degrees. The results of these trajectory calculations were consistent with the relatively rigid single-conformation models derived for these oligosaccharides from 1H-nmr data.     

Filtering molecular dynamics trajectories to reveal low-frequency collective motions: phospholipase A2

A novel method for analysing molecular dynamics trajectories has been developed, which filters out high frequencies using digital signal processing techniques and facilitates focusing on the low-frequency collective motions of proteins. These motions involve low energy slow motions, which lead to important biological phenomena such as domain closure and allosteric effects in enzymes. The filtering method treats each of the atomic trajectories obtained from the molecular dynamics simulation as a "signal". The trajectories of each of the atoms in the system (or any subset of interest) are Fourier transformed to the frequency domain, a filtering function is applied and then an inverse transformation back to the time domain yields the filtered trajectory. The filtering method has been used to study the dynamics of the enzyme phospholipase A2. In the filtered trajectory, all the high frequency bond and valence angle vibrations were eliminated, leaving only low-frequency motion, mainly fluctuations in torsions and conformational transitions. Analysis of this trajectory revealed interesting motions of the protein, including concerted movements of helices, and changes in shape of the active site cavity. Unlike normal mode analysis, which has been used to study the motion of proteins, this method does not require converged minimizations or diagonalization of a matrix of second derivatives. In addition, anharmonicity, multiple minima and conformational transitions are treated explicitly. Thus, the filtering method avoids most of the approximations implicit in other investigations of the dynamic behaviour of large systems.     

Picosecond timescale rigid-helix and side-chain motions in deoxymyoglobin

The contribution of rigidbody motions to the atomic trajectories in a 100 ps molecular dynamics simulation of deoxymyoglobin is examined. Two typesof rigid-body motions are considered: one in which the helices are rigid units and one in which the side-chains are rigid units. Using a quaternionbased algorithm, fits of the rigid reference structures are made to each time frame of the simulation to derive trajectories of the rigid-body motions. The fitted trajectories are analysed in terms of atomic position fluctuations, mean-square displacements as a function of time, velocity autocorrelation functions and densities of states. The results are compared with the corresponding quantities calculated from the full trajectory. The relative contribution of the rigid helix motions to the helix atom dynamics depends on which quantity is examined and on which subset of atoms is chosen: rigid-helix motions contribute 86% of the rms helix backbone atomic position fluctuations, but 30% of the helix,: atom (backbone and side-chain) mean square displacements and only 1.1% of total kinetic energy. Only very low-frequency motions contribute to the rigid-helix dynamics; the rigid-body analysis allows characteristic rigid-helix vibrations to be identified and described. Treating the side-chains as rigid bodies is foundto be an excellent approximation to both their diffusive and vibrationalmean-square displacements: 96% of side-chain atom mean-square displacements originate from rigid side-Chain motions. However, the errors in theside-chain atomic positional fits are not always small. An analysis is madeof factors contributing to the positional error for different types of side-chain. © Wiley-Liss, Inc.     

Extraction of the energetics of selected types of motion from molecular dynamics trajectories by filtering

A novel method for analyzing molecular dynamics trajectories has been developed which enables the study of selected motions and the corresponding energetics. In particular, it is possible to filter out the high-frequency motions and focus on the structural and energetic features of low-frequency collective motions. The trajectories of the properties of interest are Fourier transformed to the frequency domain, a filtering function is applied, and then an inverse transformation back to the time domain yields the filtered trajectory. The method is demonstrated for harmonic fluctuations and conformational transitions of acetamide and N-acetylalanine N-methylamide, as models for peptides and proteins.     

Conformational study of a nine residue fragment of the antigenic loop of foot-and-mouth disease virus.

The nine-residue peptide Ac-TASARGDLA-NHMe was selected as model peptide in order to understand the conformational features of the antigenic loop of foot-and-mouth disease virus (FMDV). A throughout exploration of the conformational space has been carried out by means of molecular dynamics (MD) and energy minimization. The calculations have been carried out using the AMBER force field. Solvent effects have been included by an effective dielectric constant of epsilon = 4r. The lowest energy conformation presents a secondary structure constituted by an alpha-helix at the N-terminal end followed by two gamma-turns in the central region. The rest of the accessible minima found present also a high tendency to form gamma-turns. Finally, a 100 ps MD trajectory calculation at 298 K suggest a stability of the secondary structure elements of the lowest energy conformation.     

Simulating nanoscale functional motions of biomolecules

We are describing efficient dynamics simulation methods for the characterization of functional motion of biomolecules on the nanometer scale. Multivariate statistical methods are widely used to extract and enhance functional collective motions from molecular dynamics (MD) simulations. A dimension reduction in MD is often realized through a principal component analysis (PCA) or a singular value decomposition (SVD) of the trajectory. Normal mode analysis (NMA) is a related collective coordinate space approach, which involves the decomposition of the motion into vibration modes based on an elastic model. Using the myosin motor protein as an example we describe a hybrid technique termed amplified collective motions (ACM) that enhances sampling of conformational space through a combination of normal modes with atomic level MD. Unfortunately, the forced orthogonalization of modes in collective coordinate space leads to complex dependencies that are not necessarily consistent with the symmetry of biological macromolecules and assemblies. In many biological molecules, such as HIV-1 protease, reflective or rotational symmetries are present that are broken using standard orthogonal basis functions. We present a method to compute the plane of reflective symmetry or the axis of rotational symmetry from the trajectory frames. Moreover, we develop an SVD that best approximates the given trajectory while respecting the symmetry. Finally, we describe a local feature analysis (LFA) to construct a topographic representation of functional dynamics in terms of local features. The LFA representations are low-dimensional, and provide a reduced basis set for collective motions, but unlike global collective modes they are sparsely distributed and spatially localized. This yields a more reliable assignment of essential dynamics modes across different MD time windows.     

Essential dynamics of reversible peptide folding: memory-free conformational dynamics governed by internal hydrogen bonds

A principal component analysis has been applied on equilibrium simulations of a beta-heptapeptide that shows reversible folding in a methanol solution. The analysis shows that the configurational space contains only three dense sub-states. These states of relatively low free energy correspond to the "native" left-handed helix, a partly helical intermediate, and a hairpin-like structure. The collection of unfolded conformations form a relatively diffuse cloud with little substructure. Internal hydrogen-bonding energies were found to correlate well with the degree of folding. The native helical structure folds from the N terminus; the transition from the major folding intermediate to the native helical structure involves the formation of the two most C-terminal backbone hydrogen bonds. A four-state Markov model was found to describe transition frequencies between the conformational states within error limits, indicating that memory-effects are negligible beyond the nanosecond time-scale. The dominant native state fluctuations were found to be very similar to unfolding motions, suggesting that unfolding pathways can be inferred from fluctuations in the native state. The low-dimensional essential subspace, describing 69% of the collective atomic fluctuations, was found to converge at time-scales of the order of one nanosecond at all temperatures investigated, whereas folding/unfolding takes place at significantly longer time-scales, even above the melting temperature.     

Characterization of the conformational behavior of peptide Contryphan Vn: a theoretical study

In this work we report the study of a peptide, the Contryphan Vn produced by Conus ventricosus, a vermivorous cone snail living in the temperate Mediterranean sea. This cyclic peptide of nine residues is a ring closed by a Cys-Cys (Cys: cysteine) disulfide bond containing two proline (Pro) residues and two tryptophans (Trp), one of them being a D-Trp. We present a statistical mechanical characterization of the peptide, simulated in water for about 200 ns with classical molecular dynamics (MD). In recent years there has been a growing interest in the study of the mechanics and dynamics of biological molecules, and in particular for proteins and peptides, about the relationship between collective motions and the active conformations which exert the biological function. To this aim we used the essential dynamics analysis on the MD trajectory and extracted, from the total fluctuations of the molecule, the dominant dynamical modes responsible of the principal conformational transitions. The Contryphan Vn small size allowed us to investigate in details the all-atoms dynamics and the corresponding thermodynamics in conformational space defined by the most significant intramolecular motions.     

Structural basis of hierarchical multiple substates of a protein. II: Monte Carlo simulation of native thermal fluctuations and energy minimization

Conformational fluctuations in a globular protein, bovine pancreatic trypsin inhibitor, in the time range between picoseconds and nanoseconds are studied by a Monte Carlo simulation method. Multiple energy minima are derived from sampled conformations by minimizing their energy. They are distributed in clusters in the conformational space. A hierarchical structure is observed in the simulated dynamics. In the time range between 10(-14) and 10(-10) seconds dynamics is well represented by a superposition of vibrational motions within an energy well with transitions among minima within each cluster. Transitions among clusters take place in the time range of nanoseconds or longer.     

Collective dynamics of EcoRI-DNA complex by elastic network model and molecular dynamics simulations

Anisotropic network model (ANM) is used to analyze the collective motions of restriction enzyme EcoRI in free form and in complex with DNA. For comparison, three independent molecular dynamics (MD) simulations, each of 1.5 ns duration, are also performed for the EcoRI-DNA complex in explicit water. Although high mobility (equilibrium fluctuations) of inner and outer loops that surround the DNA is consistent in both methods and experiments, MD runs sample different conformational subspaces from which reliable collective dynamics cannot be extracted. However, ANM employed on different conformations from MD simulations indicates very similar collective motions. The stems of the inner loops are quite immobile even in the free enzyme and form a large, almost fixed, pocket for DNA binding. As a result, the residues that make specific and non-specific interactions with the DNA exhibit very low fluctuations in the free enzyme. The vibrational entropy difference between the EcoRI complex and free protein + unkinked DNA is positive (favorable), which may partially counteract the unfavorable enthalpy difference of DNA kink formation. Dynamic domains in EcoRI complex and cross-correlations between residue fluctuations indicate possible means of communication between the distal active sites.     

Stability and fluctuations of amide hydrogen bonds in a bacterial cytochrome c: a molecular dynamics study

Molecular dynamics (MD) simulations on a bacterial cytochrome c were performed to investigate the lifetime and fluctuations of backbone hydrogen bonds and to correlate these data with protection factors for hydrogen exchange measured by NMR spectroscopy (Bartalesi et al. in Biochemistry, 42:10923–10930, 2003). The MD simulations provide a consistent pattern in that long lifetimes of hydrogen bonds go along with small amplitude fluctuations. In agreement with experiments, differences in stability were found with a rather flexible N-terminal segment as compared with a more rigid C-terminal part. Protection factors of backbone hydrogen exchange correlate strongly with the number of contacts but also with hydrogen-bond occupancy, hydrogen-bond survival times, as well as the inverse of fluctuations of backbone atoms and hydrogen-bond lengths derived from MD simulation data. We observed a conformational transition in the C-terminal loop, and significant motion in the N-terminal loop, which can be interpreted as being the structural units involved in the onset of the protein unfolding process in agreement with experimental evidence on mitochondrial cytochrome c. Electronic Supplementary Material Supplementary material is available for this article at and is accessible for authorized users. Gernot Kieseritzky and Giulia Morra both contributed equally to this work.     

Temperature dependence of the structure and dynamics of myoglobin. A simulation approach

The results of simulations of the structure and internal motions of carbonomonoxymyoglobin (MbCO) at two different temperatures (325 and 80 K) are presented and compared with experimental data. Properties calculated from the 120 ps trajectory at 325 K are used as a reference in the analysis of the motion of the protein at 80 K. Three separate 80 K molecular dynamics trajectories were calculated; they were started with different coordinate sets from the 325 K simulation and the lower temperature was achieved by scaling the velocities. The simulations yield results for the structural changes between 325 and 80 K that are in general accord with those from X-ray data. Both the experimental and calculated radii of gyration, distances from the center of mass and main-chain difference distance matrices show that there is a significant but inhomogeneous shrinkage with decreasing temperature. For the atomic fluctuations, by contrast, the calculated temperature dependence is very different from the X-ray results; i.e. the calculated root-mean-square backbone fluctuations decrease to 0.11 A at 80 K from 0.51 A at 325 K, while the fluctuations obtained from the X-ray B factors go from 0.56 A at 260 K to 0.47 A at 80 K. The smaller temperature dependence of the B factors suggests that there is significant conformational disorder in MbCO crystals at lower temperatures. This is in accord with the simulation results, which show that the protein is trapped in restricted regions of conformational space at 80 K, while at 325 K a much larger region is accessible to the protein. Analysis of the fluctuations at 325 K and 80 K shows that the room temperature flexibility of the protein is determined by the mobility of the loop regions and by side-chain torsional motions (in accord with earlier simulation results), while the low temperature fluctuations involve motion within a single well. Examination of the calculated iron atom fluctuations and comparison with Mossbauer data show good agreement. It is found that the dominant contribution to the iron motion arises from heme sliding; motion of the iron relative to the heme are much smaller.     

Molecular dynamics analysis of the engrailed homeodomain-DNA recognition

Molecular dynamics (MD) simulations were performed for investigating the role of Gln50 in the engrailed homeodomain-DNA recognition. Employing the crystal structure of free engrailed homeodomain and homeodomain-DNA complex as a starting structure, we carried out MD simulations of: (i) the complex between engrailed homeodomain and a 20 base-pair DNA containing TAATTA core sequence; (ii) the free engrailed homeodomain. The simulations show that homeodomain flexibility does not depend on its ligation state. The engrailed homeodomain shows similar flexibility, and the recognition helix-3 shows very similar characteristic of high rigidity and limited conformational space in two complexation states. At the same time, DNA structure has also no obvious conformational fluctuations. These results preclude the possibility of the side chain of Gln50 forming direct hydrogen bonds to the core DNA bases. MD simulations confirm a few well-conserved sites for water-mediated hydrogen bonds from protein to DNA are occupied by water molecules, and Gln50 interacts with corresponding core DNA bases through water-mediated hydrogen bonds. So Gln50 plays a relatively modest role in determining the affinity and specificity of the engrailed homeodomain. In addition, the electrostatic interaction between homeodomain and phosphate backbone of the DNA is a main factor for N- and C-terminal arm becoming ordered upon DNA binding.