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.     

Characterization of NMR relaxation-active motions of a partially folded A-state analogue of ubiquitin

The dominant dynamics of a partially folded A-state analogue of ubiquitin that give rise to NMR 15N spin relaxation have been investigated using molecular dynamics (MD) computer simulations and reorientational quasiharmonic analysis. Starting from the X-ray structure of native ubiquitin with a protonation state corresponding to a low pH, the A-state analogue was generated by a MD simulation of a total length of 33 ns in a 60%/40% methanol/water mixture using a variable temperature scheme to control and speed up the structural transformation. The N-terminal half of the A-state analogue consists of loosely coupled native-like secondary structural elements, while the C-terminal half is mostly irregular in structure. Analysis of dipolar N-H backbone correlation functions reveals reorientational amplitudes and time-scale distributions that are comparable to those observed experimentally. Thus, the trajectory provides a realistic picture of a partially folded protein that can be used for gaining a better understanding of the various types of reorientational motions that are manifested in spin-relaxation parameters of partially folded systems. For this purpose, a reorientational quasiharmonic reorientational analysis was performed on the final 5 ns of the trajectory of the A-state analogue, and for comparison on a 5 ns trajectory of native ubiquitin. The largest amplitude reorientational modes show a markedly distinct behavior for the two states. While for native ubiquitin, such motions have a more local character involving loops and the C-terminal end of the polypeptide chain, the A-state analogue shows highly collective motions in the nanosecond time-scale range corresponding to larger-scale movements between different segments. Changes in reorientational backbone entropy between the A-state analogue and the native state of ubiquitin, which were computed from the reorientational quasiharmonic analyses, are found to depend significantly on motional correlation effects.     

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.     

Close correspondence between the motions from principal component analysis of multiple HIV-1 protease structures and elastic network modes

The large number of available HIV-1 protease structures provides a remarkable sampling of conformations of the different conformational states, which can be viewed as direct structural information about the dynamics of the HIV-1 protease. After structure matching, we apply principal component analysis (PCA) to obtain the important apparent motions for both bound and unbound structures. There are significant similarities between the first few key motions and the first few low-frequency normal modes calculated from a static representative structure with an elastic network model (ENM), strongly suggesting that the variations among the observed structures and the corresponding conformational changes are facilitated by the low-frequency, global motions intrinsic to the structure. Similarities are also found when the approach is applied to an NMR ensemble, as well as to molecular dynamics (MD) trajectories. Thus, a sufficiently large number of experimental structures can directly provide important information about protein dynamics, but ENM can also provide similar sampling of conformations.     

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.     

Quantitative delineation of how breathing motions open ligand migration channels in myoglobin and its mutants

Ligand migration and binding are central to the biological functions of many proteins such as myoglobin (Mb) and it is widely thought that protein breathing motions open up ligand channels dynamically. However, how a protein exerts its control over the opening and closing of these channels through its intrinsic dynamics is not fully understood. Specifically, a quantitative delineation of the breathing motions that are needed to open ligand channels is lacking. In this work, we present and apply a novel normal mode‐based method to quantitatively delineate what and how breathing motions open ligand migration channels in Mb and its mutants. The motivation behind this work springs from the observation that normal mode motions are closely linked to the breathing motions that are thought to open ligand migration channels. In addition, the method provides a direct and detailed depiction of the motions of each and every residue that lines a channel and can identify key residues that play a dominating role in regulating the channel. The all‐atom model and the full force‐field employed in the method provide a realistic energetics on the work cost required to open a channel, and as a result, the method can be used to efficiently study the effects of mutations on ligand migration channels and on ligand entry rates. Our results on Mb and its mutants are in excellent agreement with MD simulation results and experimentally determined ligand entry rates. Proteins 2015; 83:757–770. © 2015 Wiley Periodicals, Inc.     

Studying the collective motions of the adenosine A2A receptor as a result of ligand binding using principal component analysis

