Computation of low-frequency normal modes in macromolecules: improvements to the method of diagonalization in a mixed basis and application to hemoglobin

Scalable computation of anisotropic vibrations for large macromolecular assemblies

The Normal Mode Analysis (NMA) is a standard approach to elucidate the anisotropic vibrations of macromolecules at their folded states, where low-frequency collective motions can reveal rearrangements of domains and changes in the exposed surface of macromolecules. Recent advances in structural biology have enabled the resolution of megascale macromolecules with millions of atoms. However, the calculation of their vibrational modes remains elusive due to the prohibitive cost associated with constructing and diagonalizing the underlying eigenproblem and the current approaches to NMA are not readily adaptable for efficient parallel computing on graphic processing unit (GPU). Here, we present eigenproblem construction and diagonalization approach that implements level-structure bandwidth-reducing algorithms to transform the sparse computation in NMA to a globally-sparse-yet-locally-dense computation, allowing batched tensor products to be most efficiently executed on GPU. We map, optimize, and compare several low-complexity Krylov-subspace eigensolvers, supplemented by techniques such as Chebyshev filtering, sum decomposition, external explicit deflation and shift-and-inverse, to allow fast GPU-resident calculations. The method allows accurate calculation of the first 1000 vibrational modes of some largest structures in PDB ( > 2.4 million atoms) at least 250 times faster than existing methods.

Unveiling mutation effects on the structural dynamics of the main protease from SARS-CoV-2 with hybrid simulation methods

The main protease of SARS-CoV-2 (called M^pro or 3CL^pro) is essential for processing polyproteins encoded by viral RNA. Several M^pro mutations were found in SARS-CoV-2 variants, which are related to higher transmissibility, pathogenicity, and resistance to neutralization antibodies. Macromolecules adopt several favored conformations in solution depending on their structure and shape, determining their dynamics and function. In this study, we used a hybrid simulation method to generate intermediate structures along the six lowest frequency normal modes and sample the conformational space and characterize the structural dynamics and global motions of WT SARS-CoV-2 M^pro and 48 mutations, including mutations found in P.1, B.1.1.7, B.1.351, B.1.525 and B.1.429+B.1.427 variants. We tried to contribute to the elucidation of the effects of mutation in the structural dynamics of SARS-CoV-2 M^pro. A machine learning analysis was performed following the investigation regarding the influence of the K90R, P99L, P108S, and N151D mutations on the dimeric interface assembling of the SARS-CoV-2 M^pro. The parameters allowed the selection of potential structurally stable dimers, which demonstrated that some single surface aa substitutions not located at the dimeric interface (K90R, P99L, P108S, and N151D) are able to induce significant quaternary changes. Furthermore, our results demonstrated, by a Quantum Mechanics method, the influence of SARS-CoV-2 M^pro mutations on the catalytic mechanism, confirming that only one of the chains of the WT and mutant SARS-CoV-2 M^pros are prone to cleave substrates. Finally, it was also possible to identify the aa residue F140 as an important factor related to the increasing enzymatic reactivity of a significant number of SARS-CoV-2 M^pro conformations generated by the normal modes-based simulations.

Effect of phosphorylation on the structural dynamics, thermal stability of human dopamine transporter: A simulation study using normal modes, molecular dynamics and Markov State Model

The Human Dopamine Transporter (hDAT) plays an essential role in modulating the Influx/Efflux of dopamine, and it is involved in the mechanism of certain neurodegenerative diseases such as Parkinson's disease. Several studies have reported important states for Dopamine transport: outward-facing open state (OFo), the outward-facing closed state (OFc), the holo-occluded state closed (holo), and the inward-facing open state (IFo). Furthermore, experimental assays have shown that different phosphorylation conditions in hDAT can affect the rate of dopamine absorption. We present a protocol using hybrid simulation methods to study the conformational dynamics and stability of states of hDAT under different phosphorylation sites. With this protocol, we explored the conformational space of hDAT, identified the states, and evaluated the free energy differences and the transition probabilities between them in each of the phosphorylation cases. We also presented the conformational changes and correlated them with those described in the literature. There is a thesis/hypothesis that the phosphorylation condition corresponding to NP-333 system (where all sites Ser/Thr from residue 2 to 62 and 254 to 613 are phosphorylated, except residue 333) would decrease the rate of dopamine transport from the extracellular medium to the intracellular medium by hDAT as previously described in the literature by Lin et al., 2003. Our results corroborated this thesis/hypothesis and the data reported. It is probably due to the affectation/changes/alteration of the conformational dynamics of this system that makes the intermediate states more likely and makes it difficult to initial states associated with the uptake of dopamine in the extracellular medium, corroborating the experimental results. Furthermore, our results showed that just single phosphorylation/dephosphorylation could alter intrinsic protein motions affecting the sampling of one or more states necessary for dopamine transport. In this sense, the modification of phosphorylation influences protein movements and conformational preferences, affecting the stability of states and the transition between them and, therefore, the transport.

Cryo-electron microscopy unveils unique structural features of the human Kir2.1 channel

We present the first structure of the human Kir2.1 channel containing both transmembrane domain (TMD) and cytoplasmic domain (CTD). Kir2.1 channels are strongly inward-rectifying potassium channels that play a key role in maintaining resting membrane potential. Their gating is modulated by phosphatidylinositol 4,5-bisphosphate (PIP₂). Genetically inherited defects in Kir2.1 channels are responsible for several rare human diseases, including Andersen's syndrome. The structural analysis (cryo-electron microscopy), surface plasmon resonance, and electrophysiological experiments revealed a well-connected network of interactions between the PIP₂-binding site and the G-loop through residues R312 and H221. In addition, molecular dynamics simulations and normal mode analysis showed the intrinsic tendency of the CTD to tether to the TMD and a movement of the secondary anionic binding site to the membrane even without PIP₂. Our results revealed structural features unique to human Kir2.1 and provided insights into the connection between G-loop and gating and the pathological mechanisms associated with this channel.

