Investigating the binding specificity of U1A-RNA by computational mutagenesis

A study on allosteric communication in U1A-snRNA binding interactions: network analysis combined with molecular dynamics data

The allosteric regulation during the binding interactions between small nuclear RNAs (snRNAs) and the associated protein factors is critical to the function of spliceosomes in alternative RNA splicing. Although network models combined with molecular dynamics simulations have shown to be powerful tools for the analysis of protein allostery, the atomic-level simulations are, however, too expensive and with limited accuracy for the large-size systems. In this work, we use a residual network model combined with a coarse-grained Gaussian network model (GNM) to investigate the binding interactions between the snRNA and the human U1A protein which is a major component of the spliceosomal U1 small nuclear ribonucleoprotein particle, and to identify the residues that play an important role in the allosteric communication in U1A during this process. We also utilize the Girvan-Newman method to detect the structural organization in U1A-snRNA recognition and interactions. Our results reveal that: (Ι) not only the residues at the binding sites that are traditionally considered to play a major role in U1A-snRNA association, but those residues that are far away from the RNA binding interface participate in the U1A's allosteric signal transmission induced by the RNA binding; (Π) the structure of U1A protein is well organized with different communities acting different roles for its RNA binding and allosteric regulation. The study demonstrates that the combination of the residual network and elastic network models is an effective and efficient method which can be readily extended to the investigation of the allosteric communication for other macromolecular interaction systems.

Interpreting the Dynamics of Binding Interactions of snRNA and U1A Using a Coarse-Grained Model

The binding interactions of small nuclear RNAs (snRNA) and the associated protein factors are critical to the function of spliceosomes in alternatively splicing primary RNA transcripts. Although molecular dynamics simulations are a powerful tool to interpret the mechanism of biological processes, the atomic-level simulations are, however, too expensive and with limited accuracy for the large-size systems, such as snRNA-protein complexes. We extend the coarse-grained Gaussian network model, which models the RNA-protein complexes as a harmonic chain of Cα, P, and O4' atoms, to investigating the impact of the snRNA-binding interaction on the dynamic stability of the human U1A protein, which is a major component of the spliceosomal U1 small nuclear ribonucleoprotein particle. The results reveal that the first and third loops and the C-terminal helix regions of the U1A domain undergo a significant loss of flexibility upon the RNA binding due to the forming of mostly electrostatic and hydrogen bond interactions with RNA 5' stem and loop. By examining the residues whose mutations significantly change the binding free energy between U1A and snRNA, the Gaussian network model-based calculations show that not only the residues at the binding sites that are traditionally considered to play a major role in U1A-RNA association but also those residues that are far away from the RNA-binding interface can participate in the long-range allosteric signal transmission; these calculations are quantitatively consistent with the data observed in the recent snRNA binding experiments. The study demonstrates a useful avenue to utilize the simplified elastic network model to investigate the dynamics characteristics of the biologically important macromolecular interactions.

ParFit: A Python-Based Object-Oriented Program for Fitting Molecular Mechanics Parameters to ab Initio Data

