Search

Covalent inhibition offers many advantages over non-covalent inhibition, but covalent warhead reactivity must be carefully balanced to maintain potency while avoiding unwanted side effects. While warhead reactivities are commonly measured with assays, a computational model to predict warhead reactivities could be useful for several aspects of the covalent inhibitor design process. Studies have shown correlations between covalent warhead reactivities and quantum mechanic (QM) properties that describe important aspects of the covalent reaction mechanism. However, the models from these studies are often linear regression equations and can have limitations associated with their usage. Applications of machine learning (ML) models to predict covalent warhead reactivities with QM descriptors are not extensively seen in the literature. This study uses QM descriptors, calculated at different levels of theory, to train ML models to predict reactivities of covalent acrylamide warheads. The QM/ML models are compared with linear regression models built upon the same QM descriptors and with ML models trained on structure-based features like Morgan fingerprints and RDKit descriptors. Experiments show that the QM/ML models outperform the linear regression models and the structure-based ML models, and literature test sets demonstrate the power of the QM/ML models to predict reactivities of unseen acrylamide warhead scaffolds. Ultimately, these QM/ML models are effective, computationally feasible tools that can expedite the design of new covalent inhibitors., (© 2024. The Author(s), under exclusive licence to Springer Nature Switzerland AG.)

Stimulator of interferon genes (STING) is a dimeric transmembrane adapter protein that plays a key role in the human innate immune response to infection and has been therapeutically exploited for its antitumor activity. The activation of STING requires its high-order oligomerization, which could be induced by binding of the endogenous ligand, cGAMP, to the cytosolic ligand-binding domain. Here we report the discovery through functional screens of a class of compounds, named NVS-STGs, that activate human STING. Our cryo-EM structures show that NVS-STG2 induces the high-order oligomerization of human STING by binding to a pocket between the transmembrane domains of the neighboring STING dimers, effectively acting as a molecular glue. Our functional assays showed that NVS-STG2 could elicit potent STING-mediated immune responses in cells and antitumor activities in animal models., (© 2023. The Author(s).)

Disruption of the YAP-TEAD protein-protein interaction is an attractive therapeutic strategy in oncology to suppress tumor progression and cancer metastasis. YAP binds to TEAD at a large flat binding interface (∼3500 Å 2 ) devoid of a well-defined druggable pocket, so it has been difficult to design low-molecular-weight compounds to abrogate this protein-protein interaction directly. Recently, work by Furet and coworkers ( ChemMedChem 2022 , DOI: 10.1002/cmdc.202200303) reported the discovery of the first class of small molecules able to efficiently disrupt the transcriptional activity of TEAD by binding to a specific interaction site of the YAP-TEAD binding interface. Using high-throughput in silico docking, they identified a virtual screening hit from a hot spot derived from their previously rationally designed peptidic inhibitor. Structure-based drug design efforts led to the optimization of the hit compound into a potent lead candidate. Given advances in rapid high-throughput screening and rational approaches to peptidic ligand discovery for challenging targets, we analyzed the pharmacophore features involved in transferring from the peptidic to small-molecule inhibitor that could enable small-molecule discovery for such targets. Here, we show retrospectively that pharmacophore analysis augmented by solvation analysis of molecular dynamics trajectories can guide the designs, while binding free energy calculations provide greater insight into the binding conformation and energetics accompanying the association event. The computed binding free energy estimates agree well with experimental findings and offer useful insight into structural determinants that influence ligand binding to the TEAD interaction surface, even for such a shallow binding site. Taken together, our results demonstrates the utility of advanced in silico methods in structure-based design efforts for difficult-to-drug targets such as the YAP-TEAD transcription factor complex.

Accurate estimation of the p K a 's of cysteine residues in proteins could inform targeted approaches in hit discovery. The p K a of a targetable cysteine residue in a disease-related protein is an important physiochemical parameter in covalent drug discovery, as it influences the fraction of nucleophilic thiolate amenable to chemical protein modification. Traditional structure-based in silico tools are limited in their predictive accuracy of cysteine p K a 's relative to other titratable residues. Additionally, there are limited comprehensive benchmark assessments for cysteine p K a predictive tools. This raises the need for extensive assessment and evaluation of methods for cysteine p K a prediction. Here, we report the performance of several computational p K a methods, including single-structure and ensemble-based approaches, on a diverse test set of experimental cysteine p K a 's retrieved from the PKAD database. The dataset consisted of 16 wildtype and 10 mutant proteins with experimentally measured cysteine p K a values. Our results highlight that these methods are varied in their overall predictive accuracies. Among the test set of wildtype proteins evaluated, the best method (MOE) yielded a mean absolute error of 2.3 p K units, highlighting the need for improvement of existing p K a methods for accurate cysteine p K a estimation. Given the limited accuracy of these methods, further development is needed before these approaches can be routinely employed to drive design decisions in early drug discovery efforts.

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased "downhill" trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand-protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

Relative binding free-energy (RBFE) calculations are experiencing resurgence in the computer-aided drug design of novel small molecules due to performance gains allowed by cutting-edge molecular mechanic force fields and computer hardware. Application of RBFE to soluble proteins is becoming a routine, while recent studies outline necessary steps to successfully apply RBFE at the orthosteric site of membrane-embedded G-protein-coupled receptors (GPCRs). In this work, we apply RBFE to a congeneric series of antagonists that bind to a lipid-exposed, extra-helical site of the P2Y1 receptor. We find promising performance of RBFE, such that it may be applied in a predictive manner on drug discovery programs targeting lipid-exposed sites. Further, by the application of the microkinetic model, binding at a lipid-exposed site can be split into (1) membrane partitioning of the drug molecule followed by (2) binding at the extra-helical site. We find that RBFE can be applied to calculate the free energy of each step, allowing the uncoupling of observed binding free energy from the influence of membrane affinity. This protocol may be used to identify binding hot spots at extra-helical sites and guide drug discovery programs toward optimizing intrinsic activity at the target.

