Molecular Dynamics Simulations of the Hypoxia-Inducible Factor PAS-B Domain Confirm That Internally Bound Water Molecules Function To Stabilize the Protein Core for Ligand Binding.

Hypoxia-inducible factors (HIFs) are heterodimeric transcription factors that are induced in many cancer types, promoting tumor-forming pathways. The HIF-2α PAS-B domain of the human HIF-2 complex contains a large buried hydrophilic cavity that is filled with eight ordered water molecules, making the domain a promising drug target. Recent work has led to the hypothesis that these buried water molecules act as an internal scaffold that stabilizes the buried cavity, holding it open to ease ligand binding. Here we show that extensive molecular dynamics simulations of the HIF-2α PAS-B domain strongly support this hypothesis. If water is prevented from entering the buried cavity, several amino acid side chains lining the cavity undergo a major conformational change that shrinks the cavity volume by an average of ∼45% relative to the hydrated cavity. The structural change within the cavity is due mainly to side chain rearrangements; the absence of bound water has only minor effects on the protein backbone conformation. Moreover, when the cavity lacks internally bound water, the side chain of Tyr307 rotates such that it blocks a major route for entry of the ligand into the cavity, greatly raising the free energy barrier for binding of an artificial ligand compared to when the cavity contains water. Finally, there is no difference in global protein stability of PAS-B when the cavity lacks water, as indicated by the same degree of protein flexibility and identical predicted melting temperatures. These findings help to explain further how water stabilizes large internal cavities, a relatively underexamined feature of proteins.

Predicting drug activity against cancer cells by random forest models based on minimal genomic information and chemical properties.

A key goal of precision medicine is predicting the best drug therapy for a specific patient from genomic information. In oncology, cancers that appear similar pathologically can vary greatly in how they respond to the same drug. Fortunately, data from high-throughput screening programs often reveal important relationships between genomic variability of cancer cells and their response to drugs. Nevertheless, many current computational methods to predict compound activity against cancer cells require large quantities of genomic, epigenomic, and additional cellular data to develop and to apply. Here we integrate recent screening data and machine learning to train classification models that predict the activity/inactivity of compounds against cancer cells based on the mutational status of only 145 oncogenes and a set of compound structural descriptors. Using IC50 values of 1 µM as activity cutoffs, our predictive models have sensitivities of 87%, specificities of 87%, and yield an area under the receiver operating characteristic curve equal to 0.94. We also develop regression models to predict log(IC50) values of compounds for cancer cells; the models achieve a Pearson correlation coefficient of 0.86 for cross-validation and up to 0.65-0.73 against blind test sets. Predictive performance remains strong when as few as 50 oncogenes are included. Finally, even when 40% of experimental IC50 values are missing from screening data, they can be imputed with sufficient reliability that classification accuracy is not diminished. The presented models are fast to generate and may serve as easily implemented screening tools for personalized oncology medicine, drug repurposing, and drug discovery.

Complete Coupled Binding-Folding Pathway of the Intrinsically Disordered Transcription Factor Protein Brinker Revealed by Molecular Dynamics Simulations and Markov State Modeling.

Intrinsically disordered proteins (IDPs) make up a large class of proteins that lack stable structures in solution, existing instead as dynamic conformational ensembles. To perform their biological functions, many IDPs bind to other proteins or nucleic acids. Although IDPs are unstructured in solution, when they interact with binding partners, they fold into defined three-dimensional structures via coupled binding-folding processes. Because they frequently underlie IDP function, the mechanisms of this coupled binding-folding process are of great interest. However, given the flexibility inherent to IDPs and the sparse populations of intermediate states, it is difficult to reveal binding-folding pathways at atomic resolution using experimental methods. Computer simulations are another tool for studying these pathways at high resolution. Accordingly, we have applied 40 µs of unbiased molecular dynamics simulations and Markov state modeling to map the complete binding-folding pathway of a model IDP, the 59-residue C-terminal portion of the DNA binding domain of Drosophila melanogaster transcription factor Brinker (BrkDBD). Our modeling indicates that BrkDBD binds to its cognate DNA and folds in ∼50 µs by an induced fit mechanism, acquiring most of its stable secondary and tertiary structure only after it reaches the final binding site on the DNA. The protein follows numerous pathways en route to its bound and folded conformation, occasionally becoming stuck in kinetic traps. Each binding-folding pathway involves weakly bound, increasingly folded intermediate states located at different sites on the DNA surface. These findings agree with experimental data and provide additional insight into the BrkDBD folding mechanism and kinetics.