A newly created object-oriented program for automating the process of fitting molecular-mechanics parameters to ab initio data, termed ParFit, is presented. ParFit uses a hybrid of deterministic and stochastic genetic algorithms. ParFit can simultaneously handle several molecular-mechanics parameters in multiple molecules and can also apply symmetric and antisymmetric constraints on the optimized parameters. The simultaneous handling of several molecules enhances the transferability of the fitted parameters. ParFit is written in Python, uses a rich set of standard and nonstandard Python libraries, and can be run in parallel on multicore computer systems. As an example, a series of phosphine oxides, important for metal extraction chemistry, are parametrized using ParFit. ParFit is in an open source program available for free on GitHub ( https://github.com/fzahari/ParFit ).

Computational alanine scanning mutagenesis--an improved methodological approach

Alanine scanning mutagenesis of protein-protein interfacial residues can be applied to a wide variety of protein complexes to understand the structural and energetic characteristics of the hot-spots. Binding free energies have been estimated with reasonable accuracy with empirical methods, such as Molecular Mechanics/Poisson-Boltzmann surface area (MM-PBSA), and with more rigorous computational approaches like Free Energy Perturbation (FEP) and Thermodynamic Integration (TI). The main objective of this work is the development of an improved methodological approach, with less computational cost, that predicts accurately differences in binding free energies between the wild-type and alanine mutated complexes (DeltaDeltaG(binding)). The method was applied to three complexes, and a mean unsigned error of 0.80 kcal/mol was obtained in a set of 46 mutations. The computational method presented here achieved an overall success rate of 80% and an 82% success rate in residues for which alanine mutation causes an increase in the binding free energy > 2.0 kcal/mol (warm- and hot-spots). This fully atomistic computational methodological approach consists in a computational Molecular Dynamics simulation protocol performed in a continuum medium using the Generalized Born model. A set of three different internal dielectric constants, to mimic the different degree of relaxation of the interface when different types of amino acids are mutated for alanine, have to be used for the proteins, depending on the type of amino acid that is mutated. This method permits a systematic scanning mutagenesis of protein-protein interfaces and it is capable of anticipating the experimental results of mutagenesis, thus guiding new experimental investigations.

Phenytoin is an estrogen receptor α-selective modulator that interacts with helix 12

RATIONALE

Phenytoin (Dilantin(®); DPH) is used to treat epilepsy but causes estrogen agonist-antagonist-like side effects. We investigated the interaction of phenytoin with estrogen receptors (ERs) α and β by computational molecular docking, ER competition binding, transcriptional assays, and biological actions, comparing outcomes with estradiol (E2), estrone (E1), and tamoxifen (TMX).

EXPERIMENTAL

(1) The DPH docking to 3-dimensional crystal structures of the ERα ligand-binding domain (LBD) showed a high degree of structural complementarity (-57.15 calculated energy units, approximating kcal/mol) with the ligand-binding pocket, including a contact at leucine (L540) in helix 12. Estrogen receptor β showed slightly less favorable interactions (-54.27 kcal/mol), without contacting L450. Estradiol, E1, and TMX contact points with ERα and ERβ do not include L450. (2) Cellular actions: Incubation of cells transfected with ERα or ERβ and a luciferase promoter phenytoin was several orders weaker than E2 as an agonist through ERα and had no effect through ERβ. However, phenytoin at clinical concentrations (10(-11) to 10(-6) mol/L) powerfully antagonized action of E2 on ERα-expressing cells. Similarly, phenytoin at clinically effective concentrations marginally induced alkaline phosphatase by ERα- and ERβ-expressing endometrial cancer cells but at doses well below clinical effectiveness blocked E2-induced alkaline phosphatase. (3) ER competition: In Scatchard plots comparing phenytoin with 17β-estradiol against endometrial cancer cell cytosol E2-alone more effectively displaced labeled E2 than phenytoin, but phenytoin was approximately equimolar effective to E2 in inhibiting E2's displacement of the radiolabel, further confirming that phenytoin is a strong E2 antagonist.

CONCLUSIONS

At clinically effective concentrations, phenytoin is a strong ERα cell antagonist but a many-fold weaker agonist. Although it interacts with ERβ LBD residues, phenytoin has no effects on ERβ-only expressing cells. Docking studies indicate phenytoin interacts with the ERα LBD at the hinge of helix 12 and could thereby interfere with the entry of other ER ligands or with the mobility of helix 12, either of which actions could explain phenytoin's antagonism of ER-mediated E2 actions. Our results suggest an explanation for the broad profile of phenytoin's actions and raise possibilities for the use of phenytoin or congeners in the clinical management of ERα-dependent conditions.

Importance of polar solvation and configurational entropy for design of antiretroviral drugs targeting HIV-1 protease

Both KNI-10033 and KNI-10075 are high affinity preclinical HIV-1 protease (PR) inhibitors with affinities in the picomolar range. In this work, the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method has been used to investigate the potency of these two HIV-1 PR inhibitors against the wild-type and mutated proteases assuming that potency correlates with the affinity of the drugs for the target protein. The decomposition of the binding free energy reveals the origin of binding affinities or mutation-induced affinity changes. Our calculations indicate that the mutation I50V causes drug resistance against both inhibitors. On the other hand, we predict that the mutant I84V causes drug resistance against KNI-10075 while KNI-10033 is more potent against the I84V mutant compared to wild-type protease. Drug resistance arises mainly from unfavorable shifts in van der Waals interactions and configurational entropy. The latter indicates that neglecting changes in configurational entropy in the computation of relative binding affinities as often done is not appropriate in general. For the bound complex PR(I50V)-KNI-10075, an increased polar solvation free energy also contributes to the drug resistance. The importance of polar solvation free energies is revealed when interactions governing the binding of KNI-10033 or KNI-10075 to the wild-type protease are compared to the inhibitors darunavir or GRL-06579A. Although the contributions from intermolecular electrostatic and van der Waals interactions as well as the nonpolar component of the solvation free energy are more favorable for PR-KNI-10033 or PR-KNI-10075 compared to PR-DRV or PR-GRL-06579A, both KNI-10033 and KNI-10075 show a similar affinity as darunavir and a lower binding affinity relative to GRL-06579A. This is because of the polar solvation free energy which is less unfavorable for darunavir or GRL-06579A relative to KNI-10033 or KNI-10075. The importance of the polar solvation as revealed here highlights that structural inspection alone is not sufficient for identifying the key contributions to binding affinities and affinity changes for the design of drugs but that solvation effects must be taken into account. A detailed understanding of the molecular forces governing binding and drug resistance might assist in the design of new inhibitors against HIV-1 PR variants that are resistant against current drugs.

U1A protein-stem loop 2 RNA recognition: prediction of structural differences from protein mutations

Molecular dynamics (MD) simulations were carried out to compare the free and bound structures of wild type U1A protein with several Phe56 mutant U1A proteins that bind the target stem loop 2 (SL2) RNA with a range of affinities. The simulations indicate the free U1A protein is more flexible than the U1A-RNA complex for both wild type and Phe56 mutant systems. A complete analysis of the hydrogen-bonding (HB) and non-bonded (VDW) interactions over the course of the MD simulations suggested that changes in the interactions in the free U1A protein caused by the Phe56Ala and Phe56Leu mutations may stabilize the helical character in loop 3, and contribute to the weak binding of these proteins to SL2 RNA. Compared with wild type, changes in HB and VDW interactions in Phe56 mutants of the free U1A protein are global, and include differences in β-sheet, loop 1 and loop 3 interactions. In the U1A-RNA complex, the Phe56Ala mutation leads to a series of differences in interactions that resonate through the complex, while the Phe56Leu and Phe56Trp mutations cause local differences around the site of mutation. The long-range networks of interactions identified in the simulations suggest that direct interactions and dynamic processes in both the free and bound forms contribute to complex stability.

Molecular recognition of RNA: challenges for modelling interactions and plasticity

There is growing interest in molecular recognition processes of RNA because of RNA's widespread involvement in biological processes. Computational approaches are increasingly used for analysing and predicting binding to RNA, fuelled by encouraging progress in developing simulation, free energy and docking methods for nucleic acids. These developments take into account challenges regarding the energetics of RNA-ligand binding, RNA plasticity, and the presence of water molecules and ions in the binding interface. Accordingly, we will detail advances in force field and scoring function development for molecular dynamics (MD) simulations, free energy computations and docking calculations of nucleic acid complexes. Furthermore, we present methods that can detect moving parts within RNA structures based on graph-theoretical approaches or normal mode analysis (NMA). As an example of the successful use of these developments, we will discuss recent structure-based drug design approaches that focus on the bacterial ribosomal A-site RNA as a drug target.

Prediction of salt and mutational effects on the association rate of U1A protein and U1 small nuclear RNA stem/loop II

We have developed a computational approach for predicting protein-protein association rates (Alsallaq and Zhou, Structure 2007, 15, 215). Here we expand the range of applicability of this approach to protein-RNA binding and report the first results for protein-RNA binding rates predicted from atomistic modeling. The system studied is the U1A protein and stem/loop II of the U1 small nuclear RNA. Experimentally it was observed that the binding rate is significantly reduced by increasing salt concentration while the dissociation changes little with salt concentration, and charges distant from the binding site make marginal contribution to the binding rate. These observations are rationalized. Moreover, predicted effects of salt and charge mutations are found to be in quantitative agreement with experimental results.

Affinity and specificity of protein U1A-RNA complex formation based on an additive component free energy model

An MM-GBSA computational protocol was used to investigate wild-type U1A-RNA and F56 U1A mutant experimental binding free energies. The trend in mutant binding free energies compared to wild-type is well-reproduced. Following application of a linear-response-like equation to scale the various energy components, the binding free energies agree quantitatively with observed experimental values. Conformational adaptation contributes to the binding free energy for both the protein and the RNA in these systems. Small differences in DeltaGs are the result of different and sometimes quite large relative contributions from various energetic components. Residual free energy decomposition indicates differences not only at the site of mutation, but throughout the entire protein. MM-GBSA and ab initio calculations performed on model systems suggest that stacking interactions may nearly, but not completely, account for observed differences in mutant binding affinities. This study indicates that there may be different underlying causes of ostensibly similar experimentally observed binding affinities of different mutants, and thus recommends caution in the interpretation of binding affinities and specificities purely by inspection.

Protein-Ligand Binding Affinity Predictions by Implicit Solvent Simulations: A Tool for Lead Optimization?

Continuum electrostatics is combined with rigorous free-energy calculations in an effort to deliver a reliable and efficient method for in silico lead optimization. The methodology is tested by calculation of the relative binding free energies of a set of inhibitors of neuraminidase, cyclooxygenase2, and cyclin-dependent kinase 2. The calculated free energies are compared to the results obtained with explicit solvent simulations and empirical scoring functions. For cyclooxygenase2, deficiencies in the continuum electrostatics theory are identified and corrected with a modified simulation protocol. For neuraminidase, it is shown that a continuum representation of the solvent leads to markedly different protein-ligand interactions compared to the explicit solvent simulations, and a reconciliation of the two protocols is problematic. Cyclin-dependent kinase 2 proves more challenging, and none of the methods employed in this study yield high quality predictions. Despite the differences observed, for these systems, the use of an implicit solvent framework to predict the ranking of congeneric inhibitors to a protein is shown to be faster, as accurate or more accurate than the explicit solvent protocol, and superior to empirical scoring schemes.

Flexible 5-guanidino-4-nitroimidazole DNA lesions: structures and thermodynamics

5-Guanidino-4-nitroimidazole (NI), derived from guanine oxidation by reactive oxygen and nitrogen species, contains an unusual flexible ring-opened structure, with nitro and guanidino groups which possess multiple hydrogen bonding capabilities. In vitro primer extension experiments with bacterial and mammalian polymerases show that NI incorporates C as well as A and G opposite the lesion, depending on the polymerase. To elucidate structural and thermodynamic properties of the mutagenic NI lesion, we have investigated the structure of the modified base itself and the NI-containing nucleoside with high-level quantum mechanical calculations and have employed molecular modeling and molecular dynamics simulations in solution for the lesion in B-DNA duplexes, with four partner bases opposite the NI. Our results show that NI adopts a planar structure at the damaged base level. However, in the nucleoside and in DNA duplexes, steric hindrance between the guanidino group and its linked sugar causes NI to be nonplanar. The NI lesion can adopt both syn and anti conformations on the DNA duplex level, with the guanidino group positioned in the DNA major and minor grooves, respectively; the specific preference depends on the partner base. On the basis of hydrogen bonding and stacking interactions, groove dimensions, and bending, we find that the least distorted NI-modified duplex contains partner C, consistent with observed incorporation of C opposite NI. However, hydrogen bonding interactions between NI and partner G or A are also found, which would be compatible with the observed mismatches.

Detailed microscopic study of the full zipA:FtsZ interface

Protein-protein interaction networks are very important for a wide range of biological processes. Crystallographic structures and mutational studies have generated a large number of information that allowed the discovery of energetically important determinants of specificity at intermolecular protein interfaces and the understanding of the structural and energetic characteristics of the binding hot spots. In this study we have used the improved MMPB/SA (molecular mechanics/Poisson-Boltzmann surface area) approach that combining molecular mechanics and continuum solvent permits to calculate the free energy differences upon alanine mutation. For a better understanding of the binding determinants of the complex formed between the FtsZ fragment and ZipA we extended the alanine scanning mutagenesis study to all interfacial residues of this complex. As a result, we present new mutations that allowed the discovery of residues for which the binding free energy differences upon alanine mutation are higher than 2.0 kcal/mol. We also observed the formation of a hydrophobic pocket with a high warm spot spatial complementarity between FtsZ and ZipA. Small molecules could be designed to bind to these amino acid residues hindering the binding of FtsZ to ZipA. Hence, these mutational data can be used to design new drugs to control more efficiently bacterial infections.

Functional analysis and molecular modeling show a preserved wild-type activity of p53C238Y

In human tumors, p53 is often disabled by mutations in its DNA-binding domain and is thus inactive as a transcription factor. Alternatively, MDM2 gene amplification or up-regulation represents a mechanism of p53 wild-type inactivation, mainly reported in soft tissue sarcomas. In a previous TP53 analysis carried out on sporadic and NF1-related malignant peripheral nerve sheath tumors, in two cases, we observed the occurrence of C238Y missense mutation, leading to p53 stabilization unexpectedly coupled with immunophenotypic MDM2 overexpression. To investigate this TP53 missense mutation not yet functionally characterized in mammalian cell, we did MDM2 Southern blot and p53(C238Y)/MDM2 biochemical and functional analyses followed by molecular modeling. The results showed a lack of MDM2 gene amplification, evidence of p53-MDM2 protein complexes, and presence of a p53 that retains the ability to become phosphorylated on Ser15 and to induce the transcription of p21(waf1). Additional molecular modeling data highlighted the structural similarities between p53(C238Y) and wild-type p53, further supporting that the p53(C238Y) mutant still retains functional wild-type p53 properties.

The role of RNA structure in the interaction of U1A protein with U1 hairpin II RNA

The N-terminal RNA Recognition Motif (RRM1) of the spliceosomal protein U1A interacting with its target U1 hairpin II (U1hpII) has been used as a paradigm for RRM-containing proteins interacting with their RNA targets. U1A binds to U1hpII via direct interactions with a 7-nucleotide (nt) consensus binding sequence at the 5' end of a 10-nt loop, and via hydrogen bonds with the closing C-G base pair at the top of the RNA stem. Using surface plasmon resonance (Biacore), we have examined the role of structural features of U1hpII in binding to U1A RRM1. Mutational analysis of the closing base pair suggests it plays a minor role in binding and mainly prevents "breathing" of the loop. Lengthening the stem and nontarget part of the loop suggests that the increased negative charge of the RNA might slightly aid association. However, this is offset by an increase in dissociation, which may be caused by attraction of the RRM to nontarget parts of the RNA. Studies of a single stranded target and RNAs with untethered loops indicate that structure is not very relevant for association but is important for complex stability. In particular, breaking the link between the stem and the 5' side of the loop greatly increases complex dissociation, presumably by hindering simultaneous contacts between the RRM and stem and loop nucleotides. While binding of U1A to a single stranded target is much weaker than to U1hpII, it occurs with nanomolar affinity, supporting recent evidence that binding of unstructured RNA by U1A has physiological significance.

Functional analyses and molecular modeling of two c-Kit mutations responsible for imatinib secondary resistance in GIST patients

Imatinib-acquired resistance related to the presence of secondary point mutations has become a frequent event in gastrointestinal stromal tumors. Here, transient transfection experiments with plasmids carrying two different KIT-acquired point mutations were performed along with immunoprecipitation of total protein extracts, derived from imatinib-treated and untreated cells. The molecular mechanics/Poisson Boltzmann surface area computational techniques were applied to study the interactions of the wild-type and mutated receptors with imatinib at the molecular level. Biochemical analyses showed KIT phosphorylation in cells transfected with vectors carrying the specific mutant genes. Imatinib treatment demonstrated that T670I was insensitive to the drug at all the applied concentrations, whereas V654A was inhibited by 6 microM of imatinib. The modeling of the mutated receptors revealed that both substitutions affect imatinib-binding site, but to a different extent: T670I substantially modifies the binding pocket, whereas V654A induces only relatively confined structural changes. We demonstrated that T670I and V654A cause indeed imatinib-acquired resistance and that the former is more resistant to imatinib than the latter. The application of molecular simulations allowed us to quantify the interactions between the mutated receptors and imatinib, and to propose a molecular rationale for this type of drug resistance.

Unraveling the Importance of Protein−Protein Interaction: Application of a Computational Alanine-Scanning Mutagenesis to the Study of the IgG1 Streptococcal Protein G (C2 Fragment) Complex

Alanine-scanning mutagenesis of protein-protein interfacial residues is a very important process for rational drug design. In this study, we have used the improved MM-PBSA approach that combining molecular mechanics and continuum solvent permits one to calculate the free energy differences through alanine mutation. To identify the binding determinants of the complex formed between the IgG1 (immunoglobulin-binding protein G) and protein G, we have extended the experimental alanine scanning mutagenesis study to both proteins of this complex and, therefore, to all interfacial residues of this binding complex. As a result, we present new residues that can be characterized as warm spots and, therefore, are important for complex formation. We have further increased the understanding of the functionality of this improved computational alanine-scanning mutagenesis approach testing its sensitivity to a protein-protein complex with an interface made up of residues mainly polar. In this study, we also have improved the method for the detection of an important amino acid residue that frequently constitutes a hot spot--tryptophan.

Role of aspartate 143 in Escherichia coli tRNA-guanine transglycosylase: alteration of heterocyclic substrate specificity

tRNA-guanine transglycosylase (TGT) is a key enzyme involved in the post-transcriptional modification of certain tRNAs in their anticodon wobble positions with queuine. To maintain the correct Watson-Crick base pairing properties of the wobble base (and hence proper translation of the genetic code), TGT must recognize its heterocyclic substrate with high specificity. The X-ray crystal structure of a eubacterial TGT bound to preQ1 [Romier, C., et al. (1996) EMBO J. 15, 2850-2857] suggested that aspartate 143 (Escherichia coli TGT numbering) was involved in heterocyclic substrate recognition. Subsequent mutagenic and computational modeling studies from our lab [Todorov, K. A., et al. (2005) Biophys. J. 89 (3), 1965-1977] provided experimental evidence supporting this hypothesis. Herein, we report further studies probing the differential heterocyclic substrate recognition properties of the aspartate 143 mutant TGTs. Our results are consistent with one of the mutants exhibiting an inversion of substrate recognition preference (xanthine vs guanine) relative to that of the wild type, as evidenced by Km values. This confirms the key role of aspartate 143 in maintaining the anticodon identities of the queuine-containing tRNAs and suggests that TGT mutants could be developed that would alter the tRNA wobble base base pairing properties.

Structural and thermodynamic features of spiroiminodihydantoin damaged DNA duplexes

Oxidation of guanine or 8-oxo-7,8-dihydroguanine can produce spiroiminodihydantoin (Sp) R and S stereoisomers. Both in vitro and in vivo experiments have shown that the Sp stereoisomers are highly mutagenic, causing G --> C and G --> T transversion mutations. Therefore, they are of interest as potential endogenous cancer causing lesions. However, their structural properties in DNA duplexes remain to be elucidated. We have employed computational methods to study the Sp lesions in 11-mer DNA duplexes with A, C, G, and T partners. Molecular dynamics simulations have been carried out to obtain ensembles of structures, and the trajectories were employed to analyze the structures and compute free energies. The structural and thermodynamic analyses reveal that the Sp stereoisomers energetically favor positioning in the B-DNA major groove, with minor groove conformers also low energy in some cases, depending on the partner base. The R and S stereoisomers adopt opposite orientations with respect to the 5' to 3' direction of the modified strand. Both syn and anti glycosidic bond conformations are energetically feasible, with partner base and stereochemistry determining the preference. The lesions adversely impact base stacking and Watson-Crick hydrogen bonding interactions in the duplex, and cause groove widening. The chemical nature of the partner base determines specific hydrogen bonding and stacking properties of the damaged duplexes. The structural characteristics may relate to observed mutagenic properties of the Sp stereoisomers, including possible stereoisomer-dependent differences.

Combining docking and molecular dynamic simulations in drug design

A rational approach is needed to maximize the chances of finding new drugs, and to exploit the opportunities of potential new drug targets emerging from genomic and proteomic initiatives, and from the large libraries of small compounds now readily available through combinatorial chemistry. Despite a shaky early history, computer-aided drug design techniques can now be effective in reducing costs and speeding up drug discovery. This happy outcome results from development of more accurate and reliable algorithms, use of more thoughtfully planned strategies to apply them, and greatly increased computer power to allow studies with the necessary reliability to be performed. Our review focuses on applications and protocols, with the main emphasis on critical analysis of recent studies where docking calculations and molecular dynamics (MD) simulations were combined to dock small molecules into protein receptors. We highlight successes to demonstrate what is possible now, but also point out drawbacks and future directions. The review is structured to lead the reader from the simpler to more compute-intensive methods. Thus, while inexpensive and fast docking algorithms can be used to scan large compound libraries and reduce their size, more accurate but expensive MD simulations can be applied when a few selected ligand candidates remain. MD simulations can be used: during the preparation of the protein receptor before docking, to optimize its structure and account for protein flexibility; for the refinement of docked complexes, to include solvent effects and account for induced fit; to calculate binding free energies, to provide an accurate ranking of the potential ligands; and in the latest developments, during the docking process itself to find the binding site and correctly dock the ligand a priori.

T315I-mutated Bcr-Abl in chronic myeloid leukemia and imatinib: insights from a computational study

The early stage of chronic myeloid leukemia is triggered by the tyrosine kinase Bcr-Abl. Imatinib mesylate, a selective inhibitor of Bcr-Abl, has been successful in chronic myeloid leukemia clinical trials, but short-lived remissions are usually observed in blast crisis patients. Sequencing of the BCR-ABL gene in relapsed patients revealed a set of mutants that mediate drug resistance. Previously reported work postulated that the missense T315I mutation both alters the three-dimensional structure of the protein binding site, thus decreasing the protein sensitivity for the drug, and does not feature a fundamental hydrogen bond that is critical for binding with imatinib. These speculations, however, were not supported by investigations at the molecular modeling level. Here, we present the results obtained from the application of molecular dynamics simulations to the study of the interactions between T315I Bcr-Abl and imatinib. For the first time, we show that, with respect to the wild-type system, the absence of the supposedly critical H-bond is not the only cause for the failure of receptor inhibition by imatinib, but also a plethora of other protein/drug interactions are drastically and unfavorably changed in the mutant protein.

Loss of G-A base pairs is insufficient for achieving a large opening of U4 snRNA K-turn motif

Upon binding to the 15.5K protein, two tandem-sheared G–A base pairs are formed in the internal loop of the kink-turn motif of U4 snRNA (Kt-U4). We have reported that the folding of Kt-U4 is assisted by protein binding. Unstable interactions that contribute to a large opening of the free RNA (‘k–e motion’) were identified using locally enhanced sampling molecular dynamics simulations, results that agree with experiments. A detailed analysis of the simulations reveals that the k–e motion in Kt-U4 is triggered both by loss of G–A base pairs in the internal loop and backbone flexibility in the stems. Essential dynamics show that the loss of G–A base pairs is correlated along the first mode but anti-correlated along the third mode with the k–e motion. Moreover, when enhanced sampling was confined to the internal loop, the RNA adopted an alternative conformation characterized by a sharper kink, opening of G–A base pairs and modified stacking interactions. Thus, loss of G–A base pairs is insufficient for achieving a large opening of the free RNA. These findings, supported by previously published RNA structure probing experiments, suggest that G–A base pair formation occurs upon protein binding, thereby stabilizing a selective orientation of the stems.

A binding mechanism in protein-nucleotide interactions: implication for U1A RNA binding

We present a close electronic view of the protein-base interface for the N-terminal domain of the human protein U1A. Combining accurate mixed quantum mechanics/molecular mechanics techniques and protein structure prediction methods, we provide a detailed electronic structure description of the protein-RNA stacking interactions. Our analysis indicates the evolution of the protein structure optimizing the interaction between Asp-92 and the RNA bases. The results show a direct coupling of the C-terminal tail and Asp-92, providing a direct rationalization of the experimentally determined role of the C-terminal domain in RNA binding. Here, we propose a mechanism where a protein side chain, with a delocalized electronic pi system, assists in the nucleotide binding. The binding mechanism involves a short-range interaction of the side chain with the nucleotide base and an electronic long-range interaction through a sandwich-stacking motif. The structural motif of the binding mechanism is observed in similar protein-RNA interactions and in various protein-ATP-binding sites.

Stochastic roadmap simulation: an efficient representation and algorithm for analyzing molecular motion

Classic molecular motion simulation techniques, such as Monte Carlo (MC) simulation, generate motion pathways one at a time and spend most of their time in the local minima of the energy landscape defined over a molecular conformation space. Their high computational cost prevents them from being used to compute ensemble properties (properties requiring the analysis of many pathways). This paper introduces stochastic roadmap simulation (SRS) as a new computational approach for exploring the kinetics of molecular motion by simultaneously examining multiple pathways. These pathways are compactly encoded in a graph, which is constructed by sampling a molecular conformation space at random. This computation, which does not trace any particular pathway explicitly, circumvents the local-minima problem. Each edge in the graph represents a potential transition of the molecule and is associated with a probability indicating the likelihood of this transition. By viewing the graph as a Markov chain, ensemble properties can be efficiently computed over the entire molecular energy landscape. Furthermore, SRS converges to the same distribution as MC simulation. SRS is applied to two biological problems: computing the probability of folding, an important order parameter that measures the "kinetic distance" of a protein's conformation from its native state; and estimating the expected time to escape from a ligand-protein binding site. Comparison with MC simulations on protein folding shows that SRS produces arguably more accurate results, while reducing computation time by several orders of magnitude. Computational studies on ligand-protein binding also demonstrate SRS as a promising approach to study ligand-protein interactions.

Free energy simulations and MM-PBSA analyses on the affinity and specificity of steroid binding to antiestradiol antibody

Antiestradiol antibody 57-2 binds 17beta-estradiol (E2) with moderately high affinity (K(a) = 5 x 10(8) M(-1)). The structurally related natural estrogens estrone and estriol as well synthetic 17-deoxy-estradiol and 17alpha-estradiol are bound to the antibody with 3.7-4.9 kcal mol(-1) lower binding free energies than E2. Free energy perturbation (FEP) simulations and the molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method were applied to investigate the factors responsible for the relatively low cross-reactivity of the antibody with these four steroids, differing from E2 by the substituents of the steroid D-ring. In addition, computational alanine scanning of the binding site residues was carried out with the MM-PBSA method. Both the FEP and MM-PBSA methods reproduced the experimental relative affinities of the five steroids in good agreement with experiment. On the basis of FEP simulations, the number of hydrogen bonds formed between the antibody and steroids, which varied from 0 to 3 in the steroids studied, determined directly the magnitude of the steroid-antibody interaction free energies. One hydrogen bond was calculated to contribute about 3 kcal mol(-1) to the interaction energy. Because the relative binding free energies of estrone (two antibody-steroid hydrogen bonds), estriol (three hydrogen bonds), 17-deoxy-estradiol (no hydrogen bonds), and 17alpha-estradiol (two hydrogen bonds) are close to each other and clearly lower than that of E2 (three hydrogen bonds), the water-steroid interactions lost upon binding to the antibody make an important contribution to the binding free energies. The MM-PBSA calculations showed that the binding of steroids to the antiestradiol antibody is driven by van der Waals interactions, whereas specificity is solely due to electrostatic interactions. In addition, binding of steroids to the antiestradiol antibody 57-2 was compared to the binding to the antiprogesterone antibody DB3 and antitestosterone antibody 3-C4F5, studied earlier with the MM-PBSA method.

Molecular recognition and docking algorithms

Molecular docking is an invaluable tool in modern drug discovery. This review focuses on methodological developments relevant to the field of molecular docking. The forces important in molecular recognition are reviewed and followed by a discussion of how different scoring functions account for these forces. More recent applications of computational chemistry tools involve library design and database screening. Last, we summarize several critical methodological issues that must be addressed in future developments.

Phenotypic expression of the C282Y/Q283P compound heterozygosity in HFE and molecular modeling of the Q283P mutation effect

In Caucasians, from 4 to 35% of hereditary hemochromatosis (HH) patients carry a least one chromosome without a common assigned HFE mutation (i.e., C282Y, H63D, and S65C). We have undertaken a D-HPLC scanning of the HFE coding region in such patients in order to identify uncommon mutations liable to explain their high transferrin saturation level. Twenty HH patients from Brittany carrying at least one chromosome without an assigned mutation were selected on the basis of a transferrin saturation level with the following threshold: > or = 60% in men and > or = 50% in women, in the absence of other known causes of iron disorders. This strategy allowed us to detect a heterozygous sequence variant in exon 4 of the HFE gene from one individual who was also heterozygous for C282Y. Subsequent DNA sequencing analysis identified an adenine to cytosine transversion at position 848 which changes amino acid 283 from glutamine to proline (Q283P). Family study revealed a clear association between the C282Y/Q283P compound heterozygote genotype and the development of HH. Molecular modeling studies are in favor of a destabilizing effect of the Q283P mutation on the tertiary structure of the HFE protein. This is the first report of a natural protein variant describing the introduction of a proline in a central beta-strand position. Our approach may have practical implications in screening strategies for hereditary hemochromatosis, molecular diagnosis, and HFE structure-function relationships.

Assembly of the U1 snRNP involves interactions with the backbone of the terminal stem of U1 snRNA

Nucleotide analog interference mapping (NAIM) is a powerful method for identifying RNA functional groups involved in protein-RNA interactions. We examined particles assembled on modified U1 small nuclear RNAs (snRNAs) in vitro and detected two categories of interferences. The first class affects the stability of two higher-order complexes and comprises changes in two adenosines, A65 and A70, in the loop region previously identified as the binding site for the U1 small nuclear ribonucleoprotein (snRNP)-specific U1A protein. Addition of an exocyclic amine to position 2 of A65 interferes strongly with protein binding, whereas removal or modification of the exocyclic amine at position 6 makes little difference. Modifications of A70 exhibit the opposite effects: Additions at position 2 are permitted, but modification of the exocyclic amine at position 6 significantly inhibits protein binding. These interactions, critical for U1A-U1 snRNA recognition in the context of in vitro snRNP assembly, are consistent with previous structural studies of the isolated protein with the RNA hairpin containing the U1A binding site. The second category of interferences affects all partially assembled U1-protein complexes by decreasing the stability of Sm core protein associations. Interestingly, most strong interferences occur at phosphates in the terminal stem-loop region of U1, rather than in the Sm binding site. These data argue that interactions with the phosphate backbone of the terminal stem loop are essential for the stable association of Sm core proteins with the U1 snRNA. We suggest that the stem loop of all Sm snRNAs may act as a clamp to hold the ring of Sm proteins in place.

Analysis of the binding energies of testosterone, 5alpha-dihydrotestosterone, androstenedione and dehydroepiandrosterone sulfate with an antitestosterone antibody

Molecular dynamics simulations and molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) free energy calculations were used to study the binding of testosterone (TES), 5alpha-dihydrotestosterone (5ADHT), androstenedione (AND), and dehydroepiandrosterone sulfate (DHEAS) to the monoclonal antitestosterone antibody 3-C(4)F(5). The relative binding free energy of TES and AND was also calculated with free energy perturbation (FEP) simulations. The antibody 3-C(4)F(5) has a relatively high affinity (3 x 10(8) M(-1)) and on overall good binding profile for testosterone but its cross-reactivity with DHEAS has been the main reason for the failure to use this antibody in clinical immunoassays. The relative binding free energies obtained with the MM-PBSA method were 1.5 kcal/mol for 5ADHT, 3.8 kcal/mol for AND, and 4.3 kcal/mol for DHEAS, as compared to TES. When a water molecule of the ligand binding site, observed in the antibody-TES crystal structure, was explicitly included in MM-PBSA calculations, the relative binding energies were 3.4, 4.9, and 5.4 kcal/mol for 5ADHT, AND, and DHEAS, respectively. The calculated numbers are in correct order but larger than the corresponding experimental energies of 1.3, 1.5, and 2.6 kcal/mol, respectively. The fact that the MM-PBSA method reproduced the relative binding free energies of DHEAS, a steroid having a negatively charged sulfate group, and the neutrally charged TES, 5ADHT, and AND in satisfactory agreement with experiment shows the robustness of the method in predicting relative binding affinities. The 800-ps FEP simulations predicted that the antibody 3-C(4)F(5) binds TES 1.3 kcal/mol tighter than AND. Computational mutagenesis of selected amino acid residues of the ligand binding site revealed that the lower affinities of AND and DHEAS as compared to TES are due to a combined effect of several residues, each contributing a small fraction to the tighter binding of TES. An exception to this is Tyr99H, whose mutation to Ala lowered the binding of DHEAS 0.7 kcal/mol more than the binding of TES. This is probably due to the hydrogen bonding interaction formed between the OH group of Tyr99H and the sulfate group of DHEAS. Computational mutagensis data also showed that the affinity of the steroids to the antitestosterone antibody 3-C(4)F(5) would be enhanced if Trp47H were repositioned so that it would make more extensive contacts with the bound ligands. In addition, the binding of steroids to antitestosterone, antiprogesterone, and antiestradiol antibodies is discussed.

