A systematic methodology for defining coarse-grained sites in large biomolecules

Coarse-grained (CG) models of biomolecules have recently attracted considerable interest because they enable the simulation of complex biological systems on length-scales and timescales that are inaccessible for atomistic molecular dynamics simulation. A CG model is defined by a map that transforms an atomically detailed configuration into a CG configuration. For CG models of relatively small biomolecules or in cases that the CG and all-atom models have similar resolution, the construction of this map is relatively straightforward and can be guided by chemical intuition. However, it is more challenging to construct a CG map when large and complex domains of biomolecules have to be represented by relatively few CG sites. This work introduces a new and systematic methodology called essential dynamics coarse-graining (ED-CG). This approach constructs a CG map of the primary sequence at a chosen resolution for an arbitrarily complex biomolecule. In particular, the resulting ED-CG method variationally determines the CG sites that reflect the essential dynamics characterized by principal component analysis of an atomistic molecular dynamics trajectory. Numerical calculations illustrate this approach for the HIV-1 CA protein dimer and ATP-bound G-actin. Importantly, since the CG sites are constructed from the primary sequence of the biomolecule, the resulting ED-CG model may be better suited to appropriately explore protein conformational space than those from other CG methods at the same degree of resolution.     

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.     

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.     

An introduction to biomolecular simulations and docking

The biomolecules in and around a living cell – proteins, nucleic acids, lipids and carbohydrates – continuously sample myriad conformational states that are thermally accessible at physiological temperatures. Simultaneously, a given biomolecule also samples (and is sampled by) a rapidly fluctuating local environment comprising other biopolymers, small molecules, water, ions, etc. that diffuse to within a few nanometres, leading to inter-molecular contacts that stitch together large supramolecular assemblies. Indeed, all biological systems can be viewed as dynamic networks of molecular interactions. As a complement to experimentation, molecular simulation offers a uniquely powerful approach to analyse biomolecular structure, mechanism and dynamics; this is possible because the molecular contacts that define a complicated biomolecular system are governed by the same physical principles (forces and energetics) that characterise individual small molecules, and these simpler systems are relatively well-understood. With modern algorithms and computing capabilities, simulations are now an indispensable tool for examining biomolecular assemblies in atomic detail, from the conformational motion in an individual protein to the diffusional dynamics and inter-molecular collisions in the early stages of formation of cellular-scale assemblies such as the ribosome. This text introduces the physicochemical foundations of molecular simulations and docking, largely from the perspective of biomolecular interactions.     

Large‐scale comparison of protein essential dynamics from molecular dynamics simulations and coarse‐grained normal mode analyses

A large‐scale comparison of essential dynamics (ED) modes from molecular dynamic simulations and normal modes from coarse‐grained normal mode methods (CGNM) was performed on a dataset of 335 proteins. As CGNM methods, the elastic network model (ENM) and the rigid cluster normal mode analysis (RCNMA) were used. Low‐frequency normal modes from ENM correlate very well with ED modes in terms of directions of motions and relative amplitudes of motions. Notably, a similar performance was found if normal modes from RCNMA were used, despite a higher level of coarse graining. On average, the space spanned by the first quarter of ENM modes describes 84% of the space spanned by the five ED modes. Furthermore, no prominent differences for ED and CGNM modes among different protein structure classes (CATH classification) were found. This demonstrates the general potential of CGNM approaches for describing intrinsic motions of proteins with little computational cost. For selected cases, CGNM modes were found to be more robust among proteins that have the same topology or are of the same homologous superfamily than ED modes. In view of recent evidence regarding evolutionary conservation of vibrational dynamics, this suggests that ED modes, in some cases, might not be representative of the underlying dynamics that are characteristic of a whole family, probably due to insufficient sampling of some of the family members by MD. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Coarse-grained modeling of the actin filament derived from atomistic-scale simulations

A coarse-grained (CG) procedure that incorporates the information obtained from all-atom molecular dynamics (MD) simulations is presented and applied to actin filaments (F-actin). This procedure matches the averaged values and fluctuations of the effective internal coordinates that are used to define a CG model to the values extracted from atomistic MD simulations. The fluctuations of effective internal coordinates in a CG model are computed via normal-mode analysis (NMA), and the computed fluctuations are matched with the atomistic MD results in a self-consistent manner. Each actin monomer (G-actin) is coarse-grained into four sites, and each site corresponds to one of the subdomains of G-actin. The potential energy of a CG G-actin contains three bonds, two angles, and one dihedral angle; effective harmonic bonds are used to describe the intermonomer interactions in a CG F-actin. The persistence length of a CG F-actin was found to be sensitive to the cut-off distance of assigning intermonomer bonds. Effective harmonic bonds for a monomer with its third nearest neighboring monomers are found to be necessary to reproduce the values of persistence length obtained from all-atom MD simulations. Compared to the elastic network model, incorporating the information of internal coordinate fluctuations enhances the accuracy and robustness for a CG model to describe the shapes of low-frequency vibrational modes. Combining the fluctuation-matching CG procedure and NMA, the achievable time- and length scales of modeling actin filaments can be greatly enhanced. In particular, a method is described to compute the force-extension curve using the CG model developed in this work and NMA. It was found that F-actin is easily buckled under compressive deformation, and a writhing mode is developed as a result. In addition to the bending and twisting modes, this novel writhing mode of F-actin could also play important roles in the interactions of F-actin with actin-binding proteins and in the force-generation process via polymerization.     

