Binding free energies and free energy components from molecular dynamics and Poisson-Boltzmann calculations. Application to amino acid recognition by aspartyl-tRNA synthetase

Sucrose and Gibberellic Acid Binding Stabilize the Inward-Open Conformation of AtSWEET13: A Molecular Dynamics Study

In plants, sugar will eventually be exported transporters (SWEETs) facilitate the translocation of mono- and disaccharides across membranes and play a critical role in modulating responses to gibberellin (GA3), a key growth hormone. However, the dynamic mechanisms underlying sucrose and GA3 binding and transport remain elusive. Here, we employed microsecond-scale molecular dynamics (MD) simulations to investigate the influence of sucrose and GA3 binding on SWEET13 transporter motions. While sucrose exhibits high flexibility within the binding pocket, GA3 remains firmly anchored in the central cavity. Binding of both ligands increases the average channel radius along the transporter's principal axis. In contrast to the apo form, which retains its initial conformation throughout the simulation, ligand-bound complexes undergo a significant conformational transition characterized by further opening of the intracellular gate relative to the inward-open crystal structure (5XPD). This opening is driven by ligand-induced bending of helix V, stabilizing the inward-open state. Sucrose binding notably enhances the flexibility of the intracellular gate and amplifies anticorrelated motions between the N- and C-terminal domains at the intracellular side, suggesting an opening-closing motion of these domains. Principal component analysis revealed that this gating motion is most pronounced in the sucrose complex and minimal in the apo form, highlighting sucrose's ability to induce high-amplitude gating. Our binding free energy calculations indicate that SWEET13 has lower binding affinity for sucrose compared to GA3, consistent with its role in sugar transport. These results provide insight into key residues involved in sucrose and GA3 binding and transport, advancing our understanding of SWEET13 dynamics.

Uncovering substrate specificity determinants of class IIb aminoacyl-tRNA synthetases with machine learning

Specific amino acid (AA) binding by aminoacyl-tRNA synthetases (aaRSs) is necessary for correct translation of the genetic code. Sequence and structure analyses have revealed the main specificity determinants and allowed a partitioning of aaRSs into two classes and several subclasses. However, the information contributed by each determinant has not been precisely quantified, and other, minor determinants may still be unidentified. Growth of genomic data and development of machine learning classification methods allow us to revisit these questions. This work considered the subclass IIb, formed by the three enzymes aspartyl-, asparaginyl-, and lysyl-tRNA synthetase (LysRS). Over 35,000 sequences from the Pfam database were considered, and used to train a machine-learning model based on ensembles of decision trees. The model was trained to reproduce the existing classification of each sequence as AspRS, AsnRS, or LysRS, and to identify which sequence positions were most important for the classification. A few positions (5-8 depending on the AA substrate) sufficed for accurate classification. Most but not all of them were well-known specificity determinants. The machine learning models thus identified sets of mutations that distinguish the three subclass members, which might be targeted in engineering efforts to alter or swap the AA specificities for biotechnology applications.

Estrogen receptor modulation of some polyphenols extracted from Daucus carota as a probable mechanism for antifertility effect: An in silico study

The seeds of Daucus carota, traditionally used by women in many countries to prevent conception, were proved to have negative impact on reproductive hormone levels as well as on the estrous cycles in albino mice affecting the fertility status. This study is an attempt to investigate the possible role of polyphenols present in the seeds in hampering the reproductive processes. The Molecular Docking, Molecular Dynamics (MD) simulation and binding free energy calculation studies reveal that six polyphenols present in the seeds can bind with the active sites of human Estrogen Receptor (ER) and may interfere in the estrogen signaling in human. These polyphenols were found to bind to a conservative pocket of ER[Formula: see text], which is comprised of residues 343–388, 421–428 and 525–540. Docking studies indicated the presence of strong hydrogen bonding, pi–pi interactions and numerous hydrophobic interactions that stabilize the ER[Formula: see text]-polyphenol complexes.

The docked complexes were further subjected to MM/GBSA analysis to calculate binding free energies. Molecular dynamic simulation studies carried out for a period of 20[Formula: see text]ns revealed low RMS deviation values suggesting high accuracy of the docking poses and stability of the complexes. Out of the six polyphenols, catechin and epicatechin have shown highest binding affinity towards the ER[Formula: see text] receptor. These findings will help in identifying ER modulators of plant origin targeting ER alpha and predicting their effects on the reproductive hormone homeostasis. Moreover, this study may form preliminary basis for further identification of potential herbal antifertility agents.

Molecular dynamics simulations of cognate and non-cognate AspRS-tRNAAsp complexes

Aspartyl tRNA synthetase (AspRS), one of the 20 aminoacyl-tRNA synthetases, plays an important role in protein synthesis by catalyzing the aminoacylation reaction and synthesises Aspartyl-tRNA (tRNA^Asp). A typical three-dimensional structure of AspRS comprises three distinct domains for the recognition of cognate tRNA and catalysis, namely, anti-codon binding domain/N-terminal domain, hinge domain and catalytic domain through their interactions with anti-codon loop, D-stem and acceptor arm of cognate tRNA, respectively. In this work, we have studied the structural characteristics of each domain of AspRS to understand the recognition mechanism of tRNA^Asp using molecular dynamics simulations. The dynamics of AspRS-tRNA^Asp complexes from E.coli (cognate and non-cognate), S.cerevisiae (cognate) and T.thermophilus (non-cognate) were compared to understand the differences in recognition of cognate and non-cognate tRNAs. Our results explain that the conformational changes associated with the recognition of tRNA occur only in the cognate complexes. Among the cognate complexes, the conformational changes in yeast AspRS are highly controlled during tRNA^Asp recognition than that of in the E. coli AspRS. Moreover, the functional motions required for the tRNA recognition are observed only in the cognate complexes, and the conformational changes in AspRS and their recognition of tRNA^Asp are organism specific.Communicated by Ramaswamy H. Sarma.

Oxocarbon-functionalized graphene as a lithium-ion battery cathode: a first-principles investigation