Molecular dynamics simulation studies of induced fit and conformational capture in U1A-RNA binding: do molecular substates code for specificity?

Molecular dynamics (MD) simulations on stem loop 2 of U1 small nuclear RNA and a construct of the U1A protein were carried out to obtain predictions of the structures for the unbound forms in solution and to elucidate dynamical aspects of induced fit upon binding. A crystal structure of the complex between the U1A protein and stem loop 2 RNA and an NMR structure for the uncomplexed form of the U1A protein are available from Oubridge et al. (Nature, 1994, Vol. 372, pp. 432-438) and Avis et al. (Journal of Molecular Biology, 1996, Vol. 257, pp. 398-411), respectively. As a consequence, U1A-RNA binding is a particularly attractive case for investigations of induced fit in protein-nucleic acid complexation. When combined with the available structural data, the results from simulations indicate that structural adaptation of U1A protein and RNA define distinct mechanisms for induced fit. For the protein, the calculations indicate that induced fit upon binding involves a non-native thermodynamic substate in which the structure is preorganized for binding. In contrast, induced fit of the RNA involves a distortion of the native structure in solution to an unstable form. However, the RNA solution structures predicted from simulation show evidence that structures in which groups of bases are favorably oriented for binding the U1A protein are thermally accessible. These results, which quantify with computational modeling recent proposals on induced fit and conformational capture by Leuillot and Varani (Biochemistry, 2001, Vol. 40, pp. 7947-7956) and by Williamson (Nature Structural Biology, 2000, Vol. 7, pp. 834-837) suggest an important role for intrinsic molecular architecture and substates other than the native form in the specificity of protein-RNA interactions.