Molecular simulations and Markov state modeling reveal the structural diversity and dynamics of a theophylline-binding RNA aptamer in its unbound state.

RNA aptamers are oligonucleotides that bind with high specificity and affinity to target ligands. In the absence of bound ligand, secondary structures of RNA aptamers are generally stable, but single-stranded and loop regions, including ligand binding sites, lack defined structures and exist as ensembles of conformations. For example, the well-characterized theophylline-binding aptamer forms a highly stable binding site when bound to theophylline, but the binding site is unstable and disordered when theophylline is absent. Experimental methods have not revealed at atomic resolution the conformations that the theophylline aptamer explores in its unbound state. Consequently, in the present study we applied 21 microseconds of molecular dynamics simulations to structurally characterize the ensemble of conformations that the aptamer adopts in the absence of theophylline. Moreover, we apply Markov state modeling to predict the kinetics of transitions between unbound conformational states. Our simulation results agree with experimental observations that the theophylline binding site is found in many distinct binding-incompetent states and show that these states lack a binding pocket that can accommodate theophylline. The binding-incompetent states interconvert with binding-competent states through structural rearrangement of the binding site on the nanosecond to microsecond timescale. Moreover, we have simulated the complete theophylline binding pathway. Our binding simulations supplement prior experimental observations of slow theophylline binding kinetics by showing that the binding site must undergo a large conformational rearrangement after the aptamer and theophylline form an initial complex, most notably, a major rearrangement of the C27 base from a buried to solvent-exposed orientation. Theophylline appears to bind by a combination of conformational selection and induced fit mechanisms. Finally, our modeling indicates that when Mg2+ ions are present the population of binding-competent aptamer states increases more than twofold. This population change, rather than direct interactions between Mg2+ and theophylline, accounts for altered theophylline binding kinetics.

Structural consequences of mutations to the α-tocopherol transfer protein associated with the neurodegenerative disease ataxia with vitamin E deficiency.

The α-tocopherol transfer protein (α-TTP) is a liver protein that transfers α-tocopherol (vitamin E) to very-low-density lipoproteins (VLDLs). These VLDLs are then circulated throughout the body to maintain blood α-tocopherol levels. Mutations to the α-TTP gene are associated with ataxia with vitamin E deficiency, a disease characterized by peripheral nerve degeneration. In this study, molecular dynamics simulations of the E141K and R59W disease-associated mutants were performed. The mutants displayed disruptions in and around the ligand-binding pocket. Structural analysis and ligand docking to the mutant structures predicted a decreased affinity for α-tocopherol. To determine the detailed mechanism of the mutation-related changes, we developed a new tool called ContactWalker that analyzes contact differences between mutant and wild-type proteins and highlights pathways of altered contacts within the mutant proteins. Taken together, our findings are in agreement with experiment and suggest structural explanations for the weakened ability of the mutants to bind and carry α-tocopherol.

Identification of binding specificity-determining features in protein families.

We present a new approach for identifying features of ligand-protein binding interfaces that predict binding selectivity and demonstrate its effectiveness for predicting kinase inhibitor specificity. We analyzed a large set of human kinases and kinase inhibitors using clustering of experimentally determined inhibition constants (to define specificity classes of kinases and inhibitors) and virtual ligand docking (to extract structural and chemical features of the ligand-protein binding interfaces). We then used statistical methods to identify features characteristic of each class. Machine learning was employed to determine which combinations of characteristic features were predictive of class membership and to predict binding specificities and affinities of new compounds. Experiments showed predictions were 70% accurate. These results show that our method can automatically pinpoint on the three-dimensional binding interfaces pharmacophore-like features that act as "selectivity filters". The method is not restricted to kinases, requires no prior hypotheses about specific interactions, and can be applied to any protein families for which sets of structures and ligand binding data are available.

