Generalized correlation for biomolecular dynamics

Correlated motions in biomolecules are often essential for their function, e.g., allosteric signal transduction or mechanical/thermodynamic energy transport. Because correlated motions in biomolecules remain difficult to access experimentally, molecular dynamics (MD) simulations are particular useful for their analysis. The established method to quantify correlations from MD simulations via calculation of the covariance matrix, however, is restricted to linear correlations and therefore misses part of the correlations in the atomic fluctuations. Herein, we propose a general statistical mechanics approach to detect and quantify any correlated motion from MD trajectories. This generalized correlation measure is contrasted with correlations obtained from covariance matrices for the B1 domain of protein G and T4 lysozyme. The new method successfully quantifies correlations and provides a valuable global overview over the functionally relevant collective motions of lysozyme. In particular, correlated motions of helix 1 together with the two main lobes of lysozyme are detected, which are not seen by the conventional covariance matrix. Overall, the established method misses more than 50% of the correlation. This failure is attributed to both, an interfering and unnecessary dependence on mutual orientations of the atomic fluctuations and, to a lesser extent, attributed to nonlinear correlations. Our generalized correlation measure overcomes these problems and, moreover, allows for an improved understanding of the conformational dynamics by separating linear and nonlinear contributions of the correlation.     

Local feature analysis: a statistical theory for reproducible essential dynamics of large macromolecules

Multivariate statistical methods are widely used to extract functional collective motions from macromolecular molecular dynamics (MD) simulations. In principal component analysis (PCA), a covariance matrix of positional fluctuations is diagonalized to obtain orthogonal eigenvectors and corresponding eigenvalues. The first few eigenvectors usually correspond to collective modes that approximate the functional motions in the protein. However, PCA representations are globally coherent by definition and, for a large biomolecular system, do not converge on the time scales accessible to MD. Also, the forced orthogonalization of modes leads to complex dependencies that are not necessarily consistent with the symmetry of biological macromolecules and assemblies. Here, we describe for the first time the application of local feature analysis (LFA) to construct a topographic representation of functional dynamics in terms of local features. The LFA representations are low dimensional, and like PCA provide a reduced basis set for collective motions, but they are sparsely distributed and spatially localized. This yields a more reliable assignment of essential dynamics modes across different MD time windows. Also, the intrinsic dynamics of local domains is more extensively sampled than that of globally coherent PCA modes.     

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.     

Conformational changes and allosteric communications in human serum albumin due to ligand binding

It is well recognized that knowledge of structure alone is not sufficient to understand the fundamental mechanism of biomolecular recognition. Information of dynamics is necessary to describe motions involving relevant conformational states of functional importance. We carried out principal component analysis (PCA) of structural ensemble, derived from 84 crystal structures of human serum albumin (HSA) with different ligands and/or different conditions, to identify the functionally important collective motions, and compared with the motions along the low-frequency modes obtained from normal mode analysis of the elastic network model (ENM) of unliganded HSA. Significant overlap is observed in the collective motions derived from PCA and ENM. PCA and ENM analysis revealed that ligand selects the most favored conformation from accessible equilibrium structures of unliganded HSA. Further, we analyzed dynamic network obtained from molecular dynamics simulations of unliganded HSA and fatty acids- bound HSA. Our results show that fatty acids-bound HSA has more robust community network with several routes to communicate among different parts of the protein. Critical nodes (residues) identified from dynamic network analysis are in good agreement with allosteric residues obtained from sequence-based statistical coupling analysis method. This work underscores the importance of intrinsic structural dynamics of proteins in ligand recognition and can be utilized for the development of novel drugs with optimum activity.     

Protein dynamics from X-ray crystallography: anisotropic, global motion in diffuse scattering patterns