Abstract

Adenosine receptors (ARs) belong to family A of GPCRs that are involved in many diseases, including cerebral and cardiac ischemic diseases, immune and inflammatory disorders, etc. Thus, they represent important therapeutic targets to treat these conditions. Computational techniques such as molecular dynamics (MD) simulations permit researchers to obtain structural information about these proteins, and principal component analysis (PCA) allows for the identification of collective motions. There are available structures for the active form (3QAK) and the inactive form (3EML) of A2AR which permit us to gain insight about their activation/inactivation mechanism. In this work, we have proposed an inverse strategy using MD simulations where the active form was coupled to the antagonist caffeine and the inactive form was coupled to adenosine agonist. Moreover, we have included four reported thermostabilizing mutations in the inactive form to study A2AR structural differences under different conditions. Some observations stand out from the PCA studies. For instance, the apo structures showed remarkable similarities, and the principal components (PCs) were rearranged in a ligand-dependent manner. Additionally, the active conformation was less stable compared to the inactive one. Some PCs inverted their direction in the presence of a ligand, and comparison of the PCs between 3EML and 3EML_ADN showed that adenosine induced major changes in the structure of A2AR. Rearrangement of PCs precedes and drives conformational changes that occur after ligand binding. Knowledge about these conformational changes provides important insights about the activity of A2AR.     

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.     

Simulating protein motions with rigidity analysis.

