Hierarchical structure of the energy landscape of proteins revisited by time series analysis. I. Mimicking protein dynamics in different time scales

Time series models, which are constructed from the projections of the molecular-dynamics (MD) runs on principal components (modes), are used to mimic the dynamics of two proteins: tendamistat and immunity protein of colicin E7 (ImmE7). Four independent MD runs of tendamistat and three independent runs of ImmE7 protein in vacuum are used to investigate the energy landscapes of these proteins. It is found that mean-square displacements of residues along the modes in different time scales can be mimicked by time series models, which are utilized in dividing protein dynamics into different regimes with respect to the dominating motion type. The first two regimes constitute the dominance of intraminimum motions during the first 5 ps and the random walk motion in a hierarchically higher-level energy minimum, which comprise the initial time period of the trajectories up to 20-40 ps for tendamistat and 80-120 ps for ImmE7. These are also the time ranges within which the linear nonstationary time series are completely satisfactory in explaining protein dynamics. Encountering energy barriers enclosing higher-level energy minima constrains the random walk motion of the proteins, and pseudorelaxation processes at different levels of minima are detected in tendamistat, depending on the sampling window size. Correlation (relaxation) times of 30-40 ps and 150-200 ps are detected for two energy envelopes of successive levels for tendamistat, which gives an overall idea about the hierarchical structure of the energy landscape. However, it should be stressed that correlation times of the modes are highly variable with respect to conformational subspaces and sampling window sizes, indicating the absence of an actual relaxation. The random-walk step sizes and the time length of the second regime are used to illuminate an important difference between the dynamics of the two proteins, which cannot be clarified by the investigation of relaxation times alone: ImmE7 has lower-energy barriers enclosing the higher-level energy minimum, preventing the protein to relax and letting it move in a random-walk fashion for a longer period of time.     

Hierarchical structure of the energy landscape of proteins revisited by time series analysis. II. Investigation of explicit solvent effects

Time series analysis tools are employed on the principal modes obtained from the C(alpha) trajectories from two independent molecular-dynamics simulations of alpha-amylase inhibitor (tendamistat). Fluctuations inside an energy minimum (intraminimum motions), transitions between minima (interminimum motions), and relaxations in different hierarchical energy levels are investigated and compared with those encountered in vacuum by using different sampling window sizes and intervals. The low-frequency low-indexed mode relationship, established in vacuum, is also encountered in water, which shows the reliability of the important dynamics information offered by principal components analysis in water. It has been shown that examining a short data collection period (100 ps) may result in a high population of overdamped modes, while some of the low-frequency oscillations (<10 cm(-1)) can be captured in water by using a longer data collection period (1200 ps). Simultaneous analysis of short and long sampling window sizes gives the following picture of the effect of water on protein dynamics. Water makes the protein lose its memory: future conformations are less dependent on previous conformations due to the lowering of energy barriers in hierarchical levels of the energy landscape. In short-time dynamics (<10 ps), damping factors extracted from time series model parameters are lowered. For tendamistat, the friction coefficient in the Langevin equation is found to be around 40-60 cm(-1) for the low-indexed modes, compatible with literature. The fact that water has increased the friction and that on the other hand has lubrication effect at first sight contradicts. However, this comes about because water enhances the transitions between minima and forces the protein to reduce its already inherent inability to maintain oscillations observed in vacuum. Some of the frequencies lower than 10 cm(-1) are found to be overdamped, while those higher than 20 cm(-1) are slightly increased. As for the long-time dynamics in water, it is found that random-walk motion is maintained for approximately 200 ps (about five times of that in vacuum) in the low-indexed modes, showing the lowering of energy barriers between the higher-level minima.     

Slow dynamics in protein fluctuations revealed by time-structure based independent component analysis: the case of domain motions