Understanding X-ray crystallographic diffuse scattering is likely to improve our comprehension of equilibrium collective protein dynamics. Here, using molecular dynamics (MD) simulation, a detailed analysis is performed of the origins of diffuse scattering in crystalline Staphylococcal nuclease, for which the complete diffuse scattering pattern has been determined experimentally. The hydrogen-atom contribution and the scattering range over which the scattering can be considered to be a sum of solvent and protein scattering are determined. Two models of correlated protein motion are investigated by calculating the model-derived diffuse scattering and comparing with the scattering calculated directly from MD trajectories. In one model, previously used in diffuse scattering interpretation, the atomic displacement correlations decay isotropically with increasing separation. Model correlation lengths are obtained by refining the model scattering against the simulation-derived scattering pattern, and are found to be significantly different from those correlation lengths derived directly from the MD trajectories. Furthermore, the convergence between the model-derived and MD-derived scattering is poor. The second model, in which the displacement correlations are calculated from the principal components of the MD trajectories, is capable of fully reproducing the MD-derived diffuse scattering if the approximately 50% lowest-frequency modes are included. However, a small number ( approximately 10) of lowest-frequency and largest-amplitude modes dominates the diffuse scattering and thus the correlated protein motions. A detailed analysis of the principal components is performed. In particular, the effective free energy profile associated with each principle mode is analyzed and the eigenfrequency and damping coefficient computed using a model of Brownian dynamics. Those collective modes with effective frequencies below approximately 0.5 THz, including those that determine the diffuse scattering, are overdamped.     

Low-frequency vibrational modes and infrared absorbance of red, blue and green opsin

Vibrational excitations of low-frequency collective modes are essential for functionally important conformational transitions in proteins. We carried out an analysis of the low-frequency modes in the G protein coupled receptors (GPCR) family of cone opsins based on both normal-mode analysis and molecular dynamics (MD) simulations. Power spectra obtained by MD can be compared directly with normal modes. In agreement with existing experimental evidence related to transmembrane proteins, cone opsins have functionally important transitions that correspond to approximately 950 modes and are found below 80 cm⁻¹. This is in contrast to bacteriorhodopsin and rhodopsin, where the important low-frequency transition modes are below 50 cm⁻¹. We find that the density of states (DOS) profile of blue opsin in a solvent (e.g. water) has increased populations in the very lowest frequency modes (<15 cm⁻¹); this is indicative of the increased thermostability of blue opsin. From our work we found that, although light absorption behaves differently in blue, green and red opsins, their low-frequency vibrational motions are similar. The similarities and differences in the domain motions of blue, red and green opsins are discussed for several representative modes. In addition, the influence of the presence of a solvent is reported and compared with vacuum spectra. We thus demonstrate that terahertz spectroscopy of low-frequency modes might be relevant for identifying those vibrational degrees of freedom that correlate to known conformational changes in opsins. An erratum to this article can be found at     

Projection of Monte Carlo and molecular dynamics trajectories onto the normal mode axes: human lysozyme

A method is presented to describe the internal motions of proteins obtained from molecular dynamics or Monte Carlo simulations as motions of normal mode variables. This method calculates normal mode variables by projecting trajectories of these simulations onto the axes of normal modes and expresses the trajectories as a linear combination of normal mode variables. This method is applied to the result of the molecular dynamics and the Monte Carlo simulations of human lysozyme. The motion of the lowest frequency mode extracted from the simulations represents the hinge bending motion very faithfully. Analysis of the obtained motions of the normal mode variables provides an explanation of the anharmonic aspects of protein dynamics as due first to the anharmonicity of the actual potential energy surface near a minimum and second to trans-minimum conformational changes.     

Model-free methods of analyzing domain motions in proteins from simulation: A comparison of normal mode analysis and molecular dynamics simulation of lysozyme

Model-free methods are introduced to determine quantities pertaining to protein domain motions from normal mode analyses and molecular dynamics simulations. For the normal mode analysis, the methods are based on the assumption that in low frequency modes, domain motions can be well approximated by modes of motion external to the domains. To analyze the molecular dynamics trajectory, a principal component analysis tailored specifically to analyze interdomain motions is applied. A method based on the curl of the atomic displacements is described, which yields a sharp discrimination of domains, and which defines a unique interdomain screw-axis. Hinge axes are defined and classified as twist or closure axes depending on their direction. The methods have been tested on lysozyme. A remarkable correspondence was found between the first normal mode axis and the first principal mode axis, with both axes passing within 3 Å of the alpha-carbon atoms of residues 2, 39, and 56 of human lysozyme, and near the interdomain helix. The axes of the first modes are overwhelmingly closure axes. A lesser degree of correspondence is found for the second modes, but in both cases they are more twist axes than closure axes. Both analyses reveal that the interdomain connections allow only these two degrees of freedom, one more than provided by a pure mechanical hinge. Proteins 27:425–437, 1997. © 1997 Wiley-Liss, Inc.     