Dynameomics: a comprehensive database of protein dynamics.

The dynamic behavior of proteins is important for an understanding of their function and folding. We have performed molecular dynamics simulations of the native state and unfolding pathways of over 2000 protein/peptide systems (approximately 11,000 independent simulations) representing the majority of folds in globular proteins. These data are stored and organized using an innovative database approach, which can be mined to obtain both general and specific information about the dynamics and folding/unfolding of proteins, relevant subsets thereof, and individual proteins. Here we describe the project in general terms and the type of information contained in the database. Then we provide examples of mining the database for information relevant to protein folding, structure building, the effect of single-nucleotide polymorphisms, and drug design. The native state simulation data and corresponding analyses for the 100 most populated metafolds, together with related resources, are publicly accessible through http://www.dynameomics.org.

Principles of ligand binding within a completely buried cavity in HIF2alpha PAS-B.

Hypoxia-inducible factors (HIFs) are heterodimeric transcription factors responsible for the metazoan hypoxia response and promote tumor growth, metastasis, and resistance to cancer treatment. The C-terminal Per-ARNT-Sim (PAS) domain of HIF2alpha (HIF2alpha PAS-B) contains a preformed solvent-inaccessible cavity that binds artificial ligands that allosterically perturb the formation of the HIF heterodimer. To better understand how small molecules bind within this domain, we examined the structures and equilibrium and transition-state thermodynamics of HIF2alpha PAS-B with several artificial ligands using isothermal titration calorimetry, NMR exchange spectroscopy, and X-ray crystallography. Rapid association rates reveal that ligand binding is not dependent upon a slow conformational change in the protein to permit ligand access, despite the closed conformation observed in the NMR and crystal structures. Compensating enthalpic and entropic contributions to the thermodynamic barrier for ligand binding suggest a binding-competent transition state characterized by increased structural disorder. Finally, molecular dynamics simulations reveal conversion between open and closed conformations of the protein and pathways of ligand entry into the binding pocket.

The R46Q, R131Q and R154H polymorphs of human DNA glycosylase/beta-lyase hOgg1 severely distort the active site and DNA recognition site but do not cause unfolding.

Reactive oxygen species can cause widespread cellular damage, including base alterations and strand breaks in DNA. An array of DNA-repair enzymes constitutes an essential part of the line of defense that cells use against oxidative damage to the genome. A DNA glycosylase/beta-lyase enzyme, Ogg1, scavenges the genome for 8-oxoguanine, a major mutagenic DNA adduct induced by reactive oxygen species, and catalyzes its excision and subsequent cleavage of the DNA phosphate backbone. Several polymorphisms of Ogg1, including the single amino-acid substitutions R46Q, R131Q and R154H, are associated with a variety of human cancers. These three mutations have previously been characterized experimentally but no structural data have been published. We have performed multiple molecular dynamics simulations of R46Q, R131Q and R154H human Ogg1 to predict the structural and dynamical effects of the substitutions throughout the protein and specifically within the active site and substrate recognition site. None of the substitutions induced unfolding or global structural changes, instead their effects were confined principally to the active and recognition sites. Although the enzyme active site is located 18-21 A from the three investigated mutation sites, these mutations' structural effects propagate through space and cause a major change in the orientation and chemical environment of the active site side chains. This change appears likely to compromise the ability of the Lys 249 side chain to undergo a necessary deprotonation step prior to its nucleophilic attack of the DNA. The mutations also cause an expansion of the active site cavity, which may explain the experimentally observed decreases in substrate specificity.

Molecular basis for the structural instability of human DJ-1 induced by the L166P mutation associated with Parkinson's disease.

