Conformational sampling and nucleotide-dependent transitions of the GroEL subunit probed by unbiased molecular dynamics simulations

GroEL is an ATP dependent molecular chaperone that promotes the folding of a large number of substrate proteins in E. coli. Large-scale conformational transitions occurring during the reaction cycle have been characterized from extensive crystallographic studies. However, the link between the observed conformations and the mechanisms involved in the allosteric response to ATP and the nucleotide-driven reaction cycle are not completely established. Here we describe extensive (in total long) unbiased molecular dynamics (MD) simulations that probe the response of GroEL subunits to ATP binding. We observe nucleotide dependent conformational transitions, and show with multiple 100 ns long simulations that the ligand-induced shift in the conformational populations are intrinsically coded in the structure-dynamics relationship of the protein subunit. Thus, these simulations reveal a stabilization of the equatorial domain upon nucleotide binding and a concomitant "opening" of the subunit, which reaches a conformation close to that observed in the crystal structure of the subunits within the ADP-bound oligomer. Moreover, we identify changes in a set of unique intrasubunit interactions potentially important for the conformational transition.     

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.     

Going beyond Clustering in MD Trajectory Analysis: An Application to Villin Headpiece Folding

Recent advances in computing technology have enabled microsecond long all-atom molecular dynamics (MD) simulations of biological systems. Methods that can distill the salient features of such large trajectories are now urgently needed. Conventional clustering methods used to analyze MD trajectories suffer from various setbacks, namely (i) they are not data driven, (ii) they are unstable to noise and changes in cut-off parameters such as cluster radius and cluster number, and (iii) they do not reduce the dimensionality of the trajectories, and hence are unsuitable for finding collective coordinates. We advocate the application of principal component analysis (PCA) and a non-metric multidimensional scaling (nMDS) method to reduce MD trajectories and overcome the drawbacks of clustering. To illustrate the superiority of nMDS over other methods in reducing data and reproducing salient features, we analyze three complete villin headpiece folding trajectories. Our analysis suggests that the folding process of the villin headpiece is structurally heterogeneous.     

Hybrid approaches to molecular simulation

Molecular dynamics (MD) simulation is an established method for studying the conformational changes that are important for protein function. Recent advances in hardware and software have allowed MD simulations over the same timescales as experiment, improving the agreement between theory and experiment to a large extent. However, running such simulations are costly, in terms of resources, storage, and trajectory analysis. There is still a place for techniques that involve short MD simulations. In order to overcome the sampling paucity of short time-scales, hybrid methods that include some form of MD simulation can exploit certain features of the system of interest, often combining experimental information in surprising ways. Here, we review some recent hybrid approaches to the simulation of proteins.     

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.     

Molecular dynamics simulations of lignin peroxidase in solution

The dynamical and structural properties of lignin peroxidase and its Trp171Ala mutant have been investigated in aqueous solution using molecular dynamics (MD) simulations. In both cases, the enzyme retained its overall backbone structure and all its noncovalent interactions in the course of the MD simulations. Very interestingly, the analysis of the MD trajectories showed the presence of large fluctuations in correspondence of the residues forming the heme access channel; these movements enlarge the opening and facilitate the access of substrates to the enzyme active site. Moreover, steered molecular dynamics docking simulations have shown that lignin peroxidase natural substrate (veratryl alcohol) can easily approach the heme edge through the access channel.     

DROIDS 3.0—Detecting Genetic and Drug Class Variant Impact on Conserved Protein Binding Dynamics