Simulations of nucleic acids and their complexes

Recent years have seen considerable progress in simulations of nucleic acids. Improvements in force fields, simulation techniques and protocols, and increasing computer power have all contributed to making nanosecond-scale simulations of both DNA and RNA commonplace. The results are already helping to explain how nucleic acids respond to their environment and to their base sequence and to reveal the factors underlying recognition processes by probing biologically important nucleic acid-protein interactions and medically important nucleic acid-drug complexation. This Account summarizes methodological progress and applications of molecular dynamics to nucleic acids over the past few years and tries to identify remaining challenges.

Free energy simulations come of age: protein-ligand recognition

In recent years, molecular dynamics simulations of biomolecular free energy differences have benefited from significant methodological advances and increased computer power. Applications to molecular recognition provide an understanding of the interactions involved that goes beyond, and is an important complement to, experimental studies. Poisson-Boltzmann electrostatic models provide a faster and simpler free energy method in cases where electrostatic interactions are important. We illustrate both molecular dynamics and Poisson-Boltzmann methods with a detailed study of amino acid recognition by aspartyl-tRNA synthetase, whose specificity is important for maintaining the integrity of the genetic code.

Structural clusters of evolutionary trace residues are statistically significant and common in proteins

Given the massive increase in the number of new sequences and structures, a critical problem is how to integrate these raw data into meaningful biological information. One approach, the Evolutionary Trace, or ET, uses phylogenetic information to rank the residues in a protein sequence by evolutionary importance and then maps those ranked at the top onto a representative structure. If these residues form structural clusters, they can identify functional surfaces such as those involved in molecular recognition. Now that a number of examples have shown that ET can identify binding sites and focus mutational studies on their relevant functional determinants, we ask whether the method can be improved so as to be applicable on a large scale. To address this question, we introduce a new treatment of gaps resulting from insertions and deletions, which streamlines the selection of sequences used as input. We also introduce objective statistics to assess the significance of the total number of clusters and of the size of the largest one. As a result of the novel treatment of gaps, ET performance improves measurably. We find evolutionarily privileged clusters that are significant at the 5% level in 45 out of 46 (98%) proteins drawn from a variety of structural classes and biological functions. In 37 of the 38 proteins for which a protein-ligand complex is available, the dominant cluster contacts the ligand. We conclude that spatial clustering of evolutionarily important residues is a general phenomenon, consistent with the cooperative nature of residues that determine structure and function. In practice, these results suggest that ET can be applied on a large scale to identify functional sites in a significant fraction of the structures in the protein databank (PDB). This approach to combining raw sequences and structure to obtain detailed insights into the molecular basis of function should prove valuable in the context of the Structural Genomics Initiative.