Protein dynamics on a long time scale was investigated using all-atom molecular dynamics (MD) simulation and time-structure based independent component analysis (tICA). We selected the lysine-, arginine-, ornithine-binding protein (LAO) as a target protein and focused on its domain motions in the open state. A MD simulation of the LAO in explicit water was performed for 600 ns, in which slow and large-amplitude domain motions of the LAO were observed. After extracting domain motions by rigid-body domain analysis, the tICA was applied to the obtained rigid-body trajectory, yielding slow modes of the LAO's domain motions in order of decreasing time scale. The slowest mode detected by the tICA represented not a closure motion described by a largest-amplitude mode determined by the principal component analysis but a twist motion with a time scale of tens of nanoseconds. The slow dynamics of the LAO were well described by only the slowest mode and were characterized by transitions between two basins. The results show that tICA is promising for describing and analyzing slow dynamics of proteins.     

Coupling of the overall molecular motion with the conformational transitions. II. The full rotational problem

The method proposed in the companion paper for analysing the coupling between overall and internal dynamics is applied to the study of the full rotational motion of a molecule with one internal degree of freedom. For systems characterized by a finite set of stable conformers determined by the minima of the intramolecular potential, a simplified time evolution operator of mixed type is derived, with the continuous diffusion equation and the generalized random walk operator representing the overall rotation and the internal dynamics, respectively. The dependence on the conformational state of the rotational diffusion tensor is one source of coupling between these two types of motion. Another source is represented by the recoil rotations acting on each subunit during a conformational transition. Both conformational-dependent rotational diffusion tensors and recoil rotations can be calculated from a model for the friction exerted by the solvent. Some applications of the theory are presented in relation to the butane molecule and the molecules having the structure of biphenyl, with particular emphasis on the calculation of the experimental observables in NMR and dielectric relaxation measurements.     

Application of time series analysis on molecular dynamics simulations of proteins: a study of different conformational spaces by principal component analysis

Time series analysis is applied on the collective coordinates obtained from principal component analysis of independent molecular dynamics simulations of alpha-amylase inhibitor tendamistat and immunity protein of colicin E7 based on the Calpha coordinates history. Even though the principal component directions obtained for each run are considerably different, the dynamics information obtained from these runs are surprisingly similar in terms of time series models and parameters. There are two main differences in the dynamics of the two proteins: the higher density of low frequencies and the larger step sizes for the interminima motions of colicin E7 than those of alpha-amylase inhibitor, which may be attributed to the higher number of residues of colicin E7 and/or the structural differences of the two proteins. The cumulative density function of the low frequencies in each run conforms to the expectations from the normal mode analysis. When different runs of alpha-amylase inhibitor are projected on the same set of eigenvectors, it is found that principal components obtained from a certain conformational region of a protein has a moderate explanation power in other conformational regions and the local minima are similar to a certain extent, while the height of the energy barriers in between the minima significantly change. As a final remark, time series analysis tools are further exploited in this study with the motive of explaining the equilibrium fluctuations of proteins.     

Strong correlations and Fickian water diffusion in narrow carbon nanotubes

The authors have used atomistic molecular dynamics (MD) simulations to study the structure and dynamics of water molecules inside an open ended carbon nanotube placed in a bath of water molecules. The size of the nanotube allows only a single file of water molecules inside the nanotube. The water molecules inside the nanotube show solidlike ordering at room temperature, which they quantify by calculating the pair correlation function. It is shown that even for the longest observation times, the mode of diffusion of the water molecules inside the nanotube is Fickian and not subdiffusive. They also propose a one-dimensional random walk model for the diffusion of the water molecules inside the nanotube. They find good agreement between the mean-square displacements calculated from the random walk model and from MD simulations, thereby confirming that the water molecules undergo normal mode diffusion inside the nanotube. They attribute this behavior to strong positional correlations that cause all the water molecules inside the nanotube to move collectively as a single object. The average residence time of the water molecules inside the nanotube is shown to scale quadratically with the nanotube length.     

Spatially heterogeneous dynamics and the Adam-Gibbs relation in the Dzugutov liquid