The application of statistical methods to comparatively framed questions about the molecular dynamics (MD) of proteins can potentially enable investigations of biomolecular function beyond the current sequence and structural methods in bioinformatics. However, the chaotic behavior in single MD trajectories requires statistical inference that is derived from large ensembles of simulations representing the comparative functional states of a protein under investigation. Meaningful interpretation of such complex forms of big data poses serious challenges to users of MD. Here, we announce Detecting Relative Outlier Impacts from Molecular Dynamic Simulation (DROIDS) 3.0, a method and software package for comparative protein dynamics that includes maxDemon 1.0, a multimethod machine learning application that trains on large ensemble comparisons of concerted protein motions in opposing functional states generated by DROIDS and deploys learned classifications of these states onto newly generated MD simulations. Local canonical correlations in learning patterns generated from independent, yet identically prepared, MD validation runs are used to identify regions of functionally conserved protein dynamics. The subsequent impacts of genetic and/or drug class variants on conserved dynamics can also be analyzed by deploying the classifiers on variant MD simulations and quantifying how often these altered protein systems display opposing functional states. Here, we present several case studies of complex changes in functional protein dynamics caused by temperature, genetic mutation, and binding interactions with nucleic acids and small molecules. We demonstrate that our machine learning algorithm can properly identify regions of functionally conserved dynamics in ubiquitin and TATA-binding protein (TBP). We quantify the impact of genetic variation in TBP and drug class variation targeting the ATP-binding region of Hsp90 on conserved dynamics. We identify regions of conserved dynamics in Hsp90 that connect the ATP binding pocket to other functional regions. We also demonstrate that dynamic impacts of various Hsp90 inhibitors rank accordingly with how closely they mimic natural ATP binding.     

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.     

Toward a detailed description of the pathways of allosteric communication in the GroEL chaperonin through atomistic simulation

GroEL, along with its coprotein GroES, is essential for ensuring the correct folding of unfolded or newly synthesized proteins in bacteria. GroEL is a complex, allosteric molecule, composed of two heptameric rings stacked back to back, that undergoes large structural changes during its reaction cycle. These structural changes are driven by the cooperative binding and subsequent hydrolysis of ATP, by GroEL. Despite numerous previous studies, the precise mechanisms of allosteric communication and the associated structural changes remain elusive. In this paper, we describe a series of all-atom, unbiased, molecular dynamics simulations over relatively long (50-100 ns) time scales of a single, isolated GroEL subunit and also a heptameric GroEL ring, in the presence and absence of ATP. Combined with results from a distance restraint-biased simulation of the single ring, the atomistic details of the earliest stages of ATP-driven structural changes within this complex molecule are illuminated. Our results are in broad agreement with previous modeling studies of isolated subunits and with a coarse-grained, forcing simulation of the single ring. These are the first reported all-atom simulations of the GroEL single-ring complex and provide a unique insight into the role of charged residues K80, K277, R284, R285, and E388 at the subunit interface in transmission of the allosteric signal. These simulations also demonstrate the feasibility of performing all-atom simulations of very large systems on sufficiently long time scales on typical high performance computing facilities to show the origins of the earliest events in biologically relevant processes.     

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.     

Multiscale modelling of material failure using particles

In this paper, we explore a prototypical model of material fragmentation – the penetration of a plate by a fast moving ball in two dimensions using both smooth particle applied mechanics (SPAM) and molecular dynamics (MD). Using a family of short-ranged pair potentials, we investigate the fracture as a function of impact velocity and plate strength. The MD provides data from which we make a direct comparison with SPAM simulations of the same system. The parameters for the SPAM simulations: equation of state, bulk and shear moduli, yield strength and tensile strength are obtained from the pair potentials and via non-equilibrium MD (tension test). By careful parameterisation of the SPAM model, we are able to observe qualitatively similar behaviour to that shown in the MD simulations of a single crystal (brittle) specimen; bending at low penetrator speed, cracking at higher speeds and eventually, complete penetration. The behaviour depends on the plate strength, which is directly related to details of the force law. For plates characterised as high strength, the agreement between MD and SPAM is evidently poor. We discuss the possible causes for this discrepancy. Nevertheless, SPAM, when used in conjunction with accurate atomistic simulation, is shown to be a highly promising tool for studying failure in materials.     

Flexible protein docking refinement using pose-dependent normal mode analysis