Correlative motions and memory effects in molecular dynamics simulations of molecules: principal components and rescaled range analysis suggest that the motions of native BPTI are more correlated than those of its mutants

In this work MD simulations of the native bovine pancreatic trypsin inhibitor (BPTI) and 16 mutants were done in vacuum in order to study memory effects in the mutants using principal component analysis (PCA) and the rescaled range analysis (Hurst exponents). Both PCA and the rescaled range analysis support our previous proposition, based on PCA of lysozyme, that the motions of a native protein are more correlated than those of mutants. The methods are compared, the nature and applications of the rule and the role of the long-range correlations in MD time series (i.e. memory) are discussed in the context of collective motions.     

Domain motions in bacteriophage T4 lysozyme: A comparison between molecular dynamics and crystallographic data

A comparison of a series of extended molecular dynamics (MD) simulations of bacteriophage T4 lysozyme in solvent with X-ray data is presented. Essential dynamics analyses were used to derive collective fluctuations from both the simulated trajectories and a distribution of crystallographic conformations. In both cases the main collective fluctuations describe domain motions. The protein consists of an N- and C-terminal domain connected by a long helix. The analysis of the distribution of crystallographic conformations reveals that the N-terminal helix rotates together with either of these two domains. The main domain fluctuation describes a closure mode of the two domains in which the N-terminal helix rotates concertedly with the C-terminal domain, while the domain fluctuation with second largest amplitude corresponds to a twisting mode of the two domains, with the N-terminal helix rotating concertedly with the N-terminal domain. For the closure mode, the difference in hinge-bending angle between the most open and most closed X-ray structure along this mode is 49 degrees. In the MD simulation that shows the largest fluctuation along this mode, a rotation of 45 degrees was observed. Although the twisting mode has much less freedom than the closure mode in the distribution of crystallographic conformations, experimental results suggest that it might be functionally important. Interestingly, the twisting mode is sampled more extensively in all MD simulations than it is in the distribution of X-ray conformations. Proteins 31:116–127, 1998. © 1998 Wiley-Liss, Inc.     

Event detection and sub‐state discovery from biomolecular simulations using higher‐order statistics: Application to enzyme adenylate kinase

Biomolecular simulations at millisecond and longer time‐scales can provide vital insights into functional mechanisms. Because post‐simulation analyses of such large trajectory datasets can be a limiting factor in obtaining biological insights, there is an emerging need to identify key dynamical events and relating these events to the biological function online, that is, as simulations are progressing. Recently, we have introduced a novel computational technique, quasi‐anharmonic analysis (QAA) (Ramanathan et al., PLoS One 2011;6:e15827), for partitioning the conformational landscape into a hierarchy of functionally relevant sub‐states. The unique capabilities of QAA are enabled by exploiting anharmonicity in the form of fourth‐order statistics for characterizing atomic fluctuations. In this article, we extend QAA for analyzing long time‐scale simulations online. In particular, we present HOST4MD—a higher‐order statistical toolbox for molecular dynamics simulations, which (1) identifies key dynamical events as simulations are in progress, (2) explores potential sub‐states, and (3) identifies conformational transitions that enable the protein to access those sub‐states. We demonstrate HOST4MD on microsecond timescale simulations of the enzyme adenylate kinase in its apo state. HOST4MD identifies several conformational events in these simulations, revealing how the intrinsic coupling between the three subdomains (LID, CORE, and NMP) changes during the simulations. Further, it also identifies an inherent asymmetry in the opening/closing of the two binding sites. We anticipate that HOST4MD will provide a powerful and extensible framework for detecting biophysically relevant conformational coordinates from long time‐scale simulations. Proteins 2012. © 2012 Wiley Periodicals, Inc.     

Bias-free separation of internal and overall motion of biomolecules