The serine protease factor XI (FXI) is a prominent drug target as it holds promise to deliver efficacious anticoagulation without an enhanced risk of major bleeds. Several efforts have been described targeting the active form of the enzyme, FXIa. Herein, we disclose our efforts to identify potent, selective, and orally bioavailable inhibitors of FXIa. Compound 1 , identified from a diverse library of internal serine protease inhibitors, was originally designed as a complement factor D inhibitor and exhibited submicromolar FXIa activity and an encouraging absorption, distribution, metabolism, and excretion (ADME) profile while being devoid of a peptidomimetic architecture. Optimization of interactions in the S1, S1β, and S1' pockets of FXIa through a combination of structure-based drug design and traditional medicinal chemistry led to the discovery of compound 23 with subnanomolar potency on FXIa, enhanced selectivity over other coagulation proteases, and a preclinical pharmacokinetics (PK) profile consistent with bid dosing in patients.

The anticancer agent indisulam inhibits cell proliferation by causing degradation of RBM39, an essential mRNA splicing factor. Indisulam promotes an interaction between RBM39 and the DCAF15 E3 ligase substrate receptor, leading to RBM39 ubiquitination and proteasome-mediated degradation. To delineate the precise mechanism by which indisulam mediates the DCAF15-RBM39 interaction, we solved the DCAF15-DDB1-DDA1-indisulam-RBM39(RRM2) complex structure to a resolution of 2.3 Å. DCAF15 has a distinct topology that embraces the RBM39(RRM2) domain largely via non-polar interactions, and indisulam binds between DCAF15 and RBM39(RRM2), coordinating additional interactions between the two proteins. Studies with RBM39 point mutants and indisulam analogs validated the structural model and defined the RBM39 α-helical degron motif. The degron is found only in RBM23 and RBM39, and only these proteins were detectably downregulated in indisulam-treated HCT116 cells. This work further explains how indisulam induces RBM39 degradation and defines the challenge of harnessing DCAF15 to degrade additional targets.

Recently much effort has been invested in using convolutional neural network (CNN) models trained on 3D structural images of protein-ligand complexes to distinguish binding from non-binding ligands for virtual screening. However, the dearth of reliable protein-ligand x-ray structures and binding affinity data has required the use of constructed datasets for the training and evaluation of CNN molecular recognition models. Here, we outline various sources of bias in one such widely-used dataset, the Directory of Useful Decoys: Enhanced (DUD-E). We have constructed and performed tests to investigate whether CNN models developed using DUD-E are properly learning the underlying physics of molecular recognition, as intended, or are instead learning biases inherent in the dataset itself. We find that superior enrichment efficiency in CNN models can be attributed to the analogue and decoy bias hidden in the DUD-E dataset rather than successful generalization of the pattern of protein-ligand interactions. Comparing additional deep learning models trained on PDBbind datasets, we found that their enrichment performances using DUD-E are not superior to the performance of the docking program AutoDock Vina. Together, these results suggest that biases that could be present in constructed datasets should be thoroughly evaluated before applying them to machine learning based methodology development., Competing Interests: The authors’ affiliation with Deep Waters, LLC (T.K., A.C.) and Novartis (C.J.D., V.H., and J.S.D.) does not alter our adherence to PLOS ONE policies on sharing data and materials.

A simple descriptor calculated from molecular dynamics simulations of the membrane partitioning event is found to correlate well with experimental measurements of passive membrane permeation from the high-throughput MDCK-LE assay using a data set of 49 drug-like molecules. This descriptor approximates the energy cost of translocation across the hydrophobic membrane core (flip-flop), which for many molecules limits permeability. Performance is found to be superior in comparison to calculated properties such as clogP, clogD, or polar surface area. Furthermore, the atomistic simulations provide a structural understanding of the partitioned drug-membrane complex, facilitating medicinal chemistry optimization of membrane permeability.

We present a simple approach to calculate the kinetic properties of lipid membrane crossing processes from biased molecular dynamics simulations. We demonstrate that by using biased simulations, one can obtain highly accurate kinetic information with significantly reduced computational time with respect to unbiased simulations. We describe how to conveniently calculate the transition rates to enter, cross, and exit the membrane in terms of the mean first passage times. To obtain free energy barriers and relaxation times from biased simulations only, we constructed Markov models using the dynamic histogram analysis method (DHAM). The permeability coefficients that are calculated from the relaxation times are found to correlate highly with experimentally evaluated values. We show that more generally, certain calculated kinetic properties linked to the crossing of the membrane layer (e.g., barrier height and barrier crossing rates) are good indicators of ordering drugs by permeability. Extending the analysis to a 2D Markov model provides a physical description of the membrane crossing mechanism.

The paper describes the SAR/SPR studies that led to the discovery of phenoxy cyclopropyl phenyl acetamide derivatives as potent and selective GPR119 agonists. Based on a cis cyclopropane scaffold discovered previously, phenyl acetamides such as compound 17 were found to have excellent GPR119 potency and improved physicochemical properties. Pharmacokinetic data of compound 17 in rat, dog and rhesus will be described. Compound 17 was suitable for QD dosing based on its predicted human half-life, and its projected human dose was much lower than that of the recently reported structurally-related benzyloxy compound 2. Compound 17 was selected as a tool compound candidate for NHP (Non-Human Primate) efficacy studies., (Copyright © 2017 Elsevier Ltd. All rights reserved.)