In recent years, organic-based, especially carbonyl-based, Li-ion battery electrode materials have attracted great attention due to their low-cost, environmentally friendly nature and strong Li-ion bonding abilities.

Electrostatic component of binding energy: Interpreting predictions from poisson-boltzmann equation and modeling protocols

Macromolecular interactions are essential for understanding numerous biological processes and are typically characterized by the binding free energy. Important component of the binding free energy is the electrostatics, which is frequently modeled via the solutions of the Poisson-Boltzmann Equations (PBE). However, numerous works have shown that the electrostatic component (ΔΔGelec ) of binding free energy is very sensitive to the parameters used and modeling protocol. This prompted some researchers to question the robustness of PBE in predicting ΔΔGelec . We argue that the sensitivity of the absolute ΔΔGelec calculated with PBE using different input parameters and definitions does not indicate PBE deficiency, rather this is what should be expected. We show how the apparent sensitivity should be interpreted in terms of the underlying changes in several numerous and physical parameters. We demonstrate that PBE approach is robust within each considered force field (CHARMM-27, AMBER-94, and OPLS-AA) once the corresponding structures are energy minimized. This observation holds despite of using two different molecular surface definitions, pointing again that PBE delivers consistent results within particular force field. The fact that PBE delivered ΔΔGelec values may differ if calculated with different modeling protocols is not a deficiency of PBE, but natural results of the differences of the force field parameters and potential functions for energy minimization. In addition, while the absolute ΔΔGelec values calculated with different force field differ, their ordering remains practically the same allowing for consistent ranking despite of the force field used. © 2016 Wiley Periodicals, Inc.

Direct Calculation of Protein Fitness Landscapes through Computational Protein Design

Naturally selected amino-acid sequences or experimentally derived ones are often the basis for understanding how protein three-dimensional conformation and function are determined by primary structure. Such sequences for a protein family comprise only a small fraction of all possible variants, however, representing the fitness landscape with limited scope. Explicitly sampling and characterizing alternative, unexplored protein sequences would directly identify fundamental reasons for sequence robustness (or variability), and we demonstrate that computational methods offer an efficient mechanism toward this end, on a large scale. The dead-end elimination and A(∗) search algorithms were used here to find all low-energy single mutant variants, and corresponding structures of a G-protein heterotrimer, to measure changes in structural stability and binding interactions to define a protein fitness landscape. We established consistency between these algorithms with known biophysical and evolutionary trends for amino-acid substitutions, and could thus recapitulate known protein side-chain interactions and predict novel ones.

Protein:Ligand binding free energies: A stringent test for computational protein design

A computational protein design method is extended to allow Monte Carlo simulations where two ligands are titrated into a protein binding pocket, yielding binding free energy differences. These provide a stringent test of the physical model, including the energy surface and sidechain rotamer definition. As a test, we consider tyrosyl-tRNA synthetase (TyrRS), which has been extensively redesigned experimentally. We consider its specificity for its substrate l-tyrosine (l-Tyr), compared to the analogs d-Tyr, p-acetyl-, and p-azido-phenylalanine (ac-Phe, az-Phe). We simulate l- and d-Tyr binding to TyrRS and six mutants, and compare the structures and binding free energies to a more rigorous "MD/GBSA" procedure: molecular dynamics with explicit solvent for structures and a Generalized Born + Surface Area model for binding free energies. Next, we consider l-Tyr, ac- and az-Phe binding to six other TyrRS variants. The titration results are sensitive to the precise rotamer definition, which involves a short energy minimization for each sidechain pair to help relax bad contacts induced by the discrete rotamer set. However, when designed mutant structures are rescored with a standard GBSA energy model, results agree well with the more rigorous MD/GBSA. As a third test, we redesign three amino acid positions in the substrate coordination sphere, with either l-Tyr or d-Tyr as the ligand. For two, we obtain good agreement with experiment, recovering the wildtype residue when l-Tyr is the ligand and a d-Tyr specific mutant when d-Tyr is the ligand. For the third, we recover His with either ligand, instead of wildtype Gln.

MD Simulations of tRNA and Aminoacyl-tRNA Synthetases: Dynamics, Folding, Binding, and Allostery

While tRNA and aminoacyl-tRNA synthetases are classes of biomolecules that have been extensively studied for decades, the finer details of how they carry out their fundamental biological functions in protein synthesis remain a challenge. Recent molecular dynamics (MD) simulations are verifying experimental observations and providing new insight that cannot be addressed from experiments alone. Throughout the review, we briefly discuss important historical events to provide a context for how far the field has progressed over the past few decades. We then review the background of tRNA molecules, aminoacyl-tRNA synthetases, and current state of the art MD simulation techniques for those who may be unfamiliar with any of those fields. Recent MD simulations of tRNA dynamics and folding and of aminoacyl-tRNA synthetase dynamics and mechanistic characterizations are discussed. We highlight the recent successes and discuss how important questions can be addressed using current MD simulations techniques. We also outline several natural next steps for computational studies of AARS:tRNA complexes.

Resolving hot spots in the C-terminal dimerization domain that determine the stability of the molecular chaperone Hsp90

Human heat shock protein of 90 kDa (hHsp90) is a homodimer that has an essential role in facilitating malignant transformation at the molecular level. Inhibiting hHsp90 function is a validated approach for treating different types of tumors. Inhibiting the dimerization of hHsp90 via its C-terminal domain (CTD) should provide a novel way to therapeutically interfere with hHsp90 function. Here, we predicted hot spot residues that cluster in the CTD dimerization interface by a structural decomposition of the effective energy of binding computed by the MM-GBSA approach and confirmed these predictions using in silico alanine scanning with DrugScore(PPI). Mutation of these residues to alanine caused a significant decrease in the melting temperature according to differential scanning fluorimetry experiments, indicating a reduced stability of the mutant hHsp90 complexes. Size exclusion chromatography and multi-angle light scattering studies demonstrate that the reduced stability of the mutant hHsp90 correlates with a lower complex stoichiometry due to the disruption of the dimerization interface. These results suggest that the identified hot spot residues can be used as a pharmacophoric template for identifying and designing small-molecule inhibitors of hHsp90 dimerization.