DJ-1 is a dimeric protein of unknown function in vivo. A mutation in the human DJ-1 gene causing substitution of proline for leucine at residue 166 (L166P) has been linked to early onset Parkinson's disease. Lack of structural stability has precluded experimental determination of atomic-resolution structures of the L166P DJ-1 polymorph. We have performed multiple molecular dynamics (MD) simulations ( approximately 1/3 mus) of the wild-type and L166P DJ-1 polymorph at physiological temperature to predict specific structural effects of the L166P substitution. L166P disrupted helices alpha1, alpha5, alpha6 and alpha8 with alpha8 undergoing particularly severe disruption. Secondary structural elements critical for protein stability and dimerization were significantly disrupted across the entire dimer interface, as were extended hydrophobic surfaces involved in dimer formation. Relative to wild-type DJ-1, L166P DJ-1 populated a broader ensemble of structures, many of which corresponded to distorted conformations. In a L166P dimer model the substitution significantly destabilized the dimer interface, interrupting >100 intermolecular contacts that are important for dimer formation. The L166P substitution also led to major perturbations in the region of a highly conserved cysteine residue (Cys-106) that participates in dimerization and that is critical for a proposed chaperone function of DJ-1. Cys-106 is located approximately 16 A from the substitution site, demonstrating that structural disruptions propagate throughout the whole protein. Furthermore, L166P DJ-1 showed a significant increase in hydrophobic surface area relative to wild-type protein, possibly explaining the tendency of the mutant protein to aggregate. These simulations provide details about specific structural disturbances throughout L166P DJ-1 that previous studies have not revealed.

Minimum sequence requirements for the binding of paromomycin to the rRNA decoding site A.

We have recently introduced a computational methodology that combines molecular dynamics (MD) simulations, free-energy calculations, and in vitro binding assays to predict the minimum RNA structural requirements for selective, high-affinity RNA binding to small-molecule ligands. Here, we show that this methodology can be applied to the conformationally flexible aminoglycoside antibiotic paromomycin. A RNA consisting of an 11-mer:10-mer duplex that contains one 16S ribosome RNA decoding A-site bound to paromomycin was simulated for 4 ns. The methodology predicts that the 11-mer:10-mer duplex binds to paromomycin with high affinity, whereas smaller RNA duplexes lose complex stability and the ability to bind paromomycin. The predicted high-affinity binding to paromomycin of the 11-mer:10-mer duplex was confirmed experimentally (EC(50) = 0.28 microM), as well as the inability of smaller complexes to bind. Our simulations show good agreement with experiment for dynamic and structural properties of the isolated A-site, including hydrogen-bonding networks and RNA structural rearrangements upon ligand binding. The results suggest that MD simulations can supplement in vitro methods as a tool for predicting minimum RNA-binding motifs for both small, rigid ligands, and large, flexible ligands when structural information is available.

Minimum sequence requirements for selective RNA-ligand binding: a molecular mechanics algorithm using molecular dynamics and free-energy techniques.

In vitro evolution techniques allow RNA molecules with unique functions to be developed. However, these techniques do not necessarily identify the simplest RNA structures for performing their functions. Determining the simplest RNA that binds to a particular ligand is currently limited to experimental protocols. Here, we introduce a molecular-mechanics based algorithm employing molecular dynamics simulations and free-energy methods to predict the minimum sequence requirements for selective ligand binding to RNA. The algorithm involves iteratively deleting nucleotides from an experimentally determined structure of an RNA-ligand complex, performing energy minimizations and molecular dynamics on each truncated structure, and assessing which truncations do not prohibit RNA binding to the ligand. The algorithm allows prediction of the effects of sequence modifications on RNA structural stability and ligand-binding energy. We have implemented the algorithm in the AMBER suite of programs, but it could be implemented in any molecular mechanics force field parameterized for nucleic acids. Test cases are presented to show the utility and accuracy of the methodology.

Identification of a 14mer RNA that recognizes and binds flavin mononucleotide with high affinity.