The membrane dipole potential (Ψ d ) constitutes one of three electrical potentials generated by cell membranes. Ψ d arises from the unfavorable parallel alignment of phospholipid and water dipoles, and varies in magnitude both longitudinally and laterally across the bilayer according to membrane composition and phospholipid packing density. In this work, we propose that dynamic counter-balancing between Ψ d and the transmembrane potential (ΔΨ m ) governs the conformational state transitions of voltage-gated ion channels. Ψ d consists of 1) static outer, and dynamic inner leaflet components (Ψ d(extra) and Ψ d(intra) , respectively); and 2) a transmembrane component (ΔΨ d(inner-outer) ), ariing from differences in intra- and extracellular leaflet composition. Ψ d(intra) , which transitions between high and low energy states (Ψ d(intra, high) and Ψ d(intra, low) ) as a function of channel conformation, is transduced by the pore domain. ΔΨ d(inner-outer) is transduced by the voltage-sensing (VS) domain in summation with ΔΨ m . Potentiation of voltage-gated ion channels is of interest for the treatment of cardiac, neuronal, and other disorders arising from inherited/acquired ion channel dysfunction. Potentiators are widely believed to alter the rates and voltage-dependencies of channel gating transitions by binding to pockets in the membrane-facing and other regions of ion channel targets. Here, we propose that potentiators alter Ψ d(intra) and/or Ψ d(extra) , thereby increasing or decreasing the energy barriers governing channel gating transitions. We used quantum mechanical and molecular dynamics (MD) simulations to predict the overall Ψ d -modulating effects of a series of published positive hERG potentiators partitioned into model DOPC bilayers. Our findings suggest a strong correlation between the magnitude of Ψ d -lowering and positive hERG potentiation across the series., (Copyright © 2016 Elsevier B.V. All rights reserved.)

Passive membrane permeation of small molecules is essential to achieve the required absorption, distribution, metabolism, and excretion (ADME) profiles of drug candidates, in particular intestinal absorption and transport across the blood-brain barrier. Computational investigations of this process typically involve either building QSAR models or performing free energy calculations of the permeation event. Although insightful, these methods rarely bridge the gap between computation and experiment in a quantitative manner, and identifying structural insights to apply toward the design of compounds with improved permeability can be difficult. In this work, we combine molecular dynamics simulations capturing the kinetic steps of permeation at the atomistic level with a dynamic mechanistic model describing permeation at the in vitro level, finding a high level of agreement with experimental permeation measurements. Calculation of the kinetic rate constants determining each step in the permeation event allows derivation of structure-kinetic relationships of permeation. We use these relationships to probe the structural determinants of membrane permeation, finding that the desolvation/loss of hydrogen bonding required to leave the membrane partitioned position controls the membrane flip-flop rate, whereas membrane partitioning determines the rate of leaving the membrane.

Cellular drug targets exist within networked function-generating systems whose constituent molecular species undergo dynamic interdependent non-equilibrium state transitions in response to specific perturbations (i.e.. inputs). Cellular phenotypic behaviors are manifested through the integrated behaviors of such networks. However, in vitro data are frequently measured and/or interpreted with empirical equilibrium or steady state models (e.g. Hill, Michaelis-Menten, Briggs-Haldane) relevant to isolated target populations. We propose that cells act as analog computers, "solving" sets of coupled "molecular differential equations" (i.e. represented by populations of interacting species)via "integration" of the dynamic state probability distributions among those populations. Disconnects between biochemical and functional/phenotypic assays (cellular/in vivo) may arise with targetcontaining systems that operate far from equilibrium, and/or when coupled contributions (including target-cognate partner binding and drug pharmacokinetics) are neglected in the analysis of biochemical results. The transformation of drug discovery from a trial-and-error endeavor to one based on reliable design criteria depends on improved understanding of the dynamic mechanisms powering cellular function/dysfunction at the systems level. Here, we address the general mechanisms of molecular and cellular function and pharmacological modulation thereof. We outline a first principles theory on the mechanisms by which free energy is stored and transduced into biological function, and by which biological function is modulated by drug-target binding. We propose that cellular function depends on dynamic counter-balanced molecular systems necessitated by the exponential behavior of molecular state transitions under non-equilibrium conditions, including positive versus negative mass action kinetics and solute-induced perturbations to the hydrogen bonds of solvating water versus kT., (Copyright© Bentham Science Publishers; For any queries, please email at epub@benthamscience.org.)

Ligand binding to membrane proteins may be significantly influenced by the interaction of ligands with the membrane. In particular, the microscopic ligand concentration within the membrane surface solvation layer may exceed that in bulk solvent, resulting in overestimation of the intrinsic protein-ligand binding contribution to the apparent/measured affinity. Using published binding data for a set of small molecules with the β2 adrenergic receptor, we demonstrate that deconvolution of membrane and protein binding contributions allows for improved structure-activity relationship analysis and structure-based drug design. Molecular dynamics simulations of ligand bound membrane protein complexes were used to validate binding poses, allowing analysis of key interactions and binding site solvation to develop structure-activity relationships of β2 ligand binding. The resulting relationships are consistent with intrinsic binding affinity (corrected for membrane interaction). The successful structure-based design of ligands targeting membrane proteins may require an assessment of membrane affinity to uncouple protein binding from membrane interactions.

Bruton's tyrosine kinase (BTK) is a Tec family kinase with a well-defined role in the B cell receptor (BCR) pathway. It has become an attractive kinase target for selective B cell inhibition and for the treatment of B cell related diseases. We report a series of compounds based on 8-amino-imidazo[1,5-a]pyrazine that are potent reversible BTK inhibitors with excellent kinase selectivity. Selectivity is achieved through specific interactions of the ligand with the kinase hinge and driven by aminopyridine hydrogen bondings with Ser538 and Asp539, and by hydrophobic interaction of trifluoropyridine in the back pocket. These interactions are evident in the X-ray crystal structure of the lead compounds 1 and 3 in the complex with the BTK enzyme. Our lead compounds show desirable PK profiles and efficacy in the preclinical rat collagen induced arthritis model.