Molecular recognition mechanism of peptide chain bound to the tRNA(Lys3) anticodon loop in silico

The mechanism by which proteins recognize and bind the post-transcriptional modifications of RNAs is unknown, yet these interactions play important functions in biology. Atomistic molecular dynamics simulations were performed to examine the folding of the model peptide chain -RVTHHAFLGAHRTVG- and the complex formed by the folded peptide with the native anticodon stem and loop of the human tRNA(Lys3) (hASL(Lys3)) in order to explore the binding mechanism. By analyzing and comparing two folded conformations of this peptide obtained from the folding simulation, we found that the van der Waals (VDW) energy is necessary for the thermal stability of the peptide, and the charge-charge (ELE + EGB) energy is crucial for determining the three-dimensional folded structure of the peptide backbone. Subsequently, two conformations of the peptide were employed to investigate their binding behaviors to hASL(Lys3). The metastable folded peptide was found to bind to hASL(Lys3) much easier than the stable folded peptide in the binding simulations. An energetic analysis reveals that the VDW energy favors the binding, whereas the ELE + EGB energies disfavor the binding. Arginines on the peptide preferentially attract the phosphate backbone via the inter-chain ELE + EGB interaction, significantly contributing to the binding affinity. The hydrophobic phenylalanine interacts with the anticodon loop of hASL(Lys3) via the inter-chain VDW interaction, significantly contributing to the binding specificity.

Combination of theoretical and experimental approaches for the design and study of fibril-forming peptides

Self-assembling peptides that can form supramolecular structures such as fibrils, ribbons, and nanotubes are of particular interest to modern bionanotechnology and materials science. Their ability to form biocompatible nanostructures under mild conditions through non-covalent interactions offers a big biofabrication advantage. Structural motifs extracted from natural proteins are an important source of inspiration for the rational design of such peptides. Examples include designer self-assembling peptides that correspond to natural coiled-coil motifs, amyloid-forming proteins, and natural fibrous proteins. In this chapter, we focus on the exploitation of structural information from beta-structured natural fibers. We review a case study of short peptides that correspond to sequences from the adenovirus fiber shaft. We describe both theoretical methods for the study of their self-assembly potential and basic experimental protocols for the assessment of fibril-forming assembly.

Computational protein design: the Proteus software and selected applications

We describe an automated procedure for protein design, implemented in a flexible software package, called Proteus. System setup and calculation of an energy matrix are done with the XPLOR modeling program and its sophisticated command language, supporting several force fields and solvent models. A second program provides algorithms to search sequence space. It allows a decomposition of the system into groups, which can be combined in different ways in the energy function, for both positive and negative design. The whole procedure can be controlled by editing 2-4 scripts. Two applications consider the tyrosyl-tRNA synthetase enzyme and its successful redesign to bind both O-methyl-tyrosine and D-tyrosine. For the latter, we present Monte Carlo simulations where the D-tyrosine concentration is gradually increased, displacing L-tyrosine from the binding pocket and yielding the binding free energy difference, in good agreement with experiment. Complete redesign of the Crk SH3 domain is presented. The top 10000 sequences are all assigned to the correct fold by the SUPERFAMILY library of Hidden Markov Models. Finally, we report the acid/base behavior of the SNase protein. Sidechain protonation is treated as a form of mutation; it is then straightforward to perform constant-pH Monte Carlo simulations, which yield good agreement with experiment. Overall, the software can be used for a wide range of application, producing not only native-like sequences but also thermodynamic properties with errors that appear comparable to other current software packages.

Carbohydrate recognition by the antiviral lectin cyanovirin-N

Cyanovirin-N (CVN) is a cyanobacterial lectin with potent antiviral activity and has been the focus of extensive preclinical investigation as a potential prophylactic for the prevention of the sexual transmission of the human immunodeficiency virus (HIV). Here we present a detailed analysis of carbohydrate recognition by this important protein, using a combination of computational methods, including extensive molecular dynamics simulations and molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) energetic analysis. The simulation results strongly suggest that the observed tendency of wild-type CVN to form domain-swapped dimers is the result of a previously unidentified cis-peptide bond present in the monomeric state. The energetic analysis additionally indicates that the highest-affinity ligand for CVN characterized to date (α-Man-(1,2)-α-Man-(1,2)-α-Man) is recognized asymmetrically by the two binding sites. Finally, we are able to provide a detailed map of the role of all binding site functional groups (both backbone and side chain) to various aspects of molecular recognition: general affinity for cognate ligands, specificity for distinct oligosaccharide targets, and the asymmetric recognition of α-Man-(1,2)-α-Man-(1,2)-α-Man. Taken as a whole, these results complement past experimental characterization (both structural and thermodynamic) to provide the most complete understanding of carbohydrate recognition by CVN to date. The results also provide strong support for the application of similar approaches to the understanding of other protein-carbohydrate complexes.

Nucleotide recognition by the initiation factor aIF5B: free energy simulations of a neoclassical GTPase

The GTPase aIF5B is a universally conserved initiation factor that assists ribosome assembly. Crystal structures of its nucleotide complexes, X-ray(GTP) and X-ray(GDP), are similar in the nucleotide vicinity, but differ in the orientation of a distant domain IV. This has led to two, contradictory, mechanistic models. One postulates that X-ray(GTP) and X-ray(GDP) are, respectively, the active, "ON" and the inactive, "OFF" states; the other postulates that both structures are OFF, whereas the ON state is still uncharacterized. We study GTP/GDP binding using molecular dynamics and a continuum electrostatic free energy method. We predict that X-ray(GTP) has a ≈ 3 kcal/mol preference to bind GDP, apparently contradicting its assignment as ON. However, the preference arises mainly from a single, nearby residue from the switch 2 motif: Glu81, which becomes protonated upon GTP binding, with a free energy cost of about 4 kcal/mol. We then propose a different model, where Glu81 protonation/deprotonation defines the ON/OFF states. With this model, the X-ray(GTP):GTP complex, with its protonated Glu81, is ON, whereas X-ray(GTP):GDP is OFF. The model postulates that distant conformational changes such as domain IV rotation are "uncoupled" from GTP/GDP exchange and do not affect the relative GTP/GDP binding affinities. We analyze the model using a general thermodynamic framework for GTPases. It yields rather precise predictions for the nucleotide specificities of each state, and the state specificities of each nucleotide, which are roughly comparable to the homologues IF2 and aIF2, despite the lack of any conformational switching in the model.