Aptamers are nucleic acids developed by in vitro evolution techniques that bind to specific ligands with high affinity and selectivity. Despite such high affinity and selectivity, however, in vitro evolution does not necessarily reveal the minimum structure of the nucleic acid required for selective ligand binding. Here, we show that a 35mer RNA aptamer for the cofactor flavin mononucleotide (FMN) identified by in vitro evolution can be computationally evolved to a mere 14mer structure containing the original binding pocket and eight scaffolding nucleotides while maintaining its ability to bind in vitro selectively to FMN. Using experimental and computational methodologies, we found that the 14mer binds with higher affinity to FMN (K(D) approximately 4 microM) than to flavin adenine dinucleotide (K(D) approximately 12 microM) or to riboflavin (K(D) approximately 13 microM),despite the negative charge of FMN. Different hydrogen-bond strengths resulting from differing ring-system electron densities associated with the aliphatic-chain charges appear to contribute to the selectivity observed for the binding of the 14mer to FMN and riboflavin. Our results suggest that high affinity and selectivity in ligand binding is not restricted to large RNAs, but can also be a property of extraordinarily short RNAs.

Unusually short RNA sequences: design of a 13-mer RNA that selectively binds and recognizes theophylline.

RNA plays critical roles in numerous biological processes and constitutes valuable therapeutic targets. RNA is significant not only for its roles in transmitting the genetic code but also for its enzymatic functions in ribozymes and in peptide bond formation in ribosomes. Recent studies have shown that RNAs containing as few as 22 nucleotides can be key elements in cellular functions. This suggests the possibility of using short RNAs as regulatory elements. Here, we show that ligand recognition and selectivity by RNA molecules can occur with only the presence of a binding pocket and as few as six additional scaffolding nucleotides holding the binding pocket in place. A 13-mer RNA truncation of a 33-mer aptamer for theophylline preserves the ability to bind to theophylline and to discriminate against the structurally similar compound caffeine. The truncated aptamer retains nearly all of the same structural elements in its binding site as those present in the original aptamer. This is the first demonstration of selective ligand binding by a 13-mer RNA.

Nrf2, not the estrogen receptor, mediates catechol estrogen-induced activation of the antioxidant responsive element.

The antioxidant responsive element (ARE) plays an important role in the gene expression of phase II detoxification enzymes, such as NAD(P)H:quinone oxidoreductase 1 (NQO1), and NF-E2-related factor2 (Nrf2) is the transcription factor for the ARE-driven genes. Interestingly, estrogen receptor (ER) was reported to increase NQO1 gene expression through the ARE. In this study, we investigated the role of ER and Nrf2 in ARE activation using IMR-32 cells and mouse primary astrocytes. Among tested estrogen-related compounds, only catechol estrogens (i.e. 4-hydroxyestradiol) activated the ARE. Since 4-hydroxyestradiol-induced ARE activation was not inhibited by either 17beta-estradiol or tamoxifen, and overexpression of ER-alpha decreased 4-hydroxyestradiol-induced ARE activation, ARE activation by catechol estrogen was independent of ER. Nrf2, however, was very important in the 4-hydroxyestradiol-induced ARE activation. 4-Hydroxyestradiol did not activate the ARE in Nrf2 knockout (-/-) primary astrocytes, but did activate the ARE when Nrf2 was transfected into Nrf2-/- astrocytes. In addition, dominant negative Nrf2 completely blocked 4-hydroxyestradiol-induced ARE activation in Nrf2+/+ astrocytes, and only 4-hydroxyestradiol induced Nrf2 nuclear translocation in IMR-32 cells. A selective phosphatidylinositol 3-kinase (PI3-kinase) inhibitor (LY294002) blocked 4-hydroxyestradiol-induced Nrf2 nuclear translocation and NQO1 activity induction in IMR-32 cells. Taken together, these observations suggest that 4-hydroxyestradiol activates the ARE by a PI3-kinase-Nrf2 dependent mechanism, not involving ER.