We report herein the design and synthesis of a series of potent and selective GPR119 agonists. Our objective was to develop a GPR119 agonist with properties that were suitable for fixed-dose combination with a DPP4 inhibitor. Starting from a phenoxy analogue (1), medicinal chemistry efforts directed toward reducing half-life and increasing solubility led to the synthesis of a series of benzyloxy analogues. Compound 28 was chosen for further profiling because of its favorable physicochemical properties and excellent GPR119 potency across species. This compound exhibited a clean off-target profile in counterscreens and good in vivo efficacy in mouse oGTT.

Ilicicolin H is a broad spectrum antifungal agent showing sub micro g/mL MICs against Candida spp., Aspergillus fumigatus and Cryptococcus spp. It is a potent inhibitor (C50 2-3ng/mL) of the mitochondrial cytochrome bc1 reductase with over 1000-fold selectivity against rat liver cytochrome bc1 reductase. Structure-activity relationship of semisynthetic derivatives by chemical modification of ilicicolin H and its 19-hydroxy derivative produced by biotransformation have been described. Basic 4'-esters and moderately polar N- and O-alkyl derivatives retained antifungal and the cytochrome bc1 reductase activities. 4',19-Diacetate and 19-cyclopropyl acetate retained antifungal and enzyme activity and selectivity with over 20-fold improvement of plasma protein binding., (Copyright © 2013 Elsevier Ltd. All rights reserved.)

Ilicicolin H is a polyketide-nonribosomal peptide synthase (NRPS)-natural product isolated from Gliocadium roseum, which exhibits potent and broad spectrum antifungal activity, with sub-μg/mL MICs against Candida spp., Aspergillus fumigatus, and Cryptococcus spp. It showed a novel mode of action, potent inhibition (IC50 = 2-3 ng/mL) of the mitochondrial cytochrome bc1 reductase, and over 1000-fold selectivity relative to rat liver cytochrome bc1 reductase. Ilicicolin H exhibited in vivo efficacy in murine models of Candida albicans and Cryptococcus neoformans infections, but efficacy may have been limited by high plasma protein binding. Systematic structural modification of ilicicolin H was undertaken to understand the structural requirement for the antifungal activity. The details of the biological activity of ilicicolin H and structural modification of some of the key parts of the molecule and resulting activity of the derivatives are discussed. These data suggest that the β-keto group is critical for the antifungal activity.

Bombesin receptor subtype 3 (BRS-3) is an orphan G-protein coupled receptor expressed primarily in the hypothalamus which plays a role in the onset of both diabetes and obesity. We report herein our progress made towards identifying a potent, selective bombesin receptor subtype-3 (BRS-3) agonist related to the previously described MK-7725(1) Chobanian et al. (2012) that would prevent atropisomerization through the increase of steric bulk at the C-2 position. This would thereby make clinical development of this class of compounds more cost effective by inhibiting racemization which can occur over long periods of time at room/elevated temperature., (Copyright © 2012 Elsevier Ltd. All rights reserved.)

Random forest is currently considered one of the best QSAR methods available in terms of accuracy of prediction. However, it is computationally intensive. Naïve Bayes is a simple, robust classification method. The Laplacian-modified Naïve Bayes implementation is the preferred QSAR method in the widely used commercial chemoinformatics platform Pipeline Pilot. We made a comparison of the ability of Pipeline Pilot Naïve Bayes (PLPNB) and random forest to make accurate predictions on 18 large, diverse in-house QSAR data sets. These include on-target and ADME-related activities. These data sets were set up as classification problems with either binary or multicategory activities. We used a time-split method of dividing training and test sets, as we feel this is a realistic way of simulating prospective prediction. PLPNB is computationally efficient. However, random forest predictions are at least as good and in many cases significantly better than those of PLPNB on our data sets. PLPNB performs better with ECFP4 and ECFP6 descriptors, which are native to Pipeline Pilot, and more poorly with other descriptors we tried., (© 2012 American Chemical Society)

We report herein the discovery of benzodiazepine sulfonamide-based bombesin receptor subtype 3 (BRS-3) agonists and their unusual chirality. Starting from a high-throughput screening lead, we prepared a series of BRS-3 agonists with improved potency and pharmacokinetic properties, of which compound 8a caused mechanism-based, dose-dependent food intake reduction and body weight loss after oral dosing in diet-induced obese mice. This effort also led to the discovery of a novel family of chiral molecules originated from the conformationally constrained seven-membered diazepine ring.

Structural restraints provided by solid-state NMR measurements of the metarhodopsin II intermediate are combined with molecular dynamics simulations to help visualize structural changes in the light activation of rhodopsin. Since the timescale for the formation of the metarhodopsin II intermediate (>1 ms) is beyond that readily accessible by molecular dynamics, we use NMR distance restraints derived from 13C dipolar recoupling measurements to guide the simulations. The simulations yield a working model for how photoisomerization of the 11-cis retinylidene chromophore bound within the interior of rhodopsin is coupled to transmembrane helix motion and receptor activation. The mechanism of activation that emerges is that multiple switches on the extracellular (or intradiscal) side of rhodopsin trigger structural changes that converge to disrupt the ionic lock between helices H3 and H6 on the intracellular side of the receptor., (Copyright (c) 2009. Elsevier Ltd. All rights reserved.)