Computational simulations of TNF receptor oligomerization on plasma membrane

The interactions between tumor necrosis factors (TNFs) and their corresponding receptors (TNFRs) play a pivotal role in inflammatory responses. Upon ligand binding, TNFR receptors were found to form oligomers on cell surfaces. However, the underlying mechanism of oligomerization is not fully understood. In order to tackle this problem, molecular dynamics (MD) simulations have been applied to the complex between TNF receptor-1 (TNFR1) and its ligand TNF-α as a specific test system. The simulations on both all-atom (AA) and coarse-grained (CG) levels achieved the similar results that the extracellular domains of TNFR1 can undergo large fluctuations on plasma membrane, while the dynamics of TNFα-TNFR1 complex is much more constrained. Using the CG model with the Martini force field, we are able to simulate the systems that contain multiple TNFα-TNFR1 complexes with the timescale of microseconds. We found that complexes can aggregate into oligomers on the plasma membrane through the lateral interactions between receptors at the end of the CG simulations. We suggest that this spatial organization is essential to the efficiency of signal transduction for ligands that belong to the TNF superfamily. We further show that the aggregation of two complexes is initiated by the association between the N-terminal domains of TNFR1 receptors. Interestingly, the cis-interfaces between N-terminal regions of two TNF receptors have been observed in the previous X-ray crystallographic experiment. Therefore, we provide supportive evidence that cis-interface is of functional importance in triggering the receptor oligomerization. Taken together, our study brings insights to understand the molecular mechanism of TNF signaling.     

REACH coarse-grained biomolecular simulation: transferability between different protein structural classes

Coarse graining of protein interactions provides a means of simulating large biological systems. The REACH (Realistic Extension Algorithm via Covariance Hessian) coarse-graining method, in which the force constants of a residue-scale elastic network model are calculated from the variance-covariance matrix obtained from atomistic molecular dynamics (MD) simulation, involves direct mapping between scales without the need for iterative optimization. Here, the transferability of the REACH force field is examined between protein molecules of different structural classes. As test cases, myoglobin (all α), plastocyanin (all β), and dihydrofolate reductase (α/β) are taken. The force constants derived are found to be closely similar in all three proteins. An MD version of REACH is presented, and low-temperature coarse-grained (CG) REACH MD simulations of the three proteins are compared with atomistic MD results. The mean-square fluctuations of the atomistic MD are well reproduced by the CGMD. Model functions for the CG interactions, derived by averaging over the three proteins, are also shown to produce fluctuations in good agreement with the atomistic MD. The results indicate that, similarly to the use of atomistic force fields, it is now possible to use a single, generic REACH force field for all protein studies, without having first to derive parameters from atomistic MD simulation for each individual system studied. The REACH method is thus likely to be a reliable way of determining spatiotemporal motion of a variety of proteins without the need for expensive computation of long atomistic MD simulations.     

RedMDStream: Parameterization and Simulation Toolbox for Coarse-Grained Molecular Dynamics Models

Coarse-grained (CG) models in molecular dynamics (MD) are powerful tools to simulate the dynamics of large biomolecular systems on micro- to millisecond timescales. However, the CG model, potential energy terms, and parameters are typically not transferable between different molecules and problems. So parameterizing CG force fields, which is both tedious and time-consuming, is often necessary. We present RedMDStream, a software for developing, testing, and simulating biomolecules with CG MD models. Development includes an automatic procedure for the optimization of potential energy parameters based on metaheuristic methods. As an example we describe the parameterization of a simple CG MD model of an RNA hairpin.     

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.     

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.     

Principal component and normal mode analysis of proteins; a quantitative comparison using the GroEL subunit