We perform molecular dynamics simulations of a one-component glass-forming liquid and use the inherent structure formalism to test the predictions of the Adam-Gibbs (AG) theory and to explore the possible connection between these predictions and spatially heterogeneous dynamics. We calculate the temperature dependence of the average potential energy of the equilibrium liquid and show that it obeys the Rosenfeld-Tarazona T(3/5) law for low temperature T, while the average inherent structure energy is found to be inversely proportional to temperature at low T, consistent with a Gaussian distribution of potential energy minima. We investigate the shape of the basins around the local minima in configuration space via the average basin vibrational frequency and show that the basins become slightly broader upon cooling. We evaluate the configurational entropy S(conf), a measure of the multiplicity of potential energy minima sampled by the system, and test the validity of the AG relation between S(conf) and the bulk dynamics. We quantify the dynamically heterogeneous motion by analyzing the motion of particles that are mobile on short and intermediate time scales relative to the characteristic bulk relaxation time. These mobile particles move in one-dimensional "strings", and these strings form clusters with a well-defined average cluster size. The AG approach predicts that the minimum size of cooperatively rearranging regions (CRR) of molecules is inversely proportional to S(conf), and recently (Phys. Rev. Lett. 2003, 90, 085506) it has been shown that the mobile-particle clusters are consistent with the CRR envisaged by Adam and Gibbs. We test the possibility that the mobile-particle strings, rather than clusters, may describe the CRR of the Adam-Gibbs approach. We find that the strings also follow a nearly inverse relation with S(conf). We further show that the logarithm of the time when the strings and clusters are maximum, which occurs in the late-beta-relaxation regime of the intermediate scattering function, follows a linear relationship with 1/TS(conf), in agreement with the AG prediction for the relationship between the configurational entropy and the characteristic time for the liquid to undergo a transition to a new configuration. Since strings are the basic elements of the clusters, we propose that strings are a more appropriate measure of the minimum size of a CRR that leads to configurational transitions. That the cluster size also has an inverse relationship with S(conf) may be a consequence of the fact that the clusters are composed of strings.     

Are Two‐Station Biased Random Walkers Always Potential Molecular Motors?

The short answer to the title question is no. Despite their tremendous complexity, many nanomachines are simply one‐dimensional systems undergoing a biased, that is, unidirectional, walk on a two‐minima potential energy curve. The initially prepared state, or station, is higher in energy than the final equilibrium state that is reached after overcoming an energy barrier. All chemical reactions comply with this scheme, which does not necessarily imply that a generic chemical reaction is a potential molecular motor. If the barrier is low, the system may walk back and the motion will have a large purely Brownian component. Alternatively, a large distance from the barrier of either of the two stations may introduce a Brownian component. Starting from a general inequality that leverages on the idea that the amount of heat dissipated along the potential energy curve is a good indication of the effectiveness of the biased walk, we provide guidelines for the selection of the features of artificial molecular motors.     

Generalized essential energy space random walks to more effectively accelerate solute sampling in aqueous environment

Molecular dynamics sampling can be enhanced via the promoting of potential energy fluctuations, for instance, based on a Hamiltonian modified with the addition of a potential-energy-dependent biasing term. To overcome the diffusion sampling issue, which reveals the fact that enlargement of event-irrelevant energy fluctuations may abolish sampling efficiency, the essential energy space random walk (EESRW) approach was proposed earlier. To more effectively accelerate the sampling of solute conformations in aqueous environment, in the current work, we generalized the EESRW method to a two-dimension-EESRW (2D-EESRW) strategy. Specifically, the essential internal energy component of a focused region and the essential interaction energy component between the focused region and the environmental region are employed to define the two-dimensional essential energy space. This proposal is motivated by the general observation that in different conformational events, the two essential energy components have distinctive interplays. Model studies on the alanine dipeptide and the aspartate-arginine peptide demonstrate sampling improvement over the original one-dimension-EESRW strategy; with the same biasing level, the present generalization allows more effective acceleration of the sampling of conformational transitions in aqueous solution. The 2D-EESRW generalization is readily extended to higher dimension schemes and employed in more advanced enhanced-sampling schemes, such as the recent orthogonal space random walk method.     