Rhodopsin is a highly specialized G protein-coupled receptor (GPCR) that is activated by the rapid photochemical isomerization of its covalently bound 11-cis-retinal chromophore. Using two-dimensional solid-state NMR spectroscopy, we defined the position of the retinal in the active metarhodopsin II intermediate. Distance constraints were obtained between amino acids in the retinal binding site and specific (13)C-labeled sites located on the beta-ionone ring, polyene chain, and Schiff base end of the retinal. We show that the retinal C20 methyl group rotates toward the second extracellular loop (EL2), which forms a cap on the retinal binding site in the inactive receptor. Despite the trajectory of the methyl group, we observed an increase in the C20-Gly(188) (EL2) distance consistent with an increase in separation between the retinal and EL2 upon activation. NMR distance constraints showed that the beta-ionone ring moves to a position between Met(207) and Phe(208) on transmembrane helix H5. Movement of the ring toward H5 was also reflected in increased separation between the Cepsilon carbons of Lys(296) (H7) and Met(44) (H1) and between Gly(121) (H3) and the retinal C18 methyl group. Helix-helix interactions involving the H3-H5 and H4-H5 interfaces were also found to change in the formation of metarhodopsin II reflecting increased retinal-protein interactions in the region of Glu(122) (H3) and His(211) (H5). We discuss the location of the retinal in metarhodopsin II and its interaction with sequence motifs, which are highly conserved across the pharmaceutically important class A GPCR family, with respect to the mechanism of receptor activation.

The second extracellular loop (EL2) of rhodopsin forms a cap over the binding site of its photoreactive 11-cis retinylidene chromophore. A crucial question has been whether EL2 forms a reversible gate that opens upon activation or acts as a rigid barrier. Distance measurements using solid-state (13)C NMR spectroscopy between the retinal chromophore and the beta4 strand of EL2 show that the loop is displaced from the retinal binding site upon activation, and there is a rearrangement in the hydrogen-bonding networks connecting EL2 with the extracellular ends of transmembrane helices H4, H5 and H6. NMR measurements further reveal that structural changes in EL2 are coupled to the motion of helix H5 and breaking of the ionic lock that regulates activation. These results provide a comprehensive view of how retinal isomerization triggers helix motion and activation in this prototypical G protein-coupled receptor.

FapydG is a common oxidative DNA lesion involving opening of the imidazole ring. It shares the same precursor as 8-oxodG and can be excised by the same enzymes as 8-oxodG. However, the loss of the aromatic imidazole in FapydG results in a reduction of the double bond character between C5 and N7, with an accompanying increase in conformational flexibility. Experimental characterization of FapydG is hampered by high reactivity, and thus it is desirable to investigate structural details through computer simulation. We show that the existing Amber force field parameters for FapydG do not reproduce X-ray structural data. We employed quantum mechanics energy profile calculations to derive new molecular mechanics parameters for the rotation of the dihedral angles in the eximidazole moiety. Using these parameters, all-atom simulations in explicit water reproduce the nonplanar conformation of cFapydG in the crystal structure of the complex with L. lactis glycosylase Fpg. We note that the nonplanar structure is stabilized by an acidic residue that is not present in most Fpg sequences. Simulations of the E-->S mutant, as present in E. coli, resulted in a more planar conformation, suggesting that the highly nonplanar form observed in the crystal structure may not have direct biological relevance for FapydG.

G protein-coupled receptors (GPCRs) belonging to class A contain several highly conserved (>90%) amino acids in their transmembrane helices. Results of mutational studies of these highly conserved residues suggest a common mechanism for locking GPCRs in an inactive conformation and for their subsequent activation upon ligand binding. Recently, a second set of sites in the transmembrane helices has been identified in which amino acids with small side chains, such as Gly, Ala, Ser, Thr, and Cys, are highly conserved (>90%) when considered as a group. These group-conserved residues have not been recognized as having essential structural or functional roles. To determine the role of group-conserved residues in the beta(2)-adrenergic receptor (beta(2)-AR), amino acid replacements guided by molecular modeling were carried out at key positions in transmembrane helices H2-H4. The most significant changes in receptor expression and activity were observed upon replacement of the amino acids Ser-161 and Ser-165 in H4. Substitution at these sites by larger residues lowered the expression and activity of the receptor but did not affect specific binding to the antagonist ligand dihydroalprenolol. A second site mutation, V114A, rescued the low expression of the S165V mutant. Substitution of other group-conserved residues in H2-H4 by larger amino acids lowered receptor activity in the order Ala-128, Ala-76, Ser-120, and Ala-78. Together these data provide comprehensive analysis of group-conserved residues in a class A GPCR and allow insights into the roles of these residues in GPCR structure and function.

The 36-residue helical subdomain of the villin headpiece, HP36, is one of the smallest cooperatively folded proteins, folding on the microsecond time scale. The domain is an extraordinarily popular model system for both experimental and computational studies of protein folding. The structure of HP36 has been determined using X-ray crystallography and NMR spectroscopy, with the resulting structures exhibiting differences in helix packing, van der Waals contacts, and hydrogen bonding. It is important to determine the solution structure of HP36 with as much accuracy as possible since this structure is widely used as a reference for simulations and experiments. We complement the existing data by using all-atom molecular dynamics simulations with explicit solvent to evaluate which of the experimental models is the better representation of HP36 in solution. After simulation for 50 ns initiated with the NMR structure, we observed that the protein spontaneously adopts structures with a backbone conformation, core packing, and C-capping motif on the third helix that are more consistent with the crystal structure. We also examined hydrogen bonding and side chain packing interactions between D44 and R55 and between F47 and R55, respectively, which were observed in the crystal structure but not in the NMR-based solution structure. Simulations showed large fluctuations in the distance between D44 and R55, while the distance between F47 and R55 remained stable, suggesting the formation of a cation-pi interaction between those residues. Experimental double mutant cycles confirmed that the F47-R55 pair has a larger energetic coupling than the D44-R55 interaction. Overall, these combined experimental and computational studies show that the X-ray crystal structure is the better reference structure for HP36 in solution at neutral pH. Our analysis also shows how detailed molecular dynamics simulations combined with experimental validation can help bridge the gap between NMR and crystallographic methods.