ABCG2/BCRP transport mechanism revealed through kinetically excited targeted molecular dynamics simulations

ABCG2/BCRP is an ABC transporter that plays an important role in tissue protection by exporting endogenous substrates and xenobiotics. ABCG2 is of major interest due to its involvement in multidrug resistance (MDR), and understanding its complex efflux mechanism is essential to preventing MDR and drug-drug interactions (DDI). ABCG2 export is characterized by two major conformational transitions between inward- and outward-facing states, the structures of which have been resolved. Yet, the entire transport cycle has not been characterized to date. Our study bridges the gap between the two extreme conformations by studying connecting pathways. We developed an innovative approach to enhance molecular dynamics simulations, ‘kinetically excited targeted molecular dynamics’, and successfully simulated the transitions between inward- and outward-facing states in both directions and the transport of the endogenous substrate estrone 3-sulfate. We discovered an additional pocket between the two substrate-binding cavities and found that the presence of the substrate in the first cavity is essential to couple the movements between the nucleotide-binding and transmembrane domains. Our study shed new light on the complex efflux mechanism, and we provided transition pathways that can help to identify novel substrates and inhibitors of ABCG2 and probe new drug candidates for MDR and DDI.

Normal-mode driven exploration of protein domain motions

Domain motions involved in the function of proteins can often be well described as a combination of motions along a handfull of low-frequency modes, that is, with the values of a few normal coordinates. This means that, when the functional motion of a protein is unknown, it should prove possible to predict it, since it amounts to guess a few values. However, without the help of additional experimental data, using normal coordinates for generating accurate conformers far away from the initial one is not so straightforward. To do so, a new approach is proposed: instead of building conformers directly with the values of a subset of normal coordinates, they are built in two steps, the conformer built with normal coordinates being just used for defining a set of distance constraints, the final conformer being built so as to match them. Note that this approach amounts to transform the problem of generating accurate protein conformers using normal coordinates into a better known one: the distance-geometry problem, which is herein solved with the help of the ROSETTA software. In the present study, this approach allowed to rebuild accurately six large amplitude conformational changes, using at most six low-frequency normal coordinates. As a consequence of the low-dimensionality of the corresponding subspace, random exploration also proved enough for generating low-energy conformers close to the known end-point of the conformational change of the LAO binding protein, lysozyme T4 and adenylate kinase.

Unexpected Gating Behaviour of an Engineered Potassium Channel Kir

In this study, we investigated the dynamics and functional characteristics of the KirBac3.1 S129R, a mutated bacterial potassium channel for which the inner pore-lining helix (TM2) was engineered so that the bundle crossing is trapped in an open conformation. The structure of this channel has been previously determined at high atomic resolution. We explored the dynamical characteristics of this open state channel using an in silico method MDeNM that combines molecular dynamics simulations and normal modes. We captured the global and local motions at the mutation level and compared these data with HDX-MS experiments. MDeNM provided also an estimation of the probability of the different opening states that are in agreement with our electrophysiological experiments. In the S129R mutant, the Arg129 mutation releases the two constriction points in the channel that existed in the wild type but interestingly creates another restriction point.

Deciphering the Mechanism of Inhibition of SERCA1a by Sarcolipin Using Molecular Simulations

SERCA1a is an ATPase calcium pump that transports Ca²⁺ from the cytoplasm to the sarco/endoplasmic reticulum lumen. Sarcolipin (SLN), a transmembrane peptide, regulates the activity of SERCA1a by decreasing its Ca²⁺ transport rate, but its mechanism of action is still not well-understood. To decipher this mechanism, we have performed normal mode analysis in the all-atom model, with the SERCA1a-SLN complex, or the isolated SERCA1a, embedded in an explicit membrane. The comparison of the results allowed us to provide an explanation at the atomic level for the action of SLN that is in good agreement with experimental observations. In our analyses, the presence of SLN locally perturbs the TM6 transmembrane helix and as a consequence modifies the position of D800, one of the key metal-chelating residues. Additionally, it reduces the flexibility of the gating residues, V304, and E309 in TM4, at the entrance of the Ca²⁺ binding sites, which would decrease the affinity for Ca²⁺. Unexpectedly, SLN has also an effect on the ATP binding site more than 35 Å away, due to the straightening of TM5, a long helix considered as the spine of the protein. The straightening of TM5 modifies the structure of the P-N linker that sits above it, and which comprises the ³⁵¹DKTG³⁵⁴ conserved motif, resulting in an increase of the distance between ATP and the phosphorylation site. As a consequence, the turn-over rate could be affected. All this gives SERCA1a the propensity to go toward a Ca²⁺ low-affinity E2-like state in the presence of SLN and toward a Ca²⁺ high-affinity E1-like state in the absence of SLN. In addition to a general mechanism of inhibition of SERCA1a regulatory peptides, this study also provides an insight into the conformational transition between the E2 and E1 states.

The allosteric activation mechanism of a phospholipase A2-like toxin from Bothrops jararacussu venom: a dynamic description

The activation process of phospholipase A₂-like (PLA₂-like) toxins is a key step in their molecular mechanism, which involves oligomeric changes leading to the exposure of specific sites. Few studies have focused on the characterization of allosteric activators and the features that distinguish them from inhibitors. Herein, a comprehensive study with the BthTX-I toxin from Bothrops jararacussu venom bound or unbound to α-tocopherol (αT) was carried out. The oligomerization state of BthTX-I bound or unbound to αT in solution was studied and indicated that the toxin is predominantly monomeric but tends to oligomerize when complexed with αT. In silico molecular simulations showed the toxin presents higher conformational changes in the absence of αT, which suggests that it is important to stabilize the structure of the toxin. The transition between the two states (active/inactive) was also studied, showing that only the unbound BthTX-I system could migrate to the inactive state. In contrast, the presence of αT induces the toxin to leave the inactive state, guiding it towards the active state, with more regions exposed to the solvent, particularly its active site. Finally, the structural determinants necessary for a molecule to be an inhibitor or activator were analyzed in light of the obtained results.