Collective internal motions are known to be important for the function of biological macromolecules. It has been discussed in the past whether the application of superimposing algorithms to remove the overall motion from a structural ensemble introduces artificial correlations between distant atoms. Here we present a new method to eliminate residual rotation and translation from cartesian modes derived from a normal mode analysis or from a principal component analysis. Bias-free separation is based on the idea that the addition of modes of pure rotation/translation can compensate the residual overall motion. Removal of overall motion must reduce the "total amount of motion" (TAM) in the mode. Our algorithm allows to back-calculate revised covariance matrices. The approach was applied to two model systems that show residual overall motion, when analyzed using all atoms as reference for the superimposing algorithm. In both cases, our algorithm was capable of eliminating residual covariances caused by the overall motion, while retaining internal covariances even for very distant atoms. A structural ensemble obtained for a 13-ns molecular dynamics simulation of the protein Ribonuclease T1 showed a covariance matrix of the corrected modes with significantly sharper contours after applying the bias-free separation.     

Investigation of the role of disulphide bond in modulating internal motions of BmK AGAP protein by molecular dynamics simulation

Elucidating structural determinants in the functional regions of toxins can provide useful knowledge for designing novel analgesic peptides. A series of 100 ns MD simulations were performed on the scorpion toxin BmK AGAP in native and disulphide bond broken states. The comparison of disulphide bond broken states with the native state showed the α-helix was found to be the key to the analgesic activity. Furthermore, our results revealed disulphide bonds have considerable influence on the functionally important essential modes of motions and the correlations between the motions of the Core domain and the C-terminal region which are involved in the analgesic activity. Therefore, we can conclude that disulphide bonds have a crucial role in modulating the function via adjusting the dynamics of scorpion toxin BmK AGAP molecule.     

Subtle functional collective motions in pancreatic-like ribonucleases: from ribonuclease A to angiogenin

The analysis of the dynamic behavior of enzymes is fundamental to structural biology. A direct relationship between protein flexibility and biological function has been shown for bovine pancreatic ribonuclease (RNase A) (Rasmussen et al., Nature 1992;357:423-424). More recently, crystallographic studies have shown that functional motions in RNase A involve the enzyme beta-sheet regions that move concertedly on substrate binding and release (Vitagliano et al., Proteins 2002;46:97-104). These motions have been shown to correspond to intrinsic dynamic properties of the native enzyme by molecular dynamics (MD) simulations. To unveil the occurrence of these collective motions in other members of pancreatic-like superfamily, we carried out MD simulations on human angiogenin (Ang). Essential dynamics (ED) analyses performed on the trajectories reveal that Ang exhibits collective motions similar to RNase A, despite the limited sequence identity (33%) of the two proteins. Furthermore, we show that these collective motions are also present in ensembles of experimentally determined structures of both Ang and RNase A. Finally, these subtle concerted beta-sheet motions were also observed for other two members of the pancreatic-like superfamily by comparing the ligand-bound and ligand-free structures of these enzymes. Taken together, these findings suggest that pancreatic-like ribonucleases share an evolutionary conserved dynamic behavior consisting of subtle beta-sheet motions, which are essential for substrate binding and release.     

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.     

Local dynamics of DNA probed with optical absorption spectroscopy of bound ethidium bromide.

We have studied the local dynamics of calf thymus double-helical DNA by means of an "optical labeling" technique. The study has been performed by measuring the visible absorption band of the cationic dye ethidium bromide, both free in solution and bound to DNA, in the temperature interval 360-30 K and in two different solvent conditions. The temperature dependence of the absorption line shape has been analyzed within the framework of the vibronic coupling theory, to extract information on the dynamic properties of the system; comparison of the thermal behavior of the absorption band of free and DNA-bound ethidium bromide gave information on the local dynamics of the double helix in the proximity of the chromophore. For the dye free in solution, large spectral heterogeneity and coupling to a "bath" of low-frequency (soft) modes is observed; moreover, anharmonic motions become evident at suitably high temperatures. The average frequency of the soft modes and the amplitude of anharmonic motions depend upon solvent composition. For the DNA-bound dye, at low temperatures, heterogeneity is decreased, the average frequency of the soft modes is increased, and anharmonic motions are hindered. However, a new dynamic regime characterized by a large increase in anharmonic motions is observed at temperatures higher than approximately 280 K. The DNA double helix therefore appears to provide, at low temperatures, a rather rigid environment for the bound chromophore, in which conformational heterogeneity is reduced and low-frequency motions (both harmonic vibrations and anharmonic contributions) are hindered. The system becomes anharmonic at approximately 180 K; however, above approximately 280 K, anharmonicity starts to increase much more rapidly than for the dye free in solution; this can be attributed to the onset of wobbling of the dye in its intercalation site, which is likely connected with the onset of (functionally relevant) DNA motions, involving local opening/unwinding of the double helix. As shown by parallel measurements of the melting curves, these motions precede the melting of the double helix and depend upon solvent composition much more than does the melting itself.     