Parallel tempering or replica-exchange molecular dynamics (REMD) significantly increases the efficiency of conformational sampling for complex molecular systems. However, obtaining converged data with REMD remains challenging, especially for large systems with complex topologies. We propose a new variant to REMD where the replicas are also permitted to exchange with an ensemble of structures that have been generated in advance using high-temperature MD simulations, similar in spirit to J-walking methods. We tested this approach on two non-trivial model systems, a β-hairpin and a 3-stranded β-sheet and compared the results to those obtained from very long (>100 ns) standard REMD simulations. The resulting ensembles were indistinguishable, including relative populations of different conformations on the unfolded state. The use of the reservoir is shown to significantly reduce the time required for convergence.

The effects of the use of three generalized Born (GB) implicit solvent models on the thermodynamics of a simple polyalanine peptide are studied via comparing several hundred nanoseconds of well-converged replica exchange molecular dynamics (REMD) simulations using explicit TIP3P solvent to REMD simulations with the GB solvent models. It is found that when compared to REMD simulations using TIP3P the GB REMD simulations contain significant differences in secondary structure populations, most notably an overabundance of alpha-helical secondary structure. This discrepancy is explored via comparison of the differences in the electrostatic component of the free energy of solvation (DeltaDeltaG(pol)) between TIP3P (via thermodynamic Integration calculations), the GB models, and an implicit solvent model based on the Poisson equation (PE). The electrostatic components of the solvation free energies are calculated using each solvent model for four representative conformations of Ala10. Since the PE model is found to have the best performance with respect to reproducing TIP3P DeltaDeltaG(pol) values, effective Born radii from the GB models are compared to effective Born radii calculated with PE (so-called perfect radii), and significant and numerous deviations in GB radii from perfect radii are found in all GB models. The effect of these deviations on the solvation free energy is discussed, and it is shown that even when perfect radii are used the agreement of GB with TIP3P DeltaDeltaG(pol) values does not improve. This suggests a limit to the optimization of the effective Born radius calculation and that future efforts to improve the accuracy of GB models must extend beyond such optimizations.

HIV-1 protease remains an important anti-AIDS drug target. Although it has been known that ligand binding induces large conformational changes in the protease, the dynamic aspects of binding have been largely ignored. Several computational models describing protease dynamics have been reported recently. These have reproduced experimental observations, and have also explained how ligands gain access to the binding site through dynamic behavior of the protease. Specifically, the transitions between three different conformations of the protein have been modeled in atomic detail. Two of these forms were determined by crystallography, and the third was implied by NMR experiments. Based on these computational models, it has been suggested that binding of inhibitors in allosteric sites might affect protease flexibility and disrupt its function.

The ff94 force field that is commonly associated with the Amber simulation package is one of the most widely used parameter sets for biomolecular simulation. After a decade of extensive use and testing, limitations in this force field, such as over-stabilization of alpha-helices, were reported by us and other researchers. This led to a number of attempts to improve these parameters, resulting in a variety of "Amber" force fields and significant difficulty in determining which should be used for a particular application. We show that several of these continue to suffer from inadequate balance between different secondary structure elements. In addition, the approach used in most of these studies neglected to account for the existence in Amber of two sets of backbone phi/psi dihedral terms. This led to parameter sets that provide unreasonable conformational preferences for glycine. We report here an effort to improve the phi/psi dihedral terms in the ff99 energy function. Dihedral term parameters are based on fitting the energies of multiple conformations of glycine and alanine tetrapeptides from high level ab initio quantum mechanical calculations. The new parameters for backbone dihedrals replace those in the existing ff99 force field. This parameter set, which we denote ff99SB, achieves a better balance of secondary structure elements as judged by improved distribution of backbone dihedrals for glycine and alanine with respect to PDB survey data. It also accomplishes improved agreement with published experimental data for conformational preferences of short alanine peptides and better accord with experimental NMR relaxation data of test protein systems., ((c) 2006 Wiley-Liss, Inc.)

The introduction of HIV-1 protease (HIV-PR) inhibitors has led to a dramatic increase in patient survival; however, these gains are threatened by the emergence of multi-drug-resistant strains. Design of inhibitors that overcome resistance would be greatly facilitated by deeper insight into the mechanistic events associated with binding of substrates and inhibitors, as well as an understanding of the effects of resistance mutations on the structure and dynamic behavior of HIV-PR. We previously reported a series of simulations that provide a model for HIV-PR dynamics, with spontaneous conversions between the bound and unbound crystal forms upon addition or removal of an inhibitor. Importantly, the unbound protease transiently sampled a third fully open state that permits entry to the active site, unlike both crystallographic forms. Recently, a crystal structure of unbound HIV-PR was reported for the MDR 769 isolate (PDB: 1TW7); unlike all previous experimental structures, the binding pocket is open. It is suggested that drug resistance in this strain arises at least in part from the inability of inhibitors to induce closing. We carried out simulations of the MDR 769 HIV-PR mutant and observed that the reported structure is unstable in solution and rapidly adopts the semi-open conformation observed for the unbound wild-type protease in solution. Further analysis suggests that the wide-open structure observed for MDR 769 arises not from sequence variation, but instead is an artifact from crystal packing. Thus, despite being the first experimental structure to reveal flap opening sufficient for substrate access to the active site, this structure may not be directly relevant to studies of inhibitor entry or to the cause of HIV-PR drug resistance.