Unveiling functional motions based on point mutations in biased signaling systems: A normal mode study on nerve growth factor bound to TrkA

Many receptors elicit signal transduction by activating multiple intracellular pathways. This transduction can be triggered by a non-specific ligand, which simultaneously activates all the signaling pathways of the receptors. However, the binding of one biased ligand preferentially trigger one pathway over another, in a process called biased signaling. The identification the functional motions related to each of these distinct pathways has a direct impact on the development of new effective and specific drugs. We show here how to detect specific functional motions by considering the case of the NGF/TrkA-Ig2 complex. NGF-mediated TrkA receptor activation is dependent on specific structural motions that trigger the neuronal growth, development, and survival of neurons in nervous system. The R221W mutation in the ngf gene impairs nociceptive signaling. We discuss how the large-scale structural effects of this mutation lead to the suppression of collective motions necessary to induce TrkA activation of nociceptive signaling. Our results suggest that subtle changes in the NGF interaction network due to the point mutation are sufficient to inhibit the motions of TrkA receptors putatively linked to nociception. The methodological approach presented in this article, based jointly on the normal mode analysis and the experimentally observed functional alterations due to point mutations provides an essential tool to reveal the structural changes and motions linked to the disease, which in turn could be necessary for a drug design study.

New Structural insights into Kir channel gating from molecular simulations, HDX-MS and functional studies

Inward rectifier potassium (Kir) channels play diverse and important roles in shaping action potentials in biological membranes. An increasing number of diseases are now known to be directly associated with abnormal Kir function. However, the gating of Kir still remains unknown. To increase our understanding of its gating mechanism, a dynamical view of the entire channel is essential. Here the gating activation was studied using a recent developped in silico method, MDeNM, which combines normal mode analysis and molecular dynamics simulations that showed for the very first time the importance of interrelated collective and localized conformational movements. In particular, we highlighted the role played by concerted movements of the different regions throughout the entire protein, such as the cytoplasmic and transmembrane domains and the slide helices. In addition, the HDX-MS analysis achieved in these studies provided a comprehensive and detailed view of the dynamics associated with open/closed transition of the Kir channel in coherence with the theoretical results. MDeNM gives access to the probability of the different opening states that are in agreement with our electrophysiological experiments. The investigations presented in this article are important to remedy dysfunctional channels and are of interest for designing new pharmacological compounds.

Microscopic calculation of absorption spectra of macromolecules: An analytic approach

We develop a cross-disciplinary approach to analytically compute optical response functions of open macromolecular systems by exploiting the mathematical formalism of quantum field theory (QFT). Indeed, the entries of the density matrix for the electronic excitations interacting with their open dissipative environment are mapped into vacuum-to-vacuum Green's functions in a fictitious relativistic closed quantum system. We show that by re-summing appropriate self-energy diagrams in this dual QFT, it is possible to obtain analytic expressions for the response functions in Mukamel's theory. This yields physical insight into the structure and dynamics of vibronic resonances, since their frequency and width is related to fundamental physical constants and microscopic model parameters. For illustration, we apply this scheme to compute the linear absorption spectrum of the Fenna-Matthews-Olson light harvesting complex, comparing analytic calculations, numerical Monte Carlo simulations, and experimental data.

Molecular, Biological and Structural Features of VL CDR-1 Rb44 Peptide, Which Targets the Microtubule Network in Melanoma Cells

Microtubules are important drug targets in tumor cells, owing to their role in supporting and determining the cell shape, organelle movement and cell division. The complementarity-determining regions (CDRs) of immunoglobulins have been reported to be a source of anti-tumor peptide sequences, independently of the original antibody specificity for a given antigen. We found that, the anti-Lewis B mAb light-chain CDR1 synthetic peptide Rb44, interacted with microtubules and induced depolymerization, with subsequent degradation of actin filaments, leading to depolarization of mitochondrial membrane-potential, increase of ROS, cell cycle arrest at G2/M, cleavage of caspase-9, caspase-3 and PARP, upregulation of Bax and downregulation of Bcl-2, altogether resulting in intrinsic apoptosis of melanoma cells. The in vitro inhibition of angiogenesis was also an Rb44 effect. Peritumoral injection of Rb44L1 delayed growth of subcutaneously grafted melanoma cells in a syngeneic mouse model. L1-CDRs from immunoglobulins and their interactions with tubulin-dimers were explored to interpret effects on microtubule stability. The opening motion of tubulin monomers allowed for efficient L1-CDR docking, impairment of dimer formation and microtubule dissociation. We conclude that Rb44 V_L-CDR1 is a novel peptide that acts on melanoma microtubule network causing cell apoptosis in vitro and melanoma growth inhibition in vivo.

Combining multi-scale modelling methods to decipher molecular motions of a branching sucrase from glycoside-hydrolase family 70