Computational alanine scanning of the 1:1 human growth hormone-receptor complex

The MM-PBSA (Molecular Mechanics-Poisson-Boltzmann surface area) method was applied to the human Growth Hormone (hGH) complexed with its receptor to assess both the validity and the limitations of the computational alanine scanning approach. A 400-ps dynamical trajectory of the fully solvated complex was simulated at 300 K in a 101 A x 81 A x 107 A water box using periodic boundary conditions. Long-range electrostatic interactions were treated with the particle mesh Ewald (PME) summation method. Equally spaced snapshots along the trajectory were chosen to compute the binding free energy using a continuum solvation model to calculate the electrostatic desolvation free energy and a solvent-accessible surface area approach to treat the nonpolar solvation free energy. Computational alanine scanning was performed on the same set of snapshots by mutating the residues in the structural epitope of the hormone and the receptor to alanine and recomputing the deltaGbinding. To further investigate a particular structure, a 200-ps dynamical trajectory of an R43A hormone-receptor complex was simulated. By postprocessing a single trajectory of the wild-type complex, the average unsigned error of our calculated deltadeltaGbinding is approximately1 kcal/mol for the alanine mutations of hydrophobic residues and polar/charged residues without buried salt bridges. When residues involved in buried salt bridges are mutated to alanine, it is demonstrated that a separate trajectory of the alanine mutant complex can lead to reasonable agreement with experimental results. Our approach can be extended to rapid screening of a variety of possible modifications to binding sites.