8-Oxoguanine (8OG) is the most prevalent form of oxidative DNA damage. In bacteria, 8OG is excised by formamidopyrimidine glycosylase (Fpg) as the initial step in base excision repair. To efficiently excise this lesion, Fpg must discriminate between 8OG and an excess of guanine in duplex DNA. In this study, we explore the structural basis underlying this high degree of selectivity. Two structures have been reported in which Fpg is bound to DNA, differing with respect to the position of the lesion in the active site, one structure showing 8OG bound in the syn conformation and the other in the anti conformation. Remarkably, the results of our all-atom simulations are consistent with both structures. The syn conformation observed in the crystallographic structure of Fpg obtained from Bacillus stearothermophilus is stabilized through interaction with E77, a nonconserved residue. Replacement of E77 with Ser, creating the Fpg sequence found in Escherichia coli and other bacteria, results in preferred binding of 8OG in the anti conformation. Our calculations provide novel insights into the roles of active site residues in binding and recognition of 8OG by Fpg.

The 36 residue villin headpiece helical subdomain (HP36) is one of the fastest cooperatively folding proteins, folding on the microsecond timescale. HP36's simple three helix topology, fast folding and small size have made it an attractive model system for computational and experimental studies of protein folding. Recent experimental studies have explored the denatured state of HP36 using fragment analysis coupled with relatively low-resolution spectroscopic techniques. These studies have shown that there is apparently only a small tendency to form locally stabilized secondary structure. Here, we complement the experimental studies by using replica exchange molecular dynamics with explicit solvent to investigate the structural features of these peptide models of unfolded HP36. To ensure convergence, two sets of simulations for each fragment were performed with different initial structures, and simulations were continued until these generated very similar final ensembles. These simulations reveal low populations of native-like structure and early folding events that cannot be resolved by experiment. For each fragment, calculated J-coupling constants and helical propensities are in good agreement with experimental trends. HP-1, corresponding to residues 41 to 53 and including the first alpha-helix, contains the highest helical population. HP-3, corresponding to residues 62 through 75 and including the third alpha-helix, contains a small population of helical turn residing at the N terminus while HP-2, corresponding to residues 52 through 61 and including the second alpha-helix, formed little to no structure in isolation. Overall, HP-1 was the only fragment to adopt a native-like conformation, but the low population suggests that formation of significant structure only occurs after formation of specific tertiary interactions.

Isomerization of the 11-cis retinal chromophore in the visual pigment rhodopsin is coupled to motion of transmembrane helix H6 and receptor activation. We present solid-state magic angle spinning NMR measurements of rhodopsin and the metarhodopsin II intermediate that support the proposal that interaction of Trp265(6.48) with the retinal chromophore is responsible for stabilizing an inactive conformation in the dark, and that motion of the beta-ionone ring allows Trp265(6.48) and transmembrane helix H6 to adopt active conformations in the light. Two-dimensional dipolar-assisted rotational resonance NMR measurements are made between the C19 and C20-methyl groups of the retinal and uniformly 13C-labeled Trp265(6.48). The retinal C20-Trp265(6.48) contact present in the dark-state of rhodopsin is lost in metarhodopsin II, and a new contact is formed with the C19 methyl group. We have previously shown that the retinal translates 4-5 A toward H5 in metarhodopsin II. This motion, in conjunction with the Trp-C19 contact, implies that the Trp265(6.48) side-chain moves significantly upon rhodopsin activation. NMR measurements also show that a packing interaction in rhodopsin between Trp265(6.48) and Gly121(3.36) is lost in metarhodopsin II, consistent with H6 motion away from H3. However, a close contact between Gly120(3.35) on H3 and Met86(2.53) on H2 is observed in both rhodopsin and metarhodopsin II, suggesting that H3 does not change orientation significantly upon receptor activation.

We report unrestrained, all-atom molecular dynamics simulations of HIV-1 protease (HIV-PR) with a continuum solvent model that reproducibly sample closing of the active site flaps following manual placement of a cyclic urea inhibitor into the substrate binding site of the open protease. The open form was obtained from the unbound, semi-open HIV-PR crystal structure, which we recently reported (Hornak, V.; et al. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 915-920.) to have spontaneously opened during unrestrained dynamics. In those simulations, the transiently open flaps always returned to the semi-open form that is observed in all crystal structures of the free protease. Here, we show that manual docking of the inhibitor reproducibly induces spontaneous conversion to the closed form as seen in all inhibitor-bound HIV-PR crystal structures. These simulations reproduced not only the greater degree of flap closure, but also the striking difference in flap "handedness" between bound and free enzyme. In most of the simulations, the final structures were highly accurate. Root-mean-square deviations (RMSD) from the crystal structure of the complex were approximately 1.5 A (averaged over the last 100 ps) for the inhibitor and each flap despite initial RMSD of 2-5 A for the inhibitors and 6-11 A for the flaps. Key hydrogen bonds were formed between the flap tips and between flaps and inhibitor that match those seen in the crystal structure. The results demonstrate that all-atom simulations have the ability to significantly improve poorly docked ligand conformations and reproduce large-scale receptor conformational changes that occur upon binding.

The flexibility of the conserved 57DTAGQ61 motif is essential for Ras proper cycling in response to growth factors. Here, we increase the flexibility of the 57DTAGQ61 motif by mutating Gln61 to Gly. The crystal structure of the RasQ61G mutant reveals a new conformation of switch 2 that bears remarkable structural homology to an intermediate for GTP hydrolysis revealed by targeted molecular dynamics simulations. The mutation increased retention of GTP and inhibited Ras binding to the catalytic site, but not to the distal site of Sos. Most importantly, the thermodynamics of RafRBD binding to Ras are altered even though the structure of switch 1 is not affected by the mutation. Our results suggest that interplay and transmission of structural information between the switch regions are important factors for Ras function. They propose that initiation of GTP hydrolysis sets off the separation of the Ras/effector complex even before the GDP conformation is reached.