Among α-transglucosylases from Glycoside-Hydrolase family 70, the ΔN₁₂₃-GB-CD2 enzyme derived from the bifunctional DSR-E from L. citreum NRRL B-1299 is particularly interesting as it was the first described engineered Branching Sucrase, not able to elongate glucan polymers from sucrose substrate. The previously reported overall structural organization of this multi-domain enzyme is an intricate U-shape fold conserved among GH70 enzymes which showed a certain conformational variability of the so-called domain V, assumed to play a role in the control of product structures, in available X-ray structures. Understanding the role of functional dynamics on enzyme reaction and substrate recognition is of utmost interest although it remains a challenge for biophysical methods. By combining long molecular dynamics simulation (1μs) and multiple analyses (NMA, PCA, Morelet Continuous Wavelet Transform and Cross Correlations Dynamics), we investigated here the dynamics of ΔN₁₂₃-GB-CD2 alone and in interaction with sucrose substrate. Overall, our results provide the detailed picture at atomic level of the hierarchy of motions occurring along different timescales and how they are correlated, in agreement with experimental structural data. In particular, detailed analysis of the different structural domains revealed cooperative dynamic behaviors such as twisting, bending and wobbling through anti- and correlated motions, and also two structural hinge regions, of which one was unreported. Several highly flexible loops surrounding the catalytic pocket were also highlighted, suggesting a potential role in the acceptor promiscuity of ΔN123-GBD-CD2. Normal modes and essential dynamics underlined an interesting two-fold dynamic of the catalytic domain A, pivoting about an axis splitting the catalytic gorge in two parts. The comparison of the conformational free energy landscapes using principal component analysis of the enzyme in absence or in presence of sucrose, also revealed a more harmonic basin when sucrose is bound with a shift population of the bending mode, consistent with the substrate binding event.

Jumping between protein conformers using normal modes

The relationship between the normal modes of a protein and its functional conformational change has been studied for decades. However, using this relationship in a predictive context remains a challenge. In this work, we demonstrate that, starting from a given protein conformer, it is possible to generate in a single step model conformers that are less than 1 Å (C_α -RMSD) from the conformer which is the known endpoint of the conformational change, particularly when the conformational change is collective in nature. Such accurate model conformers can be generated by following either the so-called robust or the 50 lowest-frequency modes obtained with various Elastic Network Models (ENMs). Interestingly, the quality of many of these models compares well with actual crystal structures, as assessed by the ROSETTA scoring function and PROCHECK. The most accurate and best quality conformers obtained in the present study were generated by using the 50 lowest-frequency modes of an all-atom ENM. However, with less than ten robust modes, which are identified without any prior knowledge of the nature of the conformational change, nearly 90% of the motion described by the 50 lowest-frequency modes of a protein can be captured. Such results strongly suggest that exploring the robust modes of ENMs may prove efficient for sampling the functionally relevant conformational repertoire of many proteins. © 2017 Wiley Periodicals, Inc.

Sampling of conformational ensemble for virtual screening using molecular dynamics simulations and normal mode analysis

Aim: Molecular dynamics simulations and normal mode analysis are well-established approaches to generate receptor conformational ensembles (RCEs) for ligand docking and virtual screening. Here, we report new fast molecular dynamics-based and normal mode analysis-based protocols combined with conformational pocket classifications to efficiently generate RCEs. Materials & Methods: We assessed our protocols on two well-characterized protein targets showing local active site flexibility, dihydrofolate reductase and large collective movements, CDK2. The performance of the RCEs was validated by distinguishing known ligands of dihydrofolate reductase and CDK2 among a dataset of diverse chemical decoys. Results & discussion: Our results show that different simulation protocols can be efficient for generation of RCEs depending on different kind of protein flexibility.[Formula: see text]

Conformational Equilibrium of CDK/Cyclin Complexes by Molecular Dynamics with Excited Normal Modes

Cyclin-dependent kinases (CDKs) and their associated regulatory cyclins are central for timely regulation of cell-cycle progression. They constitute attractive pharmacological targets for development of anticancer therapeutics, since they are frequently deregulated in human cancers and contribute to sustained, uncontrolled tumor proliferation. Characterization of their structural/dynamic features is essential to gain in-depth insight into structure-activity relationships. In addition, the identification of druggable pockets or key intermediate conformations yields potential targets for the development of novel classes of inhibitors. Structural studies of CDK2/cyclin A have provided a wealth of information concerning monomeric/heterodimeric forms of this kinase. There is, however, much less structural information for other CDK/cyclin complexes, including CDK4/cyclin D1, which displays an alternative (open) position of the cyclin partner relative to CDK, contrasting with the closed CDK2/cyclin A conformation. In this study, we carried out normal-mode analysis and enhanced sampling simulations with our recently developed method, molecular dynamics with excited normal modes, to understand the conformational equilibrium on these complexes. Interestingly, the lowest-frequency normal mode computed for each complex described the transition between the open and closed conformations. Exploration of these motions with an explicit-solvent representation using molecular dynamics with excited normal modes confirmed that the closed conformation is the most stable for the CDK2/cyclin A complex, in agreement with their experimentally available structures. On the other hand, we clearly show that an open↔closed equilibrium may exist in CDK4/cyclin D1, with closed conformations resembling that captured for CDK2/cyclin A. Such conformational preferences may result from the distinct distributions of frustrated contacts in each complex. Using the same approach, the putative roles of the Thr(160) phosphoryl group and the T-loop conformation were investigated. These results provide a dynamic view of CDKs revealing intermediate conformations not yet characterized for CDK members other than CDK2, which will be useful for the design of inhibitors targeting critical conformational transitions.

A robust all-atom model for LCAT generated by homology modeling