Conformational selection through electrostatics: Free energy simulations of GTP and GDP binding to archaeal initiation factor 2

Archaeal Initiation Factor 2 is a GTPase involved in protein biosynthesis. In its GTP-bound, "ON" conformation, it binds an initiator tRNA and carries it to the ribosome. In its GDP-bound, "OFF" conformation, it dissociates from tRNA. To understand the specific binding of GTP and GDP and their dependence on the conformational state, molecular dynamics free energy simulations were performed. The ON state specificity was predicted to be weak, with a GTP/GDP binding free energy difference of -1 kcal/mol, favoring GTP. The OFF state specificity is larger, 4 kcal/mol, favoring GDP. The overall effects result from a competition among many interactions in several complexes. To interpret them, we use a simpler, dielectric continuum model. Several effects are robust with respect to the model details. Both nucleotides have a net negative charge, so that removing them from solvent into the binding pocket carries a desolvation penalty, which is large for the ON state, and strongly disfavors GTP binding compared to GDP. Short-range interactions between the additional GTP phosphate group and ionized sidechains in the binding pocket offset most, but not all of the desolvation penalty; more distant groups also contribute significantly, and the switch 1 loop only slightly. The desolvation penalty is lower for the more open, wetter OFF state, and the GTP/GDP difference much smaller. Short-range interactions in the binding pocket and with more distant groups again make a significant contribution. Overall, the simulations help explain how conformational selection is achieved with a single phosphate group.

Unique example of amyloid aggregates stabilized by main chain H-bond instead of the steric zipper: molecular dynamics study of the amyloidogenic segment of amylin wild-type and mutants

Most proteins do not aggregate while in their native functional states. However, they may be disturbed from their native conformation by certain change in the environment, and form unwanted oligomeric or polymeric aggregates. Recent experimental data demonstrate that soluble oligomers of amyloidogenic proteins are responsible for amyloidosis and its cytotoxicity. Human islet amyloid polypeptide (IAPP or amylin) is a 37-residue hormone found as fibrillar deposits in pancreatic extracts of nearly all type II diabetics. In this study we performed in silico mutation analysis to examine the stability of the double layer five strand aggregates formed by heptapeptide NNFGAIL segment from amyline peptide. This segment is one of the shortest fragments that can form amyloid fibrils similar to those formed by the full length peptide. The mutants obtained by single glycine replacement were also studied to investigate the specificity of the dry self-complementary interface between the neighboring β-sheet layers. The molecular dynamics simulations of the aggregates run for 20 ns at 330 K, the degree of the aggregate disassembly was investigated using several geometry analysis tools: the root mean square deviations of the C(α) atoms, root mean square fluctuations per residue, twist angles, interstrand distances, fraction of the secondary structure elements, and number of H-bonds. The analysis shows that most mutations make the aggregates unstable, and their stabilities were dependent to a large extent on the position of replaced residues. Our mutational simulations are in agreement with the pervious experimental observations. We also used free binding energy calculations to determine the role of different components: nonpolar effects, electrostatics and entropy in binding. Nonpolar effects remained consistently more favorable in wild type and mutants reinforcing the importance of hydrophobic effects in protein-protein binding. While entropy systematically opposed binding in all cases, there was no clear trend in the entropy difference between wildtype and glycine mutants. Free energy decomposition shows residues situated at the interface were found to make favorable contributions to the peptide-peptide association. The study of the wild type and mutants in an explicit solvent could provide valuable insight into the future computer guided design efforts for the amyloid aggregation inhibitor.

Computational protein design with a generalized Born solvent model: application to Asparaginyl-tRNA synthetase

Computational Protein Design (CPD) is a promising method for high throughput protein and ligand mutagenesis. Recently, we developed a CPD method that used a polar-hydrogen energy function for protein interactions and a Coulomb/Accessible Surface Area (CASA) model for solvent effects. We applied this method to engineer aspartyl-adenylate (AspAMP) specificity into Asparaginyl-tRNA synthetase (AsnRS), whose substrate is asparaginyl-adenylate (AsnAMP). Here, we implement a more accurate function, with an all-atom energy for protein interactions and a residue-pairwise generalized Born model for solvent effects. As a first test, we compute aminoacid affinities for several point mutants of Aspartyl-tRNA synthetase (AspRS) and Tyrosyl-tRNA synthetase and stability changes for three helical peptides and compare with experiment. As a second test, we readdress the problem of AsnRS aminoacid engineering. We compare three design criteria, which optimize the folding free-energy, the absolute AspAMP affinity, and the relative (AspAMP-AsnAMP) affinity. The sequences and conformations are improved with respect to our previous, polar-hydrogen/CASA study: For several designed complexes, the AspAMP carboxylate forms three interactions with a conserved arginine and a designed lysine, as in the active site of the AspRS:AspAMP complex. The conformations and interactions are well maintained in molecular dynamics simulations and the sequences have an inverted specificity, favoring AspAMP over AsnAMP. The method is not fully successful, since experimental measurements with the seven most promising sequences show that they do not catalyze at a detectable level the adenylation of Asp (or Asn) with ATP. This may be due to weak AspAMP binding and/or disruption of transition-state stabilization.

Kinetics, in silico docking, molecular dynamics, and MM-GBSA binding studies on prototype indirubins, KT5720, and staurosporine as phosphorylase kinase ATP-binding site inhibitors: the role of water molecules examined