The use of parallel tempering or replica exchange molecular dynamics (REMD) simulations has facilitated the exploration of free energy landscapes for complex molecular systems, but application to large systems is hampered by the scaling of the number of required replicas with increasing system size. Use of continuum solvent models reduces system size and replica requirements, but these have been shown to provide poor results in many cases, including overstabilization of ion pairs and secondary structure bias. Hybrid explicit/continuum solvent models can overcome some of these problems through an explicit representation of water molecules in the first solvation shells, but these methods typically require restraints on the solvent molecules and show artifacts in water properties due to the solvation interface. We propose an REMD variant in which the simulations are performed with a fully explicit solvent, but the calculation of exchange probability is carried out using a hybrid model, with the solvation shells calculated on the fly during the fully solvated simulation. The resulting reduction in the perceived system size in the REMD exchange calculation provides a dramatic decrease in the computational cost of REMD, while maintaining a very good agreement with results obtained from the standard explicit solvent REMD. We applied several standard and hybrid REMD methods with different solvent models to alanine polymers of 1, 3, and 10 residues, obtaining ensembles that were essentially independent of the initial conformation, even with explicit solvation. Use of only a continuum model without a shell of explicit water provided poor results for Ala3 and Ala10, with a significant bias in favor of the α-helix. Likewise, using only the solvation shells and no continuum model resulted in ensembles that differed significantly from the standard explicit solvent data. Ensembles obtained from hybrid REMD are in very close agreement with explicit solvent data, predominantly populating polyproline II conformations. Inclusion of a second shell of explicit solvent was found to be unnecessary for these peptides.

We report unrestrained, all-atom molecular dynamics simulations of HIV-1 protease that sample large conformational changes of the active site flaps. In particular, the unliganded protease undergoes multiple conversions between the "closed" and "semiopen" forms observed in crystal structures of inhibitor-bound and unliganded protease, respectively, including reversal of flap "handedness." Simulations in the presence of a cyclic urea inhibitor yield stable closed flaps. Furthermore, we observe several events in which the flaps of the unliganded protease open to a much greater degree than observed in crystal structures and subsequently return to the semiopen state. Our data strongly support the hypothesis that the unliganded protease predominantly populates the semiopen conformation, with closed and fully open structures being a minor component of the overall ensemble. The results also provide a model for the flap opening and closing that is considered to be essential to enzyme function.

Potentials of mean force (PMFs) of salt bridge formation between oppositely charged amino acid side chains were calculated both in explicit solvent and in a Generalized Born (GB) continuum solvent model to quantify the potential overstabilization of side chain ion pairs in GB relative to explicit solvation. These show that salt bridges are too stable by as much as 3-4 kcal/mol in the GB solvent models that we tested, consistent with previously reported observations of significantly different structural ensembles in GB models and explicit solvent for proteins containing ionizable groups. We thus investigated a simple empirical correction, wherein the intrinsic GB radii of hydrogen atoms bound to charged nitrogen atoms are reduced, effectively increasing the desolvation penalty of the positively charged groups. The thermodynamics of salt bridge formation were considerably improved, as exemplified by the close match of the corrected GB PMF to the reference explicit solvent PMF, and more significantly by our ability to closely reproduce the experimental temperature melting profile of the TC5b Trp-cage miniprotein, which is otherwise highly distorted by prevalent non-native salt bridges when using standard GB parameters.

The process by which DNA repair enzymes recognize and selectively excise damaged bases in duplex DNA is fundamental to our mechanistic understanding of these critical biological reactions. 8-Oxoguanine (8-oxoG) is the most common form of oxidative DNA damage; unrepaired, this lesion generates a G:C-->T:A mutation. Central to the recognition and repair of DNA damage is base extrusion, a process in which the damaged base lesion or, in some cases, its partner disengages from the helix and is bound to the enzyme's active site where base excision takes place. The conformation adopted by 8-oxoG in duplex DNA is affected by the base positioned opposite this lesion; conformational changes may also take place when the damaged base binds to its cognate repair enzyme. We performed unrestrained molecular dynamics simulations for several 13-mer DNA duplexes. Oligomers containing G:C and 8oxoG:C pairs adopted Watson-Crick geometries in stable B-form duplexes; 8oxoG showed increased local and global flexibility and a reduced barrier to base extrusion. Duplexes containing the G:A mismatch showed much larger structural fluctuations and failed to adopt a well-defined structure. For the 8oxoG:A mismatch that is recognized by the DNA glycosylase MutY, the damaged nucleoside underwent spontaneous and reproducible anti-->syn transitions. The syn conformation is thermodynamically preferred. Steric hindrance and unfavorable electrostatics associated with the 8oxoG O8 atom in the anti conformation were the major driving forces for this transition. Transition events follow two qualitatively different pathways. The overall anti-->syn transition rate and relative probability of the two transition paths were dependent on local sequence context. These simulations indicate that both the dynamic and equilibrium behavior of the duplex change as a result of oxidation; these differences may provide valuable new insight into the selective action of enzymes on damaged DNA.

The thermodynamic behavior of a previously designed three-stranded beta-sheet was studied via several microseconds of standard and replica exchange molecular dynamics simulations. The system is shown to populate at least four thermodynamic minima, including two partially folded states in which only a single hairpin is formed. Simulated melting curves show different profiles for the C and N-terminal hairpins, consistent with differences in secondary structure content in published NMR and CD/FTIR measurements, which probed different regions of the chain. Individual beta-hairpins that comprise the three-stranded beta-sheet are observed to form cooperatively. Partial folding cooperativity between the component hairpins is observed, and good agreement between calculated and experimental values quantifying this cooperativity is obtained when similar analysis techniques are used. However, the structural detail in the ensemble of conformations sampled in the simulations permits a more direct analysis of this cooperativity than has been performed on the basis of experimental data. The results indicate the actual folding cooperativity perpendicular to strand direction is significantly larger than the lower bound obtained previously.