Harmonic analysis of large systems. III. Comparison with molecular dynamics

Atomic motions in bovine pancreatic trypsin inhibitor (BPTI), derived from molecular dynamics, harmonic analysis, and quasiharmonic analysis, are compared when a single protein model, energy parameters, and environment are employed. Molecular dynamics (MD) was carried out for 2 nanoseconds. An average structure was determined from the last nanosecond of the MD simulation, when no major structural changes were observed. This structure was used for several harmonic analysis calculations as well as for a reference structure for the quasiharmonic analysis, for both full basis and reduced basis sets. In contrast to the harmonic analysis results, the quasiharmonic reduced basis calculation using a spherical harmonics reduced basis provided good agreement with the full basis calculation, suggesting that when anharmonic effects are considered, BPTI can behave as a homogeneous object. An extensive analysis of the normal modes from a diverse set of 201 minimized MD simulation frames was performed. On only the sub-picosecond time scale were energy minima revisited after a transition to another state. This analysis shows that the dynamics average structure is not representative of the simulation frames in terms of energy and vibrational frequencies. For this model of BPTI, 42% of the motion (mean-squared fluctuation) can be attributed to harmonic limit behavior. A spectral analysis of the correlation function of deformation for a particular normal mode or quasiharmonic mode can be used to determine the time scales of motions which correspond to harmonic vibration, large-scale drift, or sharp transitions between local substrates. © 1995 John Wiley & Sons, Inc.     

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.     

More on the melting of Lennard-Jones clusters

The melting of 13-atom clusters interacting via Lennard-Jones potentials has been revisited using molecular dynamics coupled to steepest descent quenches. A procedure was devised to account for the fraction of times the global and local minima of the potential energy surface are accessed during a long trajectory. This quantity presents a sigmoid shape. A phenomenological model of melting is given in terms of a correlated walk that maps the short time excursions among the global and local minima in configuration space. Comparison between the simulation results and the theoretical model shows that the melting transition is well described in terms of the temperature changes of the fraction of high energy minima accessed during the cluster trajectory. Cooperativity is clear from the S shape of this quantity, i.e., the access to a local minimum favours the access to other local minima.     

Multidimensional infrared spectroscopy of the N-H bond motions in formamide

The heterodyned two-dimensional (2D) IR spectra and equilibrium dynamics of the N-H stretching motion of DCONHD in deuterated formamide, DCOND(2), were studied with 80 fs pulses at 3 microm. The time evolution of the heterodyned 2D IR spectra, pump-probe spectra, and photon echo peak shift demonstrate that interstate dynamics is occurring by relaxation of the original N-H excitation. The N-H vibrational frequency correlation function can be expressed as a sum of three exponentials with correlation times 0.24 ps, 0.8 ps, and 11 ps. The intermediate component is attributed to motions of the N-Hcdots, three dots, centeredO unit involving only slight angular variations of the N-H bond. The slow component is attributed to the structure breaking and making. The anisotropy decay confirmed that the significant angular N-H bond motion occurs on the 11 ps time scale. The fast component, which is the least well determined, might correspond to the modulation of the H-bond distance without angular motion. The correlation coefficient between the pumped and relaxed state distributions was +0.51, implying that the excited state phase memory is only slightly diminished by the relaxation of the N-H excitation. The relaxed modes are concluded to be local to the driven N-H mode.     

Dynamic light scattering in sols of a poly(N-vinylcaprolactam) and polyacrylamide mixture

The properties of semidilute solutions of a mixture of polyacrylamide and thermoresponsive polyvinylcaprolactam are studied via dynamic light scattering. It is found that the relaxation time (correlation length) distribution contains several modes. The obtained experimental data are interpreted in terms of the theory describing the dynamics of semidilute and concentrated solutions with allowance for viscoelastic interactions. The investigation of mode dispersion (the dependence of relaxation time on the wave vector of scattering) shows that the fastest and slowest modes are of the diffusion type. The fast mode is the mode of collective diffusion; it may be attributed to the fluctuation motion of polymer units inside blobs. The slowest mode corresponds to the diffusion motion of macromolecule clusters. The rate of relaxation of medium modes nonlinearly depends on the squared wave vector, and these modes are apparently viscoelastic.     