Principal component analysis (PCA) and normal mode analysis (NMA) have emerged as two invaluable tools for studying conformational changes in proteins. To compare these approaches for studying protein dynamics, we have used a subunit of the GroEL chaperone, whose dynamics is well characterized. We first show that both PCA on trajectories from molecular dynamics (MD) simulations and NMA reveal a general dynamical behavior in agreement with what has previously been described for GroEL. We thus compare the reproducibility of PCA on independent MD runs and subsequently investigate the influence of the length of the MD simulations. We show that there is a relatively poor one-to-one correspondence between eigenvectors obtained from two independent runs and conclude that caution should be taken when analyzing principal components individually. We also observe that increasing the simulation length does not improve the agreement with the experimental structural difference. In fact, relatively short MD simulations are sufficient for this purpose. We observe a rapid convergence of the eigenvectors (after ca. 6 ns). Although there is not always a clear one-to-one correspondence, there is a qualitatively good agreement between the movements described by the first five modes obtained with the three different approaches; PCA, all-atoms NMA, and coarse-grained NMA. It is particularly interesting to relate this to the computational cost of the three methods. The results we obtain on the GroEL subunit contribute to the generalization of robust and reproducible strategies for the study of protein dynamics, using either NMA or PCA of trajectories from MD simulations.     

Elastic network models for understanding biomolecular machinery: from enzymes to supramolecular assemblies

With advances in structure genomics, it is now recognized that knowledge of structure alone is insufficient to understand and control the mechanisms of biomolecular function. Additional information in the form of dynamics is needed. As demonstrated in a large number of studies, the machinery of proteins and their complexes can be understood to a good approximation by adopting Gaussian (or elastic) network models (GNM) for simplified normal mode analyses. While this approximation lacks chemical details, it provides us with a means for assessing the collective motions of large structures/assemblies and perform a comparative analysis of a series of proteins, thus providing insights into the mechanical aspects of biomolecular dynamics. In this paper, we discuss recent applications of GNM to a series of enzymes as well as large structures such as the HK97 bacteriophage viral capsids. Understanding the dynamics of large protein structures can be computationally challenging. To this end, we introduce a new approach for building a hierarchical, reduced rank representation of the protein topology and consequently the fluctuation dynamics.     

A flexible docking scheme to explore the binding selectivity of PDZ domains

Modeling of protein binding site flexibility in molecular docking is still a challenging problem due to the large conformational space that needs sampling. Here, we propose a flexible receptor docking scheme: A dihedral restrained replica exchange molecular dynamics (REMD), where we incorporate the normal modes obtained by the Elastic Network Model (ENM) as dihedral restraints to speed up the search towards correct binding site conformations. To our knowledge, this is the first approach that uses ENM modes to bias REMD simulations towards binding induced fluctuations in docking studies. In our docking scheme, we first obtain the deformed structures of the unbound protein as initial conformations by moving along the binding fluctuation mode, and perform REMD using the ENM modes as dihedral restraints. Then, we generate an ensemble of multiple receptor conformations (MRCs) by clustering the lowest replica trajectory. Using ROSETTA LIGAND , we dock ligands to the clustered conformations to predict the binding pose and affinity. We apply this method to postsynaptic density‐95/Dlg/ZO‐1 (PDZ) domains; whose dynamics govern their binding specificity. Our approach produces the lowest energy bound complexes with an average ligand root mean square deviation of 0.36 Å. We further test our method on (i) homologs and (ii) mutant structures of PDZ where mutations alter the binding selectivity. In both cases, our approach succeeds to predict the correct pose and the affinity of binding peptides. Overall, with this approach, we generate an ensemble of MRCs that leads to predict the binding poses and specificities of a protein complex accurately.     

NbIT - A New Information Theory-Based Analysis of Allosteric Mechanisms Reveals Residues that Underlie Function in the Leucine Transporter LeuT

Complex networks of interacting residues and microdomains in the structures of biomolecular systems underlie the reliable propagation of information from an input signal, such as the concentration of a ligand, to sites that generate the appropriate output signal, such as enzymatic activity. This information transduction often carries the signal across relatively large distances at the molecular scale in a form of allostery that is essential for the physiological functions performed by biomolecules. While allosteric behaviors have been documented from experiments and computation, the mechanism of this form of allostery proved difficult to identify at the molecular level. Here, we introduce a novel analysis framework, called N-body Information Theory (NbIT) analysis, which is based on information theory and uses measures of configurational entropy in a biomolecular system to identify microdomains and individual residues that act as (i)-channels for long-distance information sharing between functional sites, and (ii)-coordinators that organize dynamics within functional sites. Application of the new method to molecular dynamics (MD) trajectories of the occluded state of the bacterial leucine transporter LeuT identifies a channel of allosteric coupling between the functionally important intracellular gate and the substrate binding sites known to modulate it. NbIT analysis is shown also to differentiate residues involved primarily in stabilizing the functional sites, from those that contribute to allosteric couplings between sites. NbIT analysis of MD data thus reveals rigorous mechanistic elements of allostery underlying the dynamics of biomolecular systems.