With an aim toward glycogenolysis control in Type 2 diabetes, we have investigated via kinetic experiments and computation the potential of indirubin (IC₅₀ > 50 μM), indirubin-3'-oxime (IC₅₀ = 144 nM), KT5720 (K(i) = 18.4 nM) and staurosporine (K(i) = 0.37 nM) as phosphorylase kinase (PhKγtrnc) ATP-binding site inhibitors, with the latter two revealed as potent inhibitors in the low nM range. Because of lack of structural information, we have exploited information from homologous kinase complexes to direct in silico calculations (docking, molecular dynamics, and MMGBSA) to predict the binding characteristics of the four ligands. All inhibitors are predicted to bind in the same active site area as the ATP adenine ring, with binding dominated by hinge region hydrogen bonds to Asp104:O and Met106:O (all four ligands) and also Met106:NH (for the indirubins). The PhKγtrnc-staurosporine complex has the greatest number of receptor-ligand hydrogen bonds, while for the indirubin-3'-oxime and KT5720 complexes there is an important network of interchanging water molecules bridging inhibitor-enzyme contacts. The MM-GBSA results revealed the source of staurosporine's low nM potency to be favorable electrostatic interactions, while KT5720 has strong van der Waals contributions. KT5720 interacts with the greatest number of protein residues either by direct or 1-water bridged hydrogen bond interactions, and the potential for more selective PhK inhibition based on a KT5720 analogue has been established. Including receptor flexibility in Schrödinger induced-fit docking calculations in most cases correctly predicted the binding modes as compared with the molecular dynamics structures; the algorithm was less effective when there were key structural waters bridging receptor-ligand contacts.

Molecular Basis for Drug Resistance in HIV-1 Protease

HIV-1 protease is one of the major antiviral targets in the treatment of patients infected with HIV-1. The nine FDA approved HIV-1 protease inhibitors were developed with extensive use of structure-based drug design, thus the atomic details of how the inhibitors bind are well characterized. From this structural understanding the molecular basis for drug resistance in HIV-1 protease can be elucidated. Selected mutations in response to therapy and diversity between clades in HIV-1 protease have altered the shape of the active site, potentially altered the dynamics and even altered the sequence of the cleavage sites in the Gag polyprotein. All of these interdependent changes act in synergy to confer drug resistance while simultaneously maintaining the fitness of the virus. New strategies, such as incorporation of the substrate envelope constraint to design robust inhibitors that incorporate details of HIV-1 protease’s function and decrease the probability of drug resistance, are necessary to continue to effectively target this key protein in HIV-1 life cycle.

A Mixed QM/MM Scoring Function to Predict Protein-Ligand Binding Affinity

Computational methods for predicting protein-ligand binding free energy continue to be popular as a potential cost-cutting method in the drug discovery process. However, accurate predictions are often difficult to make as estimates must be made for certain electronic and entropic terms in conventional force field based scoring functions. Mixed quantum mechanics/molecular mechanics (QM/MM) methods allow electronic effects for a small region of the protein to be calculated, treating the remaining atoms as a fixed charge background for the active site. Such a semi-empirical QM/MM scoring function has been implemented in AMBER using DivCon and tested on a set of 23 metalloprotein-ligand complexes, where QM/MM methods provide a particular advantage in the modeling of the metal ion. The binding affinity of this set of proteins can be calculated with an R(2) of 0.64 and a standard deviation of 1.88 kcal/mol without fitting and 0.71 and a standard deviation of 1.69 kcal/mol with fitted weighting of the individual scoring terms. In this study we explore using various methods to calculate terms in the binding free energy equation, including entropy estimates and minimization standards. From these studies we found that using the rotational bond estimate to ligand entropy results in a reasonable R(2) of 0.63 without fitting. We also found that using the ESCF energy of the proteins without minimization resulted in an R(2) of 0.57, when using the rotatable bond entropy estimate.

Dimer-tetramer transition controls RUNX1/ETO leukemogenic activity

RUNX1/ETO, the fusion protein resulting from the chromosomal translocation t(8;21), is one of the most frequent translocation products in acute myeloid leukemia. Several in vitro and in vivo studies have shown that the homo-tetramerization domain of ETO, the nervy homology region 2 (NHR2), is essential for RUNX1/ETO oncogenic activity. We analyzed the energetic contribution of individual amino acids within the NHR2 to RUNX1/ETO dimer-tetramer transition and found a clustered area of 5 distinct amino acids with strong contribution to the stability of tetramers. Substitution of these amino acids abolishes tetramer formation without affecting dimer formation. Similar to RUNX1/ETO monomers, dimers failed to bind efficiently to DNA and to alter expression of RUNX1-dependent genes. RUNX1/ETO dimers do not block myeloid differentiation, are unable to enhance the self-renewal capacity of hematopoietic progenitors, and fail to induce leukemia in a murine transplantation model. Our data reveal the existence of an essential structural motif (hot spot) at the NHR2 dimer-tetramer interface, suitable for a molecular intervention in t(8;21) leukemias.

Energetic decomposition with the generalized-born and Poisson-Boltzmann solvent models: lessons from association of G-protein components

Continuum electrostatic models have been shown to be powerful tools in providing insight into the energetics of biomolecular processes. While the Poisson-Boltzmann (PB) equation provides a theoretically rigorous approach to computing electrostatic free energies of solution in such a model, computational cost makes its use for large ensembles of states impractical. The generalized-Born (GB) approximation provides a much faster alternative, although with a weaker theoretical framework. While much attention has been given to how GB recapitulates PB energetics for the overall stability of a biomolecule or the affinity of a complex, little attention has been given to how the contributions of individual functional groups are captured by the two methods. Accurately capturing these individual electrostatic components is essential both for the development of a mechanistic understanding of biomolecular processes and for the design of variant sequences and structures with desired properties. Here, we present a detailed comparison of the group-wise decomposition of both PB and GB electrostatic free energies of binding, using association of various components of the heterotrimeric-G-protein complex as a model. We find that, while net binding free energies are strongly correlated in the two models, the correlations of individual group contributions are highly variable; in some cases, strong correlation is seen, while in others, there is essentially none. Structurally, the GB model seems to capture the magnitude of direct, short-range electrostatic interactions quite well but performs more poorly with moderate-range "action-at-a-distance" interactions--GB has a tendency to overestimate solvent screening over moderate distances, and to underestimate the costs of desolvating charged groups somewhat removed from the binding interface. Despite this, however, GB does seem to be quite effective as a predictor of those groups that will be computed to be most significant in a PB-based model.