Modeling conformational changes in protein docking calculations is challenging. To make the calculations tractable, most current docking algorithms typically treat proteins as rigid bodies and use soft scoring functions that implicitly accommodate some degree of flexibility. Alternatively, ensembles of structures generated from molecular dynamics (MD) may be cross-docked. However, such combinatorial approaches can produce many thousands or even millions of docking poses, and require fast and sensitive scoring functions to distinguish them. Here, we present a novel approach called "EigenHex," which is based on normal mode analyses (NMAs) of a simple elastic network model of protein flexibility. We initially assume that the proteins to be docked are rigid, and we begin by performing conventional soft docking using the Hex polar Fourier correlation algorithm. We then apply a pose-dependent NMA to each of the top 1000 rigid body docking solutions, and we sample and re-score multiple perturbed docking conformations generated from linear combinations of up to 20 eigenvectors using a multi-threaded particle swarm optimization algorithm. When applied to the 63 "rigid body" targets of the Protein Docking Benchmark version 2.0, our results show that sampling and re-scoring from just one to three eigenvectors gives a modest but consistent improvement for these targets. Thus, pose-dependent NMA avoids the need to sample multiple eigenvectors and it offers a promising alternative to combinatorial cross-docking.     

Effects of surface water on protein dynamics studied by a novel coarse-grained normal mode approach

Normal mode analysis (NMA) has received much attention as a direct approach to extract the collective motions of macromolecules. However, the stringent requirement of computational resources by classical all-atom NMA limits the size of the macromolecules to which the method is normally applied. We implemented a novel coarse-grained normal mode approach based on partitioning the all-atom Hessian matrix into relevant and nonrelevant parts. It is interesting to note that, using classical all-atom NMA results as a reference, we found that this method generates more accurate results than do other coarse-grained approaches, including elastic network model and block normal mode approaches. Moreover, this new method is effective in incorporating the energetic contributions from the nonrelevant atoms, including surface water molecules, into the coarse-grained protein motions. The importance of such improvements is demonstrated by the effect of surface water to shift vibrational modes to higher frequencies and by an increase in overlap of the coarse-grained eigenvector space (the motion directions) with that obtained from molecular dynamics simulations of solvated protein in a water box. These results not only confirm the quality of our method but also point out the importance of incorporating surface structural water in studying protein dynamics.     

Principal component analysis for protein folding dynamics

Protein folding is considered here by studying the dynamics of the folding of the triple β-strand WW domain from the Formin-binding protein 28. Starting from the unfolded state and ending either in the native or nonnative conformational states, trajectories are generated with the coarse-grained united residue (UNRES) force field. The effectiveness of principal components analysis (PCA), an already established mathematical technique for finding global, correlated motions in atomic simulations of proteins, is evaluated here for coarse-grained trajectories. The problems related to PCA and their solutions are discussed. The folding and nonfolding of proteins are examined with free-energy landscapes. Detailed analyses of many folding and nonfolding trajectories at different temperatures show that PCA is very efficient for characterizing the general folding and nonfolding features of proteins. It is shown that the first principal component captures and describes in detail the dynamics of a system. Anomalous diffusion in the folding/nonfolding dynamics is examined by the mean-square displacement (MSD) and the fractional diffusion and fractional kinetic equations. The collisionless (or ballistic) behavior of a polypeptide undergoing Brownian motion along the first few principal components is accounted for.     

2H-NMR study and molecular dynamics simulation of the location, alignment, and mobility of pyrene in POPC bilayers

The alignment of pyrene in a 1-palmitoyl-2-oleoyl-phosphatidylcholine bilayer was investigated using two different approaches, namely solid-state (2)H-NMR spectroscopy and molecular dynamics (MD) simulations. Quadrupolar splittings from (2)H-NMR spectra of deuterated pyrene-d(10) in an oriented lipid bilayer give information about the orientation of C-D bonds with respect to the membrane normal. From MD simulations, geometric information is accessible via trajectories. By defining molecular and bond order parameters, the data from MD trajectories and NMR spectra can be compared straightforwardly. To ensure that the results from both methods are comparable, parameters of the experimental and the simulation setup were chosen to be as similar as possible. From simulations, we saw that pyrene prefers a position inside the lipid membrane near the headgroups and has no tendency to diffuse from one monolayer of the membrane to the other. The results from simulation and NMR show that the normal of the molecular plane is aligned nearly perpendicular to the bilayer normal. The long axis of pyrene lies preferentially parallel to the bilayer normal within a range of +/-30 degrees . The results from the two different methods are remarkably consistent. The good agreement can be explained by the fact that the different kind of motions of a pyrene molecule are already averaged within a few nanoseconds, which is the timescale covered by the MD simulation.     

A theory of protein dynamics to predict NMR relaxation