Calculations of free-energy contributions to protein-RNA complex stabilization

The problem of calculating binding affinities of protein-RNA complexes is addressed by analyzing a computational strategy of modeling electrostatic free energies based on a nonlinear Poisson-Boltzmann (NLPB) model and linear response approximation (LRA). The underlying idea is to treat binding as a two-step process. Solutions to the NLPB equation calculate free energies arising from electronic polarizability and the LRA is constructed from molecular dynamics simulations to model reorganization free energies due to conformational transitions. By implementing a consistency condition of requiring the NLPB model to reproduce the solute-solvent free-energy transitions determined by the LRA, a "macromolecule dielectric constant" (epsilon(m)) for treating reorganization is obtained. The applicability of this hybrid approach was evaluated by calculating the absolute free energy of binding and free-energy changes for amino acid substitutions in the complex between the U1A spliceosomal protein and its cognate RNA hairpin. Depending on the residue substitution, epsilon(m) varied from 3 to 18, and reflected dipolar reorientation not included in the polarization modeled by epsilon(m) = 2. Although the changes in binding affinities from substitutions modeled strictly at the implicit level by the NLPB equation with epsilon(m) = 4 reproduced the experimental values with good overall agreement, substitutions problematic to this simple treatment showed significant improvement when solved by the NLPB-LRA approach.