Computational design of protein-ligand binding: modifying the specificity of asparaginyl-tRNA synthetase

A method for computational design of protein-ligand interactions is implemented and tested on the asparaginyl- and aspartyl-tRNA synthetase enzymes (AsnRS, AspRS). The substrate specificity of these enzymes is crucial for the accurate translation of the genetic code. The method relies on a molecular mechanics energy function and a simple, continuum electrostatic, implicit solvent model. As test calculations, we first compute AspRS-substrate binding free energy changes due to nine point mutations, for which experimental data are available; we also perform large-scale redesign of the entire active site of each enzyme (40 amino acids) and compare to experimental sequences. We then apply the method to engineer an increased binding of aspartyl-adenylate (AspAMP) into AsnRS. Mutants are obtained using several directed evolution protocols, where four or five amino acid positions in the active site are randomized. Promising mutants are subjected to molecular dynamics simulations; Poisson-Boltzmann calculations provide an estimate of the corresponding, AspAMP, binding free energy changes, relative to the native AsnRS. Several of the mutants are predicted to have an inverted binding specificity, preferring to bind AspAMP rather than the natural substrate, AsnAMP. The computed binding affinities are significantly weaker than the native, AsnRS:AsnAMP affinity, and in most cases, the active site structure is significantly changed, compared to the native complex. This almost certainly precludes catalytic activity. One of the designed sequences has a higher affinity and more native-like structure and may represent a valid candidate for Asp activity.

Alchemical free energy simulations for biological complexes: powerful but temperamental

Free energy simulations compare multiple ligand:receptor complexes by "alchemically" transforming one into another, yielding binding free energy differences. Since their introduction in the 1980s, many technical and theoretical obstacles were surmounted, and the method ("MDFE," since molecular dynamics are often used) has matured into a powerful tool. We describe its current status, its effectiveness, and the challenges it faces. MDFE has provided chemical accuracy for many systems but remains expensive, with significant human overhead costs. The bottlenecks have shifted, partly due to increased computer power. To study diverse sets of ligands, force field availability and accuracy can be a major difficulty. Another difficulty is the frequent need to consider multiple states, related to sidechain protonation or buried waters, for example. Sophisticated, automated methods to sample these states are maturing, such as constant pH simulations. Meanwhile, combinations of MDFE and simpler approaches, like continuum dielectric models, can be very effective. As illustrations, we show how, with careful force field parameterization, MDFE accurately predicts binding specificities between complex tetracycline ligands and their targets. We describe substrate binding to the aspartyl-tRNA synthetase enzyme, where many distinct electrostatic states play a role, and a histidine and a Mg(2+) ion act as coupled switches that help enforce a strict preference for the aspartate substrate, relative to several analogs. Overall, MDFE has achieved a predictive status, where novel ligands can be studied and molecular recognition elucidated in depth. It should play an increasing role in the analysis of complex cellular processes and biomolecular engineering.

Exit strategies for charged tRNA from GluRS

For several class I aminoacyl-tRNA synthetases (aaRSs), the rate-determining step in aminoacylation is the dissociation of charged tRNA from the enzyme. In this study, the following factors affecting the release of the charged tRNA from aaRSs are computationally explored: the protonation states of amino acids and substrates present in the active site, and the presence and the absence of AMP and elongation factor Tu. Through molecular modeling, internal pK(a) calculations, and molecular dynamics simulations, distinct, mechanistically relevant post-transfer states with charged tRNA bound to glutamyl-tRNA synthetase from Thermus thermophilus (Glu-tRNA(Glu)) are considered. The behavior of these nonequilibrium states is characterized as a function of time using dynamical network analysis, local energetics, and changes in free energies to estimate transitions that occur during the release of the tRNA. The hundreds of nanoseconds of simulation time reveal system characteristics that are consistent with recent experimental studies. Energetic and network results support the previously proposed mechanism in which the transfer of amino acid to tRNA is accompanied by the protonation of AMP to H-AMP. Subsequent migration of proton to water reduces the stability of the complex and loosens the interface both in the presence and in the absence of AMP. The subsequent undocking of AMP or tRNA then proceeds along thermodynamically competitive pathways. Release of the tRNA acceptor stem is further accelerated by the deprotonation of the alpha-ammonium group on the charging amino acid. The proposed general base is Glu41, a residue binding the alpha-ammonium group that is conserved in both structure and sequence across nearly all class I aaRSs. This universal handle is predicted through pK(a) calculations to be part of a proton relay system for destabilizing the bound charging amino acid following aminoacylation. Addition of elongation factor Tu to the aaRS.tRNA complex stimulates the dissociation of the tRNA core and the tRNA acceptor stem.

Decomposing the energetic impact of drug resistant mutations in HIV-1 protease on binding DRV

Darunavir (DRV) is a high affinity (4.5×10(-12) M, ΔG = -15.2 kcal/mol) HIV-1 protease inhibitor. Two drug-resistant protease variants FLAP+ (L10I, G48V, I54V, V82A) and ACT (V82T, I84V) decrease the binding affinity with DRV by 1.0 kcal/mol and 1.6 kcal/mol respectively. In this study the absolute and relative binding free energies of DRV with wild-type protease, FLAP+ and ACT were calculated with MM-PB/GBSA and thermodynamic integration methods, respectively. Free energy decomposition elucidated that the mutations conferred resistance by distorting the active site of HIV-1 protease so that the residues that lost binding free energy were not limited to the sites of mutation. Specifically the bis-tetrahydrofuranylurethane moiety of DRV maintained interactions with the FLAP+ and ACT variants, whereas the 4 - amino phenyl group lost more binding free energy with the protease in the FLAP+ and ACT complexes than in the wild-type protease which could account for the majority of the loss in binding free energy. This suggested that replacement of the 4 - amino phenyl group might generate new inhibitors less susceptible to the drug resistant mutations.

Molecular dynamics simulations show that conformational selection governs the binding preferences of imatinib for several tyrosine kinases