LCAT is activated by apoA-I to form cholesteryl ester. We combined two structures, phospholipase A2 (PLA2) that hydrolyzes the ester bond at the sn-2 position of oxidized (short) acyl chains of phospholipid, and bacteriophage tubulin PhuZ, as C- and N-terminal templates, respectively, to create a novel homology model for human LCAT. The juxtaposition of multiple structural motifs matching experimental data is compelling evidence for the general correctness of many features of the model: i) The N-terminal 10 residues of the model, required for LCAT activity, extend the hydrophobic binding trough for the sn-2 chain 15-20 Å relative to PLA2. ii) The topography of the trough places the ester bond of the sn-2 chain less than 5 Å from the hydroxyl of the catalytic nucleophile, S181. iii) A β-hairpin resembling a lipase lid separates S181 from solvent. iv) S181 interacts with three functionally critical residues: E149, that regulates sn-2 chain specificity, and K128 and R147, whose mutations cause LCAT deficiency. Because the model provides a novel explanation for the complicated thermodynamic problem of the transfer of hydrophobic substrates from HDL to the catalytic triad of LCAT, it is an important step toward understanding the antiatherogenic role of HDL in reverse cholesterol transport.

On the relationship between low-frequency normal modes and the large-scale conformational changes of proteins

Normal mode analysis is a computational technique that allows to study the dynamics of biological macromolecules. It was first applied to small protein cases, more than thirty years ago. The interest in this technique then raised when it was realized that it can provide insights about the large-scale conformational changes a protein can experience, for instance upon ligand binding. As it was also realized that studying highly simplified protein models can provide similar insights, meaning that this kind of analysis can be both quick and simple to handle, several applications were proposed, in the context of various structural biology techniques. This review focuses on these applications, as well as on how the functional relevance of the lowest-frequency modes of proteins was established.

Differential effects of CSF-1R D802V and KIT D816V homologous mutations on receptor tertiary structure and allosteric communication

The colony stimulating factor-1 receptor (CSF-1R) and the stem cell factor receptor KIT, type III receptor tyrosine kinases (RTKs), are important mediators of signal transduction. The normal functions of these receptors can be compromised by gain-of-function mutations associated with different physiopatological impacts. Whereas KIT D816V/H mutation is a well-characterized oncogenic event and principal cause of systemic mastocytosis, the homologous CSF-1R D802V has not been identified in human cancers. The KIT D816V oncogenic mutation triggers resistance to the RTK inhibitor Imatinib used as first line treatment against chronic myeloid leukemia and gastrointestinal tumors. CSF-1R is also sensitive to Imatinib and this sensitivity is altered by mutation D802V. Previous in silico characterization of the D816V mutation in KIT evidenced that the mutation caused a structure reorganization of the juxtamembrane region (JMR) and facilitated its departure from the kinase domain (KD). In this study, we showed that the equivalent CSF-1R D802V mutation does not promote such structural effects on the JMR despite of a reduction on some key H-bonds interactions controlling the JMR binding to the KD. In addition, this mutation disrupts the allosteric communication between two essential regulatory fragments of the receptors, the JMR and the A-loop. Nevertheless, the mutation-induced shift towards an active conformation observed in KIT D816V is not observed in CSF-1R D802V. The distinct impact of equivalent mutation in two homologous RTKs could be associated with the sequence difference between both receptors in the native form, particularly in the JMR region. A local mutation-induced perturbation on the A-loop structure observed in both receptors indicates the stabilization of an inactive non-inhibited form, which Imatinib cannot bind.

Aligning experimental and theoretical anisotropic B-factors: water models, normal-mode analysis methods, and metrics

The strength of X-ray crystallography in providing the information for protein dynamics has been under appreciated. The anisotropic B-factors (ADPs) from high-resolution structures are invaluable in studying the relationship among structure, dynamics, and function. Here, starting from an in-depth evaluation of the metrics used for comparing the overlap between two ellipsoids, we applied normal-mode analysis (NMA) to predict the theoretical ADPs and then align them with experimental results. Adding an extra layer of explicitly treated water on protein surface significantly improved the energy minimization results and better reproduced the anisotropy of experimental ADPs. In comparing experimental and theoretical ADPs, we focused on the overlap in shape, the alignment of dominant directions, and the similarity in magnitude. The choices of water molecules, NMA methods, and the metrics for evaluating the overlap of ADPs determined final results. This study provides useful information for exploring the physical basis and the application potential of experimental ADPs.

An atomistic view of microtubule stabilization by GTP

A microtubule is a dynamic system formed of αβ-tubulins. The presence of nonhydrolyzable guanosine-5'-triphosphate (GTP)/guanosine diphosphate (GDP) on the β-tubulins provokes microtubule polymerization/depolymerization. Despite the large number of experimental studies of this dynamical process, its mechanism is still unclear. To provide insights into this mechanism we studied the first depolymerization steps of GDP/GTP-bound microtubules by normal-mode analysis with the all-atom model. We also constructed a depolymerizing microtubule and compared it to cryo-electron microscopy tomograms (cyro-ET). The results show that during depolymerization, the protofilaments not only curve but twist to weaken their lateral interactions. These interactions are stabilized by GTP, but not evenly. Not all of the interface residues are of equal importance: five of them, belonging to the H2-S3 loop, play a special role; acting as a lock whose key is the γ-phosphate of GTP. Sequence alignments of several tubulins confirm the importance of these residues.

Elastic network models: theoretical and empirical foundations

Fifteen years ago, Monique Tirion showed that the low-frequency normal modes of a protein are not significantly altered when nonbonded interactions are replaced by Hookean springs, for all atom pairs whose distance is smaller than a given cutoff value. Since then, it has been shown that coarse-grained versions of Tirion's model are able to provide fair insights on many dynamical properties of biological macromolecules. In this chapter, theoretical tools required for studying these so-called Elastic Network Models are described, focusing on practical issues and, in particular, on possible artifacts. Then, an overview of some typical results that have been obtained by studying such models is given.

Elastic network normal modes provide a basis for protein structure refinement