Identification of functional motions in the adenylate kinase (ADK) protein family by computational hybrid approaches

Based on integrative computational hybrid approaches that combined statistical coupling analysis (SCA), molecular dynamics (MD), and normal mode analysis (NMA), evolutionarily coupled residues involved in functionally relevant motion in the adenylate kinase protein family were identified. The hybrids identified four top-ranking site pairs that belong to a conserved hydrogen bond network that is involved in the enzyme's flexibility. A second group of top-ranking site pairs was identified in critical regions for functional dynamics, such as those related to enzymatic turnover. The high consistency of the results obtained by SCA with NMA (SCA.NMA) and by SCA.MD hybrid analyses suggests that suitable replacement of the matrix of cross-correlation analysis of atomic fluctuations (derived by using NMA) with those based on MD contributes to the identification of such sites by means of a fast computational calculation. The analysis presented here strongly supports the hypothesis that evolutionary forces, such as coevolution at the sequence level, have promoted functional dynamic properties of the adenylate kinase protein family. Finally, these hybrid approaches can be used to identify, at the residue level, protein motion coordination patterns not previously observed, such as in hinge regions.     

Computational studies on the conformations of some large‐ring cyclodextrins (CDn,n = 20, 21, 22, 23)

Computational studies were carried out on the conformations of large‐ring cyclodextrins with degree of polymerization from 20 to 23. Principal component analysis (PCA) was applied for postprocessing of trajectories from conformational search, based on 100.0 ns molecular dynamics simulations. The dominant PCA modes for concerted motions of the macroring atoms were monitored in a lower‐dimensions subspace. The first six lowest indexed principal components contribute more than 90% of the total atomic motions in all cases, with about 70% (CD21) to 83% (CD22) contribution coming from the three highest‐eigenvalue principal components. Representative average geometries of the cyclodextrin macrorings were also obtained for the whole simulation and for the ten 10.0 ns time intervals of the simulation. We concluded that resemblance exists of the representative conformations of these four cyclodextrins with the circularized three‐turn single helical structure proposed for CD21 from small‐angle X‐ray scattering, as well as with the representative conformations of CD26. Chirality, 2011. © 2011 Wiley‐Liss, Inc.     

Role of disulfide bonds in modulating internal motions of proteins to tune their function: molecular dynamics simulation of scorpion toxin Lqh III

A series of 1-ns MD simulations were performed on the scorpion toxin Lqh III in native and disulfide bond broken states. The removal of disulfide bonds has caused hydrogen bond network alteration in the five-residue turn, the long loop, the alpha-helix, the loop connecting strands II and III, and the C-terminal region. In addition and more importantly, it has influenced the amplitude of the fluctuations of five-residue turn, loops, and C-terminal region with a minor effect on the fluctuations of the cysteines in the broken bond sites. These findings suggest that disulfide bonds are not the most important factors in rigidifying their own locations, while they have more important effects at a global scale. Furthermore, our results reveal that disulfide bonds have considerable influence on the functionally important essential modes of motions and the correlations between the motions of the binding site residues. Therefore, we can conclude that disulfide bonds have a crucial role in modulating the function via adjusting the dynamics of scorpion toxin molecules. Although this conclusion cannot be generalized to all peptides and proteins, it demonstrates the importance of more investigations on this aspect of disulfide bond efficacy.