Tyrosine kinases transmit cellular signals through a complex mechanism, involving their phosphorylation and switching between inactive and active conformations. The cancer drug imatinib binds tightly to several homologous kinases, including Abl, but weakly to others, including Src. Imatinib specifically targets the inactive, so-called "DFG-out" conformation of Abl, which differs from the preferred, "DFG-in" conformation of Src in the orientation of a conserved Asp-Phe-Gly (DFG) activation loop. However, recent x-ray structures showed that Src can also adopt the DFG-out conformation and uses it to bind imatinib. The Src/Abl-binding free energy difference can thus be decomposed into two contributions. Contribution i measures the different protein-imatinib interactions when either kinase is in its DFG-out conformation. Contribution ii depends on the ability of imatinib to select or induce this conformation, i.e. on the relative stabilities of the DFG-out and DFG-in conformations of each kinase. Neither contribution has been measured experimentally. We use molecular dynamics simulations to show that contribution i is very small, 0.2 +/- 0.6 kcal/mol; imatinib interactions are very similar in the two kinases, including long range electrostatic interactions with the imatinib positive charge. Contribution ii, deduced using the experimental binding free energy difference, is much larger, 4.4 +/- 0.9 kcal/mol. Thus, conformational selection, easy in Abl, difficult in Src, underpins imatinib specificity. Contribution ii has a simple interpretation; it closely approximates the stability difference between the DFG-out and DFG-in conformations of apo-Src. Additional calculations show that conformational selection also governs the relative binding of imatinib to the kinases c-Kit and Lck. These results should help clarify the current framework for engineering kinase signaling.

Tetracycline-tet repressor binding specificity: insights from experiments and simulations

Tetracycline (Tc) antibiotics have been put to new uses in the construction of artificial gene regulation systems, where they bind to the Tet repressor protein (TetR) and modulate its affinity for DNA. Many Tc variants have been produced, both to overcome bacterial resistance and to achieve a broad range of binding strengths. To better understand TetR-Tc binding, we investigate a library of 16 tetracyclines, using fluorescence experiments and molecular dynamics free energy simulations (MDFE). The relative TetR binding free energies are computed by reversibly transforming one Tc variant into another during the simulation, with no adjustable parameters. The chemical variations involve polar and nonpolar substitutions along one entire edge of the elongated Tc structure, which provides many of the protein-ligand contacts. The binding constants span five orders of magnitude. The simulations reproduce the experimental binding free energies, when available, within the uncertainty of either method (+/-0.5 kcal/mol), and reveal many additional details. Contributions of individual Tc substituents are evaluated, along with their additivity and transferability among different positions on the Tc scaffold; differences between D- and B-class repressors are quantified. With increasing computer power, the MDFE approach provides an attractive complement to experiment and should play an increasing role in the understanding and engineering of protein-ligand recognition.

Relative free energy of binding between antimicrobial peptides and SDS or DPC micelles

We present relative binding free energy calculations for six antimicrobial peptide-micelle systems, three peptides interacting with two types of micelles. The peptides are the scorpion derived antimicrobial peptide (AMP), IsCT and two of its analogues. The micelles are dodecylphosphatidylcholine (DPC) and sodium dodecylsulphate (SDS) micelles. The interfacial electrostatic properties of DPC and SDS micelles are assumed to be similar to those of zwitterionic mammalian and anionic bacterial membrane interfaces, respectively. We test the hypothesis that the binding strength between peptides and the anionic micelle SDS can provide information on peptide antimicrobial activity, since it is widely accepted that AMPs function by binding to and disrupting the predominantly anionic lipid bilayer of the bacterial cytoplasmic membrane. We also test the hypothesis that the binding strength between peptides and the zwitterionic micelle DPC can provide information on peptide haemolytic activities, since it is accepted that they also bind to and disrupt the zwitterionic membrane of mammalian cells. Equilibrium structures of the peptides, micelles and peptide-micelle complexes are obtained from more than 300 ns of molecular dynamics simulations. A thermodynamic cycle is introduced to compute the binding free energy from electrostatic, non-electrostatic and entropic contributions. We find relative binding free energy strengths between peptides and SDS to correlate with the experimentally measured rankings for peptide antimicrobial activities, and relative free energy binding strengths between peptides and DPC to correlate with the observed rankings for peptide haemolytic toxicities. These findings point to the importance of peptide-membrane binding strength for antimicrobial activity and haemolytic activity.

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.

Calculation of free-energy differences by confinement simulations. Application to peptide conformers

Conformational free-energy differences are key quantities for understanding important phenomena in molecular biology that involve large structural changes of macromolecules. In this paper, an improved version of the confinement approach, which is based on earlier developments, is used to determine the free energy of individual molecular states by progressively restraining the corresponding molecular structures to pure harmonic basins, whose absolute free energy can be computed by normal-mode analysis. The method is used to calculate the free-energy difference between two conformational states of the alanine dipeptide in vacuo, and of the beta-hairpin from protein G with an implicit solvation model. In all cases, the confinement results are in excellent agreement with the ones obtained from converged equilibrium molecular dynamics simulations, which have a much larger computational cost. The systematic and statistical errors of the results are evaluated and the origin of the errors is identified. The sensitivity of the calculated free-energy differences to structure-based definitions of the molecular states is discussed. A variant of the method, which closes the thermodynamic cycle by a quasi-harmonic rather than harmonic analysis, is introduced. The latter is proposed for possible use with explicit solvent simulations.

Insights into affinity and specificity in the complexes of alpha-lytic protease and its inhibitor proteins: binding free energy from molecular dynamics simulation