It is well recognized that thermal motions of atoms in the protein native state, the fluctuations about the minimum of the global free energy, are well reproduced by the simple elastic network models (ENMs) such as the anisotropic network model (ANM). Elastic network models represent protein dynamics as vibrations of a network of nodes (usually represented by positions of the heavy atoms or by the C(α) atoms only for coarse-grained representations) in which the spatially close nodes are connected by harmonic springs. These models provide a reliable representation of the fluctuational dynamics of proteins and RNA, and explain various conformational changes in protein structures including those important for ligand binding. In the present paper, we study the problem of protein structure refinement by analyzing thermal motions of proteins in non-native states. We represent the conformational space close to the native state by a set of decoys generated by the I-TASSER protein structure prediction server utilizing template-free modeling. The protein substates are selected by hierarchical structure clustering. The main finding is that thermal motions for some substates, overlap significantly with the deformations necessary to reach the native state. Additionally, more mobile residues yield higher overlaps with the required deformations than do the less mobile ones. These findings suggest that structural refinement of poorly resolved protein models can be significantly enhanced by reduction of the conformational space to the motions imposed by the dominant normal modes.

Mutation D816V alters the internal structure and dynamics of c-KIT receptor cytoplasmic region: implications for dimerization and activation mechanisms

The type III receptor tyrosine kinase (RTK) KIT plays a crucial role in the transmission of cellular signals through phosphorylation events that are associated with a switching of the protein conformation between inactive and active states. D816V KIT mutation is associated with various pathologies including mastocytosis and cancers. D816V-mutated KIT is constitutively active, and resistant to treatment with the anti-cancer drug Imatinib. To elucidate the activating molecular mechanism of this mutation, we applied a multi-approach procedure combining molecular dynamics (MD) simulations, normal modes analysis (NMA) and binding site prediction. Multiple 50-ns MD simulations of wild-type KIT and its mutant D816V were recorded using the inactive auto-inhibited structure of the protein, characteristic of type III RTKs. Computed free energy differences enabled us to quantify the impact of D816V on protein stability in the inactive state. We evidenced a local structural alteration of the activation loop (A-loop) upon mutation, and a long-range structural re-organization of the juxta-membrane region (JMR) followed by a weakening of the interaction network with the kinase domain. A thorough normal mode analysis of several MD conformations led to a plausible molecular rationale to propose that JMR is able to depart its auto-inhibitory position more easily in the mutant than in wild-type KIT and is thus able to promote kinase mutant dimerization without the need for extra-cellular ligand binding. Pocket detection at the surface of NMA-displaced conformations finally revealed that detachment of JMR from the kinase domain in the mutant was sufficient to open an access to the catalytic and substrate binding sites.

Protein kinases are involved in a huge amount of cellular processes through phosphorylation, a crucial mechanism in cell signaling, and their misregulation often results in disease. The deactivation of protein tyrosine kinases (PTKs) or their oncogenic activation arises from mutations which affect the protein primary structure and the configuration of the enzymatic site apparently by stabilizing the activation loop (A-loop) extended conformation. Particularly, mutation D816V of receptor tyrosine kinase (RTK) KIT, found in patients with pediatric mastocytosis, acute leukemia or germ cell tumors, can be considered as the archetype of mutation inducing a displacement of the population equilibrium toward the active conformation. We present a comprehensive computational study of the activating mechanism(s) of this mutation. Our multi-approach in silico procedure evidenced a local alteration of the A-loop structure, and a long-range structural re-organization of the juxta-membrane region (JMR) followed by a weakening of the interaction network with the kinase domain. Our results provided a plausible conception of how the observed departure of JMR from kinase domain in the mutant promotes kinase mutant dimerization without requiring extra-cellular ligand binding. The pocket profiles we obtained suggested putative allosteric binding sites that could be targeted by ligands/modulators that trap the mutated enzyme.

A concerted mechanism for opening the GDP binding pocket and release of the nucleotide in hetero-trimeric G-proteins

G-protein hetero-trimers play a fundamental role in cell function. Their dynamic behavior at the atomic level remains to be understood. We have studied the Gi hetero-trimer through a combination of molecular dynamics simulations and normal mode analyses. We showed that these big proteins could undergo large-amplitude conformational changes, without any energy penalty and with an intrinsic dynamics centered on their GDP binding pocket. Among the computed collective motions, one of the modes (mode 17) was particularly able to significantly open both the base and the phosphate sides of the GDP binding pocket. This mode describing mainly a motion between the Ras-like and the helical domains of G(α) was in close agreement with some available X-ray data and with many other biochemical/biophysical observations including the kink of helix α5. The use of a new protocol, which allows extraction of the GDP ligand along the computed normal modes, supported that the exit of GDP was largely coupled to an opening motion along mode 17. We propose for the first time a "concerted mechanism" model in which the opening of the GDP pocket and the kink of the α5 helix occur concomitantly and favor GDP release from G(αβγ) complexes. This model is discussed in the context of the G-protein-coupled receptor/G-protein interaction close to the cell membrane.

Long-range energy transfer in proteins

Proteins are large and complex molecular machines. In order to perform their function, most of them need energy, e.g. either in the form of a photon, as in the case of the visual pigment rhodopsin, or through the breaking of a chemical bond, as in the presence of adenosine triphosphate (ATP). Such energy, in turn, has to be transmitted to specific locations, often several tens of A away from where it is initially released. Here we show, within the framework of a coarse-grained nonlinear network model, that energy in a protein can jump from site to site with high yields, covering in many instances remarkably large distances. Following single-site excitations, few specific sites are targeted, systematically within the stiffest regions. Such energy transfers mark the spontaneous formation of a localized mode of nonlinear origin at the destination site, which acts as an efficient energy-accumulating center. Interestingly, yields are found to be optimum for excitation energies in the range of biologically relevant ones.

Activation of the ghrelin receptor is described by a privileged collective motion: a model for constitutive and agonist-induced activation of a sub-class A G-protein coupled receptor (GPCR)