Protein motions, ranging from molecular flexibility to large-scale conformational change, play an essential role in many biochemical processes. Despite the explosion in our knowledge of structural and functional data, our understanding of protein movement is still very limited. In previous work, we developed and validated a motion planning based method for mapping protein folding pathways from unstructured conformations to the native state. In this paper, we propose a novel method based on rigidity theory to sample conformation space more effectively, and we describe extensions of our framework to automate the process and to map transitions between specified conformations. Our results show that these additions both improve the accuracy of our maps and enable us to study a broader range of motions for larger proteins. For example, we show that rigidity-based sampling results in maps that capture subtle folding differences between protein G and its mutants, NuG1 and NuG2, and we illustrate how our technique can be used to study large-scale conformational changes in calmodulin, a 148 residue signaling protein known to undergo conformational changes when binding to Ca(2+). Finally, we announce our web-based protein folding server which includes a publicly available archive of protein motions: (http://parasol.tamu.edu/foldingserver/).     

Domain motions of the Mip protein from Legionella pneumophila

The homodimeric 45.6 kDa (total mass) Mip protein, a virulence factor from Legionella pneumophila, was investigated with solution NMR spectroscopy and molecular dynamics (MD) simulations. Two Mip monomers are dimerized via an N-terminal helix bundle that is connected via a long alpha-helix to a C-terminal FKBP domain in each subunit. More than 85% of the amino acids were identified in triple-resonance NMR spectra. (15)N relaxation analysis showed a bimodal distribution of R(1)/R(2) values, with the lower ratio in the N-terminal domain. Relaxation dispersion measurements confirmed that these reduced ratios did not originate from conformational exchange. Thus, two different correlation times (tau(c)) can be deduced, reflecting partly uncoupled motions of both domains. Relaxation data of a Mip(77)(-)(213) monomer mutant were similar to those observed in the dimer, corroborating that the FKBP domain, including part of the connecting helix, behaves as one dynamic entity. MD simulations (18 ns) of the Mip dimer also yielded two different correlation times for the two domains and thus confirm the independence of the domain motions. Principal component analysis of the dihedral space covariance matrix calculated from the MD trajectory suggests a flexible region in the long connecting helix that acts as a hinge between the two domains. Such motion provides a possible explanation of how Mip can bind to complex molecular components of the extracellular matrix and mediate alveolar damage and bacterial spread in the lung.     

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.     

Functional domain motions in proteins on the ~1-100 ns timescale: comparison of neutron spin-echo spectroscopy of phosphoglycerate kinase with molecular-dynamics simulation

Protein function often requires large-scale domain motion. An exciting new development in the experimental characterization of domain motions in proteins is the application of neutron spin-echo spectroscopy (NSE). NSE directly probes coherent (i.e., pair correlated) scattering on the ~1-100 ns timescale. Here, we report on all-atom molecular-dynamics (MD) simulation of a protein, phosphoglycerate kinase, from which we calculate small-angle neutron scattering (SANS) and NSE scattering properties. The simulation-derived and experimental-solution SANS results are in excellent agreement. The contributions of translational and rotational whole-molecule diffusion to the simulation-derived NSE and potential problems in their estimation are examined. Principal component analysis identifies types of domain motion that dominate the internal motion's contribution to the NSE signal, with the largest being classic hinge bending. The associated free-energy profiles are quasiharmonic and the frictional properties correspond to highly overdamped motion. The amplitudes of the motions derived by MD are smaller than those derived from the experimental analysis, and possible reasons for this difference are discussed. The MD results confirm that a significant component of the NSE arises from internal dynamics. They also demonstrate that the combination of NSE with MD is potentially useful for determining the forms, potentials of mean force, and time dependence of functional domain motions in proteins.     

Domain motions of hyaluronan lyase underlying processive hyaluronan translocation

Hyaluronan lyase (Hyal) is a surface enzyme occurring in many bacterial organisms including members of Streptococcus species. Streptococcal Hyal primarily degrades hyaluronan‐substrate (HA) of the extracellular matrix. This degradation appears to facilitate the spread of this bacterium throughout host tissues. Unlike purely endolytic degradation of its other substrates, unsulfated chondroitin or some chondroitin sulfates, the degradation of HA by Hyal proceeds by processive exolytic cleavage of one disaccharide at a time following an initial endolytic cut. Molecular dynamics (MD) studies of Hyal from Streptococcus pneumoniae are presented that address the enzyme's molecular mechanism of action and the role of domain motions for processive functionality. The analysis of extensive sub‐microsecond MD simulations of this enzyme action on HA‐substrates of different lengths and the connection between the domain dynamics of Hyal and the translocation of the HA‐substrate reveals that opening/closing and twisting domain motions of the Hyal are intimately linked to processive HA degradation. Enforced simulations confirmed this finding as the domain motions in SpnHyal were found to be induced by enforced substrate translocation. These results establish the dynamic interplay between Hyal flexibility and substrate translocation and provide insight into the processive mechanism of Hyal. Proteins 2009. © 2008 Wiley‐Liss, Inc.     

Thermodynamic dissection of large‐scale domain motions coupled with ligand binding of enzyme I

Domain motions are central to the biological functions of many proteins. The energetics of the motions, however, is often difficult to characterize when motions are coupled with the ligand binding. Here, we determined the thermodynamic parameters of individual domain motions and ligand binding of enzyme I (EI) using strategic domain-deletion mutants that selectively removed particular motions. Upon ligand binding, EI employs two large-scale domain motions, the hinge motion and the swivel motion, to switch between conformational states of distinct domain−domain orientations. Calorimetric analysis of the EI mutants separated the free energy changes of the binding and motions, demonstrating that the unfavorable hinge motion (ΔG = 1.5 kcal mol⁻¹) was driven by the favorable swivel motion (ΔG = −5.2 kcal mol⁻¹). The large free energy differences could be explained by the physicochemical nature of the domain interfaces associated with the motions; the hinge motion employed much narrower interface than the swivel motion without any hydrogen bonds or salt bridges. The small heat capacity further suggested that the packing of the domain interfaces associated with the hinge motion was less compact than that commonly observed in proteins. Lastly, thermodynamic analysis of phosphorylated EI suggests that the domain motions are regulated by the ligand binding and the phosphorylation states. Taken together, the thermodynamic dissection approach illustrates how multiple motions and ligand binding are energetically connected during the functional cycle of EI.     

KCTD5 is endowed with large,functionally relevant,interdomain motions

The KCTD family is an emerging class of proteins that are involved in important biological processes whose biochemical and structural properties are rather poorly characterized or even completely undefined. We here used KCTD5, the only member of the family with a known three-dimensional structure, to gain insights into the intrinsic structural stability of the C-terminal domain (CTD) and into the mutual dynamic interplay between the two domains of the protein. Molecular dynamics (MD) simulations indicate that in the simulation timescale (120 ns), the pentameric assembly of the CTD is endowed with a significant intrinsic stability. Moreover, MD analyses also led to the identification of exposed β-strand residues. Being these regions intrinsically sticky, they could be involved in the substrate recognition. More importantly, simulations conducted on the full-length protein provide interesting information of the relative motions between the BTB domain and the CTD of the protein. Indeed, the dissection of the overall motion of the protein is indicative of a large interdomain twisting associated with limited bending movements. Notably, MD data indicate that the entire interdomain motion is pivoted by a single residue (Ser150) of the hinge region that connects the domains. The functional relevance of these motions was evaluated in the context of the functional macromolecular machinery in which KCTD5 is involved. This analysis indicates that the interdomain twisting motion here characterized may be important for the correct positioning of the substrate to be ubiquitinated with respect to the other factors of the ubiquitination machinery.     

Effects of ligand binding and oxidation on hinge-bending motions in S-adenosyl-L-homocysteine hydrolase

Domain motions of S-adenosyl-l-homocysteine (AdoHcy) hydrolase have been detected by time-resolved fluorescence anisotropy measurements. Time constants for reorientational motions in the native enzyme were compared with those for enzymes where key residues were altered by site-directed mutation. Mutations M351P, H353A, and P354A were selected in a hinge region for motion between the open and closed forms of the enzyme, as identified in a previous normal-mode study [Wang et al. (2005) Domain motions and the open-to-closed conformational transition of an enzyme: A normal-mode analysis of S-adenosyl-l-homocysteine hydrolase, Biochemistry 44, 7228-7239]. In wild-type, substrate-free AdoHcy hydrolase (NAD(+) cofactor in each subunit), reorientational motions were detected on time scales of 10-20 and 80-90 ns. The faster motion is attributed to the domain motion, and the slower motion is attributed to the tumbling of the enzyme. The domain motion was also detected for the enzyme complexes E(NADH/3'-keto-adenosine) and E(NAD(+)/3'-deoxyadenosine) but was absent for the complex E(NADH/3'-keto-neplanocin A). The results indicate that AdoHcy hydrolase exists in equilibrium of open and closed structures, with the equilibrium shifted toward the more mobile open form for the substrate-free enzyme, E(NAD(+)), and for intermediates formed early in the catalytic cycle after substrate binding or formed late prior to product release, E(NAD(+)/ligand). However, the strong inhibitor neplanocin A upon binding undergoes oxidation, forming the complex E(NADH/3'-keto-neplanocin). For this complex, which is analogous to the enzyme complex with the central catalytic intermediate, the equilibrium was shifted toward the more rigid closed form. A similar pattern was observed for M351P and P354A mutants. In contrast, the domain motion could not be detected, either in the absence or presence of ligands or with the cofactor in either the oxidized or reduced state, for the H353A protein, suggesting that this mutation changes the hinge-bending dynamics of the enzyme.     

Insight derived from molecular dynamics simulation into dynamics and molecular motions of cuticle-degrading serine protease Ver112

Cuticle-degrading serine protease Ver112, which derived from a nematophagous fungus Lecanicillium psalliotae, has been exhibited to have high cuticle-degrading and nematicidal activities. We have performed molecular dynamics (MD) simulation based on the crystal structure of Ver112 to investigate its dynamic properties and large-scale concerted motions. The results indicate that the structural core of Ver112 shows a small fluctuation amplitude, whereas the substrate binding sites, and the regions close to and opposite the substrate binding sites experience significant conformational fluctuations. The large concerted motions obtained from essential dynamics (ED) analysis of MD trajectory can lead to open or close of the substrate binding sites, which are proposed to be linked to the functional properties of Ver112, such as substrate binding, orientation, catalytic, and release. The significant motion in the loop regions that is located opposite the binding sites are considered to play an important role in modulating the dynamics of the substrate binding sites. Furthermore, the bottom of free energy landscape (FEL) of Ver112 are rugged, which is mainly caused by the fluctuations of substrate binding regions and loops located opposite the binding site. In addition, the mechanism underlying the high flexibility and catalytic activity of Ver112 was also discussed. Our simulation study complements the biochemical and structural studies, and provides insight into the dynamics-function relationship of cuticle-degrading serine protease Ver112.     

A comparative analysis of the equilibrium dynamics of a designed protein inferred from NMR,X‐ray,and computations

A detailed analysis of high‐resolution structural data and computationally predicted dynamics was carried out for a designed sugar‐binding protein. The mean‐square deviations in the positions of residues derived from nuclear magnetic resonance (NMR) models and those inferred from X‐ray crystallographic B‐factors for two different crystal forms were compared with the predictions based on the Gaussian Network Model (GNM) and the results from molecular dynamics (MD) simulations. GNM systematically yielded a higher correlation than MD, with experimental data, suggesting that the lack of atomistic details in the coarse‐grained GNM is more than compensated for by the mathematically exact evaluation of fluctuations using the native contacts topology. Evidence is provided that particular loop motions are curtailed by intermolecular contacts in the crystal environment causing a discrepancy between theory and experiments. Interestingly, the information conveyed by X‐ray crystallography becomes more consistent with NMR models and computational predictions when ensembles of X‐ray models are considered. Less precise (broadly distributed) ensembles indeed appear to describe the accessible conformational space under native state conditions better than B‐factors. Our results highlight the importance of using multiple conformations obtained by alternative experimental methods, and analyzing results from both coarse‐grained models and atomic simulations, for accurate assessment of motions accessible to proteins under native state conditions. Proteins 2009. © 2009 Wiley‐Liss, Inc.     

Comprehensive analysis of motions in molecular dynamics trajectories of the actin capping protein and its inhibitor complexes

The actin capping protein (CP) binds to actin filaments to block further elongation. The capping activity is inhibited by proteins V‐1 and CARMIL interacting with CP via steric and allosteric mechanisms, respectively. The crystal structures of free CP, CP/V‐1, and CP/CARMIL complexes suggest that the binding of CARMIL alters the flexibility of CP rather than the overall structure of CP, and this is an allosteric inhibition mechanism. Here, we performed molecular dynamics (MD) simulations of CP in the free form, and in complex with CARMIL or V‐1. The resulting trajectories were analyzed exhaustively using Motion Tree, which identifies various rigid‐body motions ranging from small local motions to large domain motions. After enumerating all the motions, CP flexibilities with different ligands were characterized by a list of frequencies for 20 dominant rigid‐body motions, some of which were not identified in previous studies. The comparative analysis highlights the influence of the binding of the CARMIL peptide to CP flexibility. In free CP and the CP/V‐1 complex, domain motions around a large crevice between the N‐stalk and the CP‐S domain occur frequently. The CARMIL peptide binds the crevice and suppresses the motions effectively. In addition, the binding of the CARMIL peptide enhances and alters local motions around the pocket that participates in V‐1 binding. These newly identified motions are likely to suppress the binding of V‐1 to CP. The observed changes in CP motion provide insights that describe the mechanism of allosteric regulation by CARMIL through modulating CP flexibility. Proteins 2016; 84:948–956. © 2016 Wiley Periodicals, Inc.