Structural and dynamic features of water and amorphous ice

Structural properties and microscopic dynamics of water and amorphous ice have been studied by the molecular dynamics method. It has been found that the distribution function of the tetrahedricity parameter exhibits two ranges, which correspond to local molecular formations with low and high degrees of tetrahedricity. The number of molecular clusters with a high degree of tetrahedricity grows as temperature decreases. It has been shown that the vibrational density of states comprises two vibrational modes. A low-frequency vibrational mode strongly depends on pressure and is almost independent of temperature, while a high-frequency mode is relevant to the pressure-independent heat motion of molecules. The geometric criterion of hydrogen bonds has been used to evaluate their continuous lifetime as depending on temperature for molecules with different coordination values. The average lifetime of a hydrogen bond substantially depends on the coordination of molecules, with the temperature dependence of the coordination obeying the activation dynamics.     

Clever and efficient method for searching optimal geometries of lennard-jones clusters

An unbiased algorithm for determining global minima of Lennard-Jones (LJ) clusters is proposed in the present study. In the algorithm, a global minimum is searched by using two operators: one modifies a cluster configuration by moving atoms to the most stable positions on the surface of a cluster and the other gives a perturbation on a cluster configuration by moving atoms near the center of mass of a cluster. The moved atoms are selected by employing contribution of the atoms to the potential energy of a cluster. It was possible to find new global minima for LJ506, LJ521, LJ536, LJ537, LJ538, and LJ541 together with putative global minima of LJ clusters of 10-561 atoms reported in the literature. This indicates that the present method is clever and efficient for cluster geometry optimization.     

Ferric and ferrous iron in nitroso-myoglobin: computer simulations of stable and metastable States and their infrared spectra.

The binding of NO to iron is involved in the biological function of many heme proteins. Contrary to ligands like CO and O(2), which only bind to ferrous (Fe(II)) iron, NO binds to both ferrous and ferric (Fe(III)) iron. In a particular protein, the natural oxidation state can therefore be expected to be tailored to the required function. Herein, we present an ab initio potential-energy surface for ferric iron interacting with NO. This potential-energy surface exhibits three minima corresponding to eta(1)-NO coordination (the global minimum), eta(1)-ON coordination and eta(2) coordination. This contrasts with the potential-energy surface for Fe(II)-NO, which exhibits only two minima (the eta(2) coordination mode for Fe(II) is a transition state, not a minimum). In addition, the binding energies of NO are substantially larger for Fe(III) than for Fe(II). We have performed molecular dynamics simulations for NO bound to ferric myoglobin (Mb(III)) and compare these with results obtained for Mb(II). Over the duration of our simulations (1.5 ns), all three binding modes are found to be stable at 200 K and transiently stable at 300 K, with eventual transformation to the eta(1)-NO global-minimum conformation. We discuss the implication of these results related to studies of rebinding processes in myoglobin.     

Effective numbers of modes applied to analysis of internal dynamics of weakly bound clusters

The dependence of the volume of the chaotic component in the internal dynamics of triatomic van der Waals clusters on the angular momentum is calculated using the Monte Carlo and molecular dynamics methods. It has been found that this dependence is nonmonotonic and that its functional form varies for different values of the total energy. The effective number of rotational modes was used to clarify why a change in the volume of chaotic component of the phase space happens for certain values of the angular momentum. We conclude that a large fraction of regular trajectories in relation to all trajectories appears only when there is a possibility for the regular motion to perform a rotation different from that for a chaotic motion. When such difference is small, the regular motion disappears. The effective number of rotational modes can be used to estimate the difference in the type of rotation and is a convenient parameter which controls changes in the dynamics of the system.