We present a theoretical, site-specific, approach to predict protein subunit correlation times, as measured by NMR experiments of ¹H-¹⁵N nuclear Overhauser effect, spin-lattice relaxation, and spin-spin relaxation. Molecular dynamics simulations are input to our equation of motion for protein dynamics, which is solved analytically to produce the eigenvalues and the eigenvectors that specify the NMR parameters. We directly compare our theoretical predictions to experiments and to simulation data for the signal transduction chemotaxis protein Y (CheY), which regulates the swimming response of motile bacteria. Our theoretical results are in good agreement with both simulations and experiments, without recourse to adjustable parameters. The theory is general, since it allows calculations of any dynamical property of interest. As an example, we present theoretical calculations of NMR order parameters and x-ray Debye-Waller temperature factors; both quantities show good agreement with experimental data.     

Molecular dynamics and binding specificity analysis of the bovine immunodeficiency virus BIV Tat-TAR complex

We have performed molecular dynamics (MD) simulations, with particle-mesh Ewald, explicit waters, and counterions, and binding specificity analyses using combined molecular mechanics and continuum solvent (MM-PBSA) on the bovine immunodeficiency virus (BIV) Tat peptide-TAR RNA complex. The solution structure for the complex was solved independently by Patel and co-workers and Puglisi and co-workers. We investigated the differences in both structures and trajectories, particularly in the formation of the U-A-U base triple, the dynamic flexibility of the Tat peptide, and the interactions at the binding interface. We observed a decrease in RMSD in comparing the final average RNA structures and initial RNA structures of both trajectories, which suggests the convergence of the RNA structures to a MD equilibrated RNA structure. We also calculated the relative binding of different Tat peptide mutants to TAR RNA and found qualitative agreement with experimental studies.     

Conformational ensemble of an intrinsically flexible loop in mitochondrial import protein Tim21 studied by modeling and molecular dynamics simulations

BackgroundTim21, a subunit of a highly dynamic translocase of the inner mitochondrial membrane (TIM23) complex, translocates proteins by interacting with subunits in the translocase of the outer membrane (TOM) complex and Tim23 channel in the TIM23 complex. A loop segment in Tim21, which is in close proximity of the binding site of Tim23, has different conformations in X-ray, NMR and new crystal contact-free space (CCFS) structures. MD simulations can provide information on the structure and dynamics of the loop in solution.MethodsThe conformational ensemble of the loop was characterized using loop modeling and molecular dynamics (MD) simulations.ResultsMD simulations confirmed mobility of the loop. Multidimensional scaling and clustering were used to characterize the dynamic conformational ensemble of the loop. Free energy landscape showed that the CCFS crystal structure occupied a low energy region as compared to the conventional X-ray crystal structure. Analysis of crystal packing indicates that the CCFS provides larger conformational space for the motions of the loop.ConclusionsOur work reported the conformational ensemble of the loop in solution, which is in agreement with the structure obtained from CCFS approach. The combination of the experimental techniques and computational methods is beneficial for studying highly flexible regions of proteins.General significanceComputational methods, such as loop modeling and MD simulations, have proved to be useful for studying conformational flexibility of proteins. These methods in integration with experimental techniques such as CCFS has the potential to transform the studies on flexible regions of proteins.     

Fuzzy cluster analysis of molecular dynamics trajectories.

We propose fuzzy clustering as a method to analyze molecular dynamics (MD) trajectories, especially of proteins and polypeptides. A fuzzy cluster analysis locates classes of similar three-dimensional conformations explored during a molecular dynamics simulation. The method can be readily applied to results from both equilibrium and nonequilibrium simulations, with clustering on either global or local structural parameters. The potential of this technique is illustrated by results from fuzzy cluster analyses of trajectories from MD simulations of various fragments of human parathyroid hormone (PTH). For large molecules, it is more efficient to analyze the clustering of root-mean-square distances between conformations comprising the trajectory. We found that the results of the clustering analysis were unambiguous, in terms of the optimal number of clusters of conformations, for the majority of the trajectories examined. The conformation closest to the cluster center can be chosen as being representative of the class of structures making up the cluster, and can be further analyzed, for example, in terms of its secondary structure. The CPU time used by the cluster analysis was negligible compared to the MD simulation time.