Three homology models of the human ghrelin receptor (GHS-R1a) have been generated from the available X-ray structures of rhodopsin (RHO model), opsin (OPS model) and beta-2 adrenergic receptor (B2 model). The latter was used as a starting point for combined molecular dynamics simulation (MDS) and full atom normal modes analysis (NMA). A low-frequency normal mode (mode 16) perfectly reproduced the intracellular motions observed between B2 and RHO models; in the opposite direction along the same mode, the generated structures are closer to the OPS model, suggesting a direct link with GHS-R1a activation. This was in agreement with motions of the seven transmembranous segments, increase of the solvent accessibility of the 140-ERY-142 sequence, and flip of the Trp276 (C WLP) residue, some features related to GPCRs activation. According to our model, His280 was proposed to stabilize Trp276 in the active state; this was verified by site-directed mutagenesis and biochemical characterization of the resulting H280A and H280S mutants, which were fully functional but sharing an important decrease of their basal activities. Docking performed with short ghrelin derivatives Gly-Ser-Ser ([octa])-Phe-NH (2) and Gly-Ser-Ser ([octa])-Phe-Leu-NH (2) allowed the identification of a robust position of these peptides in the active site of the receptor. This model was refined by MDS and validated by docking experiments performed on a set of 55 ghrelin receptor ligands based on the 1,2,4- triazole scaffold. Finally, NMA performed on the obtained peptide-receptor complex suggested stabilization of the Trp276 residue and of the whole receptor in the active state, preventing the motion observed along mode 16 computed for the unbound receptor. Our results show that NMA offers a powerful approach to study the conformational diversity and the activation mechanism of GPCRs.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Determination of phosphorylation sites for NADP-specific isocitrate dehydrogenase from mycobacterium tuberculosis

Isocitrate Dehydrogenase (ICD) catalyzes the oxidative decarboxylation reaction of 2R,3S-isocitrate to yield 2-oxoglutarate in the Tricarboxylic Acid (TCA) cycle. Two isoforms of NADP-specific ICDs with the E.C number 1.1.1.42 have been annotated in the organism Mycobacterium tuberculosis, monomeric ICD2 and dimeric ICD1. BLAST search against the Protein Data Bank (PDB) database shows a marked similarity between dimeric Mycobacterium tuberculosis ICD1 sequence and that of Sus scrofa, a cytosolic eukaryotic ICD (65% identity). Escherischia coli ICD shows less sequence similarity than the eukaryotic structure. A Homology model has thus been built for M. tuberculosis ICD1 using Sus scrofa and human ICD as templates. Inactivation of ICD1 by phosphorylation similar to E. coli ICD is important to open up the shunt pathway in the TCA cycle, which has been indicated in the case of M. tuberculosis. We therefore attempted to identify a number of likely phosphorylation sites in M. tuberculosis using pattern prediction and checked with the homology models for the accessibility of the peptides containing Serine. It was found that the homologous Serine by alignment with E. coli on M. tuberculosis ICD1 is difficult to access by specific kinases. Hence other probable sites of phosphorylation were checked and three highly probable serine-containing peptides were identified. The effect of phosphorylation at each of these sites was determined by checking the degree of conformational changes, the differences caused by the effect of phosphorylation in the active-site and other apparent motion different from that of the control, i.e., unphosphorylated M. tuberculosis ICD1 model, using molecular dynamics simulations.

Collective motions in glucosamine-6-phosphate synthase: influence of ligand binding and role in ammonia channelling and opening of the fructose-6-phosphate binding site

The large protein motions of the bacterial enzyme glucosamine-6-phosphate synthase have been addressed using full atom normal modes analysis for the empty, the glucose-6-phosphate and the glucose-6-phosphate+glutamate bound proteins. The approach that was used involving energy minimizations along the normal modes coordinates identified functional motions of the protein, some of which were characterized earlier by X-ray diffraction studies. This method made it possible for the first time to highlight significant energy differences according to whether none, only one or both of the active sites of the protein were occupied. Our data favoured a specific motion of the glutamine binding domain following the fixation of fructose-6-phosphate and suggested a rigidified structure with both sites occupied. Here, we show that most of the collective large amplitude motions of glucosamine-6-phosphate synthase that are modulated by ligand binding are crucial for the enzyme catalytic cycle, as they strongly modify the geometry of both the ammonia channel and the C-tail, demonstrating their role in ammonia transfer and ligand binding.

Discrete breathers in protein structures

Recently, using a numerical surface cooling approach, we have shown that highly energetic discrete breathers (DBs) can form in the stiffest parts of nonlinear network models of large protein structures. In the present study, using an analytical approach, we extend our previous results to low-energy discrete breathers as well as to smaller proteins. We confirm and further scrutinize the striking site selectiveness of energy localization in the presence of spatial disorder. In particular, we find that, as a sheer consequence of disorder, a non-zero energy gap for exciting a DB at a given site either exists or not. Remarkably, in the former case, the gaps arise as a result of the impossibility of exciting small-amplitude modes in the first place. In contrast, in the latter case, a small subset of linear edge modes acts as accumulation points, whereby DBs can be continued to arbitrary small energies, while unavoidably approaching one of such normal modes. In particular, the case of the edge mode seems peculiar, its dispersion relation being simple and little system dependent. Concerning the structure-dynamics relationship, we find that the regions of protein structures where DBs form easily (zero or small gaps) are unfailingly the most highly connected ones, also characterized by weak local clustering. Remarkably, a systematic analysis on a large database of enzyme structures reveals that amino-acid residues involved in catalysis tend to be located in such regions. This finding reinforces the idea that localized modes of nonlinear origin may play an important biological role, e.g., by providing a ready channel for energy storage and/or contributing to lower energy barriers of chemical reactions.

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.

Discrete breathers in nonlinear network models of proteins