Energetic analysis of binding of progesterone and 5 beta-androstane-3,17-dione to anti-progesterone antibody DB3 using molecular dynamics and free energy calculations

Molecular dynamics simulations and molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) free energy calculations were used to study the energetics of the binding of progesterone (PRG) and 5 beta-androstane-3,17-dione (5AD) to anti-PRG antibody DB3. Although the two steroids bind to DB3 in different orientations, their binding affinities are of the same magnitude, 1 nM for PRG and 8 nM for 5AD. The calculated relative binding free energy of the steroids, 8.8 kJ/mol, is in fair agreement with the experimental energy, 5.4 kJ/mol. In addition, computational alanine scanning was applied to study the role of selected amino acid residues of the ligand-binding site on the steroid cross-reactivity. The electrostatic and van der Waals components of the total binding free energies were found to favour more the binding of PRG, whereas solvation energies were more favourable for the binding of 5AD. The differences in the free energy components are due to the binding of the A rings of the steroids to different binding pockets: PRG is bound to a pocket in which electrostatic antibody-steroid interactions are dominating, whereas 5AD is bound to a pocket in which van der Waals and hydrophobic interactions dominate.

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

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

Biomolecular simulations: recent developments in force fields, simulations of enzyme catalysis, protein-ligand, protein-protein, and protein-nucleic acid noncovalent interactions

Computer modeling has been developed and widely applied in studying molecules of biological interest. The force field is the cornerstone of computer simulations, and many force fields have been developed and successfully applied in these simulations. Two interesting areas are (a) studying enzyme catalytic mechanisms using a combination of quantum mechanics and molecular mechanics, and (b) studying macromolecular dynamics and interactions using molecular dynamics (MD) and free energy (FE) calculation methods. Enzyme catalysis involves forming and breaking of covalent bonds and requires the use of quantum mechanics. Noncovalent interactions appear ubiquitously in biology, but here we confine ourselves to review only noncovalent interactions between protein and protein, protein and ligand, and protein and nucleic acids.