We report the binding free energy calculation and its decomposition for the complexes of alpha-lytic protease and its protein inhibitors using molecular dynamics simulation. Standard mechanism serine protease inhibitors eglin C and OMTKY3 are known to have strong binding affinity for many serine proteases. Their binding loops have significant similarities, including a common P1 Leu as the main anchor in the binding interface. However, recent experiments demonstrate that the two inhibitors have vastly different affinity towards alpha-lytic protease (ALP), a bacterial serine protease. OMTKY3 inhibits the enzyme much more weakly (by approximately 10(6) times) than eglin C. Moreover, a variant of OMTKY3 with five mutations, OMTKY3M, has been shown to inhibit 10(4) times more strongly than the wild-type inhibitor. The underlying mechanisms for the unusually large difference in binding affinities and the effect of mutation are not well understood. Here we use molecular dynamics simulation with molecular mechanics-Poisson Boltzmann/surface area method (MM-PB/SA) to investigate quantitatively the binding specificity. The calculated absolute binding free energies correctly differentiate the thermodynamic stabilities of these protein complexes, but the magnitudes of the binding affinities are systematically overestimated. Analysis of the binding free energy components provides insights into the molecular mechanism of binding specificity. The large DeltaDeltaG(bind) between eglin C and wild type OMTKY3 towards ALP is mainly attributable to the stronger nonpolar interactions in the ALP-eglin C complex, arising from a higher degree of structural complementarity. Here the electrostatic interaction contributes to a lesser extent. The enhanced inhibition in the penta-mutant OMTKY3M over its wild type is entirely due to an overall improvement in the solvent-mediated electrostatic interactions in the ALP-OMTKY3M complex. The results suggest that for these protein-complexes and similar enzyme-inhibitor systems (1) the binding is driven by nonpolar interactions, opposed by overall electrostatic and solute entropy contributions; (2) binding specificity can be tuned by improving the complementarity in electrostatics between two associating proteins. Binding free energy decomposition into contributions from individual protein residues provides additional detailed information on the structural determinants and subtle conformational changes responsible for the binding specificity.

Binding of tetracyclines to elongation factor Tu, the Tet repressor, and the ribosome: a molecular dynamics simulation study

Tetracycline (Tc) is a broad-spectrum antibiotic that kills bacteria by interrupting protein biosynthesis. It is thought that the bacteriostatic action of Tc is associated with its binding to the acceptor site (or A site) in the bacterial ribosome, interfering with the attachment of aminoacyl-tRNA. Recently, however, the crystal structure of a complex between Tc and trypsin-modified elongation factor Tu (tm-EF-Tu) was determined, raising the question of whether Tc binding to EF-Tu has a role in its inhibition of protein synthesis. We address this question using computer simulations. As controls, we first compute relative ribosome binding free energies for seven Tc variants for which experimental data are available, obtaining good agreement. We then consider the binding of Tc to both the trypsin-modified and unmodified EF-Tu-GDP complexes. We show that the direct contribution of EF-Tu to the binding free energy is negligible; rather, the binding can be solely attributed to interactions of Tc with a bridging Mg(2+) ion and the GDP phosphate groups. The effects of trypsin modification are modest. Further, our calculations show that EF-Tu does not exhibit any binding preference for Tc over the nonantibiotic, 4-dedimethyl-Tc, and EF-Tu does not bind the Tc analogue tigecycline, which is a potent antibiotic. In contrast, both the ribosome and the Tet Repressor protein (involved in Tc resistance) do show a binding preference for Tc over 4-dedimethyl-Tc, and the ribosome prefers to bind tigecycline over Tc. Overall, our results provide insights into the binding properties of tetracyclines and support the idea that EF-Tu is not their primary target.

Computer simulations reveal a novel nucleotide-type binding orientation for ellipticine-based anticancer c-kit kinase inhibitors

Receptor tyrosine kinase (RTK) enzymes regulate cell signaling pathways and so are an important target for cancer chemotherapy. Current inhibitors of c-kit, a key RTK stem cell factor receptor, are inactive against the most common mutated variant Asp816Val, associated with highly malignant cancers. Recent combined experimental/simulation work has highlighted the utility of the ellipticine pharmacore in inhibiting mutant c-kit, and the present simulation study applies a combination of high-level simulation tools to probe further the binding of ellipticine-based derivatives to c-kit. We find a large preference for protonation of bound ellipticine, which stabilizes the negative protein residues that coordinated ADP.Mg (2+) in the native complex. The resulting ellipticine inhibitor binding mode resembles the native nucleotide complex and serves to explain some existing experimental data on binding specificities, indicating that functionalization at the C4/C5 sites of ellipticine derivatives may be important for the design of novel nucleotide analogues that inhibit mutant c-kit.

Probing electrostatic interactions and ligand binding in aspartyl-tRNA synthetase through site-directed mutagenesis and computer simulations

Faithful genetic code translation requires that each aminoacyl-tRNA synthetase recognise its cognate amino acid ligand specifically. Aspartyl-tRNA synthetase (AspRS) distinguishes between its negatively-charged Asp substrate and two competitors, neutral Asn and di-negative succinate, using a complex network of electrostatic interactions. Here, we used molecular dynamics simulations and site-directed mutagenesis experiments to probe these interactions further. We attempt to decrease the Asp/Asn binding free energy difference via single, double and triple mutations that reduce the net positive charge in the active site of Escherichia coli AspRS. Earlier, Glutamine 199 was changed to a negatively-charged glutamate, giving a computed reduction in Asp affinity in good agreement with experiment. Here, Lysine 198 was changed to a neutral leucine; then, Lys198 and Gln199 were mutated simultaneously. Both mutants are predicted to have reduced Asp binding and improved Asn binding, but the changes are insufficient to overcome the initial, high specificity of the native enzyme, which retains a preference for Asp. Probing the aminoacyl-adenylation reaction through pyrophosphate exchange experiments, we found no detectable activity for the mutant enzymes, indicating weaker Asp binding and/or poorer transition state stabilization. The simulations show that the mutations' effect is partly offset by proton uptake by a nearby histidine. Therefore, we performed additional simulations where the nearby Histidines 448 and 449 were mutated to neutral or negative residues: (Lys198Leu, His448Gln, His449Gln), and (Lys198Leu, His448Glu, His449Gln). This led to unexpected conformational changes and loss of active site preorganization, suggesting that the AspRS active site has a limited structural tolerance for electrostatic modifications. The data give insights into the complex electrostatic network in the AspRS active site and illustrate the difficulty in engineering charged-to-neutral changes of the preferred ligand.