We introduce a topology-based nonlinear network model of protein dynamics with the aim of investigating the interplay of spatial disorder and nonlinearity. We show that spontaneous localization of energy occurs generically and is a site-dependent process. Localized modes of nonlinear origin form spontaneously in the stiffest parts of the structure and display site-dependent activation energies. Our results provide a straightforward way for understanding the recently discovered link between protein local stiffness and enzymatic activity. They strongly suggest that nonlinear phenomena may play an important role in enzyme function, allowing for energy storage during the catalytic process.

Domain motions of Argonaute, the catalytic engine of RNA interference

Background

The Argonaute protein is the core component of the RNA-induced silencing complex, playing the central role of cleaving the mRNA target. Visual inspection of static crystal structures already has enabled researchers to suggest conformational changes of Argonaute that might occur during RNA interference. We have taken the next step by performing an all-atom normal mode analysis of the Pyrococcus furiosus and Aquifex aeolicus Argonaute crystal structures, allowing us to quantitatively assess the feasibility of these conformational changes. To perform the analysis, we begin with the energy-minimized X-ray structures. Normal modes are then calculated using an all-atom molecular mechanics force field.

Results

The analysis reveals low-frequency vibrations that facilitate the accommodation of RNA duplexes – an essential step in target recognition. The Pyrococcus furiosus and Aquifex aeolicus Argonaute proteins both exhibit low-frequency torsion and hinge motions; however, differences in the overall architecture of the proteins cause the detailed dynamics to be significantly different.

Conclusion

Overall, low-frequency vibrations of Argonaute are consistent with mechanisms within the current reaction cycle model for RNA interference.

Normal-mode flexible fitting of high-resolution structure of biological molecules toward one-dimensional low-resolution data

We present a method for reconstructing a 3D structure from a pair distribution function by flexibly fitting known x-ray structures toward a conformation that agrees with the low-resolution data. This method uses a linear combination of low-frequency normal modes from elastic-network description of the molecule in an iterative manner to deform the structure optimally to conform to the target pair distribution function. A simple function, pair distance distribution function between atoms, is chosen as a test model to establish computational algorithms-optimization algorithm and scoring function-that can utilize low-resolution 1D data. To select a correct structural model based on less information, we developed a scoring function that takes into account a characteristic of pair distribution functions. In addition, we employ a new optimization algorithm, the trusted region method, that relies on both first and second derivatives of the scoring function. Illustrative results of our studies on simulated 1D data from five different proteins, for which large conformational changes are known to occur, are presented.

Coarse-grained biomolecular simulation with REACH: realistic extension algorithm via covariance Hessian

Coarse-graining of protein interactions provides a means of simulating large biological systems. Here, a coarse-graining method, REACH, is introduced, 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. In test calculations, the C(alpha)-atoms variance-covariance matrices are calculated from the ensembles of 1-ns atomistic MD trajectories in monomeric and dimeric myoglobin, and used to derive coarse-grained force constants for the local and nonbonded interactions. Construction of analytical model functions of the distance-dependence of the interresidue force constants allows rapid calculation of the REACH normal modes. The model force constants from monomeric and dimeric myoglobin are found to be similar in magnitude to each other. The MD intra- and intermolecular mean-square fluctuations and the vibrational density of states are well reproduced by the residue-scale REACH normal modes without requiring rescaling of the force constant parameters. The temperature-dependence of the myoglobin REACH force constants reveals that the dynamical transition in protein internal fluctuations arises principally from softening of the elasticity in the nonlocal interactions. The REACH method is found to be a reliable way of determining spatiotemporal protein motion without the need for expensive computations of long atomistic MD simulations.

Single-atom test of all-atom empirical potentials: Fe in myoglobin

The measured Fe vibrational density of states in deoxy-myoglobin, obtained from nuclear resonance vibrational spectroscopy, is compared to results from a normal-mode analysis using an all-atom empirical potential. Substantial disagreement reveals that for this one atom, the empirical potential does not accurately describe the actual forces. A Green function technique is developed to calculate the iron vibrational spectrum of deoxy-myoglobin by coupling the independently calculated heme and globin normal modes; nonbonded interactions between the heme molecule and the protein are essential for a good fit to the measurements. A projection of the eigenvectors from this potential onto the displacements induced by binding of CO demonstrates that normal modes over a broad range centered around 50-150 cm(-1) may drive the ligand-induced structural changes.

How well can we understand large-scale protein motions using normal modes of elastic network models?

In this article, we apply a coarse-grained elastic network model (ENM) to study conformational transitions to address the following questions: How well can a conformational change be predicted by the mode motions? Is there a way to improve the model to gain better results? To answer these questions, we use a dataset of 170 pairs having "open" and "closed" structures from Gerstein's protein motion database. Our results show that the conformational transitions fall into three categories: 1), the transitions of these proteins that can be explained well by ENM; 2), the transitions that are not explained well by ENM, but the results are significantly improved after considering the rigidity of some residue clusters and modeling them accordingly; and 3), the intrinsic nature of these transitions, specifically the low degree of collectivity, prevents their conformational changes from being represented well with the low frequency modes of any elastic network models. Our results thus indicate that the applicability of ENM for explaining conformational changes is not limited by the size of the studied protein or even the scale of the conformational change. Instead, it depends strongly on how collective the transition is.

Normal mode analysis as a prerequisite for drug design: Application to matrix metalloproteinases inhibitors

We demonstrate the utility of normal mode analysis in correctly predicting the binding modes of inhibitors in the active sites of matrix metalloproteinases (MMPs). We show the accuracy in predicting the positions of MMP-3 inhibitors is strongly dependent on which structure is used as the target, especially when it has been energy minimized. This dependency can be overcome by using intermediate structures generated along one of the normal modes previously calculated for a given target. These results may be of prime importance for further in silico drug discovery.