A Three‐Site Mechanism for Agonist/Antagonist Selective Binding to Vasopressin Receptors

Molecular‐dynamics simulations with metadynamics enhanced sampling reveal three distinct binding sites for arginine vasopressin (AVP) within its V₂‐receptor (V₂R). Two of these, the vestibule and intermediate sites, block (antagonize) the receptor, and the third is the orthosteric activation (agonist) site. The contacts found for the orthosteric site satisfy all the requirements deduced from mutagenesis experiments. Metadynamics simulations for V₂R and its V_1aR‐analog give an excellent correlation with experimental binding free energies by assuming that the most stable binding site in the simulations corresponds to the experimental binding free energy in each case. The resulting three‐site mechanism separates agonists from antagonists and explains subtype selectivity.     

DIRECT‐ID: An automated method to identify and quantify conformational variations—application to β2‐adrenergic GPCR

The conformational dynamics of a macromolecule can be modulated by a number of factors, including changes in environment, ligand binding, and interactions with other macromolecules, among others. We present a method that quantifies the differences in macromolecular conformational dynamics and automatically extracts the structural features responsible for these changes. Given a set of molecular dynamics (MD) simulations of a macromolecule, the norms of the differences in covariance matrices are calculated for each pair of trajectories. A matrix of these norms thus quantifies the differences in conformational dynamics across the set of simulations. For each pair of trajectories, covariance difference matrices are parsed to extract structural elements that undergo changes in conformational properties. As a demonstration of its applicability to biomacromolecular systems, the method, referred to as DIRECT‐ID, was used to identify relevant ligand‐modulated structural variations in the β₂‐adrenergic (β₂AR) G‐protein coupled receptor. Micro‐second MD simulations of the β₂AR in an explicit lipid bilayer were run in the apo state and complexed with the ligands: BI‐167107 (agonist), epinephrine (agonist), salbutamol (long‐acting partial agonist), or carazolol (inverse agonist). Each ligand modulated the conformational dynamics of β₂AR differently and DIRECT‐ID analysis of the inverse‐agonist vs. agonist‐modulated β₂AR identified residues known through previous studies to selectively propagate deactivation/activation information, along with some previously unidentified ligand‐specific microswitches across the GPCR. This study demonstrates the utility of DIRECT‐ID to rapidly extract functionally relevant conformational dynamics information from extended MD simulations of large and complex macromolecular systems. © 2015 Wiley Periodicals, Inc.     

Binding structures of tri‐N‐acetyl‐β‐glucosamine in hen egg white lysozyme using molecular dynamics with a polarizable force field

Lysozyme is a well‐studied enzyme that hydrolyzes the β‐(1,4)‐glycosidic linkage of N‐acetyl‐β‐glucosamine (NAG)_n oligomers. The active site of hen egg‐white lysozyme (HEWL) is believed to consist of six subsites, A‐F that can accommodate six sugar residues. We present studies exploring the use of polarizable force fields in conjunction with all‐atom molecular dynamics (MD) simulations to analyze binding structures of complexes of lysozyme and NAG trisaccharide, (NAG)₃. MD trajectories are applied to analyze structures and conformation of the complex as well as protein–ligand interactions, including the hydrogen‐bonding network in the binding pocket. Two binding modes (ABC and BCD) of (NAG)₃ are investigated independently based on a fixed‐charge model and a polarizable model. We also apply molecular mechanics with generalized born and surface area (MM‐GBSA) methods based on MD using both nonpolarizable and polarizable force fields to compute binding free energies. We also study the correlation between root‐mean‐squared deviation and binding free energies of the wildtype and W62Y mutant; we find that for this prototypical system, approaches using the MD trajectories coupled with implicit solvent models are equivalent for polarizable and fixed‐charge models. © 2012 Wiley Periodicals, Inc.     

Comparison of radii sets,entropy, QM methods,and sampling on MM‐PBSA,MM‐GBSA,and QM/MM‐GBSA ligand binding energies of F. tularensis enoyl‐ACP reductase (FabI)

To validate a method for predicting the binding affinities of FabI inhibitors, three implicit solvent methods, MM‐PBSA, MM‐GBSA, and QM/MM‐GBSA were carefully compared using 16 benzimidazole inhibitors in complex with Francisella tularensis FabI. The data suggests that the prediction results are sensitive to radii sets, GB methods, QM Hamiltonians, sampling protocols, and simulation length, if only one simulation trajectory is used for each ligand. In this case, QM/MM‐GBSA using 6 ns MD simulation trajectories together with GB^neck2, PM3, and the mbondi2 radii set, generate the closest agreement with experimental values (r² = 0.88). However, if the three implicit solvent methods are averaged from six 1 ns MD simulations for each ligand (called “multiple independent sampling”), the prediction results are relatively insensitive to all the tested parameters. Moreover, MM/GBSA together with GB^HCT and mbondi, using 600 frames extracted evenly from six 0.25 ns MD simulations, can also provide accurate prediction to experimental values (r² = 0.84). Therefore, the multiple independent sampling method can be more efficient than a single, long simulation method. Since future scaffold expansions may significantly change the benzimidazole's physiochemical properties (charges, etc.) and possibly binding modes, which may affect the sensitivities of various parameters, the relatively insensitive “multiple independent sampling method” may avoid the need of an entirely new validation study. Moreover, due to large fluctuating entropy values, (QM/)MM‐P(G)BSA were limited to inhibitors’ relative affinity prediction, but not the absolute affinity. The developed protocol will support an ongoing benzimidazole lead optimization program. © 2015 Wiley Periodicals, Inc.     

Identification of Two Distinct Sites for Antagonist and Biased Agonist Binding to the Human Chemokine Receptor CXCR3

The chemokine receptor CXCR3 is a G protein‐coupled receptor that conveys extracellular signals into cells by changing its conformation upon ligand binding. We previously hypothesized that small‐molecule allosteric CXCR3‐agonists do not bind to the same allosteric binding pocket as 8‐azaquinazolinone‐based negative allosteric modulators. We have now performed molecular‐dynamics (MD) simulations with metadynamics enhanced sampling on the CXCR3 system to refine structures and binding modes and to predict the CXCR3‐binding affinities of the biased allosteric agonist FAUC1036 and the negative allosteric modulator RAMX3. We have identified two distinct binding sites; a “shallow” and a second “deeper” pocket to which the biased allosteric agonist FAUC1036 and negative allosteric modulator RAMX3 bind, respectively.     

Structural dissimilarity sampling with dynamically self‐guiding selection

Structural dissimilarity sampling (SDS) has been proposed as an enhanced conformational sampling method for reproducing the structural transitions of a given protein. SDS consists of cycles of two steps: (1) Selections of initial structures with structural dissimilarities by referring to a measure. (2) Conformational resampling by restarting short‐time molecular dynamics (MD) simulations from the initial structures. In the present study, an efficient measure is proposed as a dynamically self‐guiding selection to accelerate the structural transitions from a reactant state to a product state as an extension to the original SDS. In the extended SDS, the inner product (IP ) between the reactant and the snapshots generated by short‐time MD simulations are evaluated and ranked according to the IP s at every cycle. Then, the snapshots with low IP s are selected as initial structures for the short‐time MD simulations. This scheme enables one to choose dissimilar and distant initial structures from the reactant, and thus the initial structures dynamically head towards the product, promoting structural transitions from the reactant. To confirm the conformational sampling efficiency, the extended SDS was applied to maltodextrin binding protein (MBP), and we successfully reproduced the structural transition from the open to closed states with submicrosecond‐order simulation times. However, a conventional long‐time MD simulation failed to reproduce the same structural transition. We also compared the performance with that obtained by the ordinary SDS and other sampling techniques that have been developed by us to characterize the possible utility of the extended SDS for actual applications. © 2017 Wiley Periodicals, Inc.     

Atomistic detailed free‐energy landscape of intrinsically disordered protein studied by multi‐scale divide‐and‐conquer molecular dynamics simulation

Calcineurin (CaN) is a eukaryotic serine/threonine protein phosphatase activated by both Ca²⁺ and calmodulin (CaM), including intrinsically disordered region (IDR). The region undergoes folding into an α‐helix form in the presence Ca²⁺‐loaded CaM. To sample the ordered structure of the IDR by conventional all atom model (AAM) molecular dynamics (MD) simulation, the IDR and Ca²⁺‐loaded CaM must be simultaneously treated. However, it is time‐consuming task because the coupled folding and binding should include repeated binding and dissociation. Then, in this study, we propose novel multi‐scale divide‐and‐conquer MD (MSDC‐MD), which combines AAM‐MD and coarse‐grained model MD (CGM‐MD). To speed up the conformation sampling, MSDC‐MD simulation first treats the IDR by CGM to sample conformations from wide conformation space; then, multiple AAM‐MD in a limited area is initiated using the resultant CGM conformation, which is reconstructed by homology modeling method. To investigate performance, we sampled the ordered conformation of the IDR using MSDC‐MD; the root‐mean‐square distance (RMSD) with respect to the experimental structure was 2.23 Å.     

Differences between G‐Protein‐Stabilized Agonist–GPCR Complexes and their Nanobody‐Stabilized Equivalents

Protein nanobodies have been used successfully as surrogates for unstable G‐proteins in order to crystallize G‐protein‐coupled receptors (GPCRs) in their active states. We used molecular dynamics (MD) simulations, including metadynamics enhanced sampling, to investigate the similarities and differences between GPCR–agonist ternary complexes with the α‐subunits of the appropriate G‐proteins and those with the protein nanobodies (intracellular binding partners, IBPs) used for crystallization. In two of the three receptors considered, the agonist‐binding mode differs significantly between the two alternative ternary complexes. The ternary‐complex model of GPCR activation entails enhancement of ligand binding by bound IBPs: Our results show that IBP‐specific changes can alter the agonist binding modes and thus also the criteria for designing GPCR agonists.     

NMR Structure of a Triangulenium‐Based Long‐Lived Fluorescence Probe Bound to a G‐Quadruplex

An NMR structural study of the interaction between a small‐molecule optical probe (DAOTA‐M2) and a G‐quadruplex from the promoter region of the c‐myc oncogene revealed that they interact at 1:2 binding stoichiometry. NMR‐restrained structural calculations show that binding of DAOTA‐M2 occurs mainly through π–π stacking between the polyaromatic core of the ligand and guanine residues of the outer G‐quartets. Interestingly, the binding affinities of DAOTA‐M2 differ by a factor of two for the outer G‐quartets of the unimolecular parallel G‐quadruplex under study. Unrestrained MD calculations indicate that DAOTA‐M2 displays significant dynamic behavior when stacked on a G‐quartet plane. These studies provide molecular guidelines for the design of triangulenium derivatives that can be used as optical probes for G‐quadruplexes.     

Insight into the structural requirements of benzimidazole derivatives as interleukin‐2 inducible t‐cell kinase inhibitors by computational explorations

In the present work, a set of ligand‐ and receptor‐based 3D‐QSAR models were developed to explore the structure–activity relationship of 109 benzimidazole‐based interleukin‐2‐inducible T‐cell kinase (ITK) inhibitors. In order to reveal the requisite 3D structural features impacting the biological activities, a variety of in silico modeling approaches including the comparative molecular field analysis (CoMFA), comparative molecular similarity indices analysis (CoMSIA), docking, and molecular dynamics were applied. The results showed that the ligand‐based CoMFA model (Q² = 0.552, R²_ncv = 0.908, R²_pred = 0.787, SEE = 0.252, SEP = 0.558) and CoMSIA model (Q² = 0.579, R²_ncv = 0.914, R²_pred = 0.893, SEE = 0.240, SEP = 0.538) were superior to other models with greater predictive power. In addition, a combined analysis between the 3D contour maps and docking results showed that: (1) Compounds with bulky or hydrophobic substituents near ring D and electropositive or hydrogen acceptor groups around rings C and D could increase the activity. (2) The key amino acids impacting the receptor–ligand interactions in the binding pocket are Met438, Asp500, Lys391, and Glu439. The results obtained from this work may provide helpful guidelines in design of novel benzimidazole analogs as inhibitors of ITK. © 2013 Wiley Periodicals, Inc.     

Calculation of free energy landscapes: A histogram reweighted metadynamics approach

We present an efficient method for the calculation of free energy landscapes. Our approach involves a history‐dependent bias potential, which is evaluated on a grid. The corresponding free energy landscape is constructed via a histogram reweighting procedure a posteriori. Because of the presence of the bias potential, it can be also used to accelerate rare events. In addition, the calculated free energy landscape is not restricted to the actual choice of collective variables and can in principle be extended to auxiliary variables of interest without further numerical effort. The applicability is shown for several examples. We present numerical results for the alanine dipeptide and the Met‐Enkephalin in explicit solution to illustrate our approach. Furthermore, we derive an empirical formula that allows the prediction of the computational cost for the ordinary metadynamics variant in comparison with our approach, which is validated by a dimensionless representation. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Drug release by pH‐responsive molecular tweezers: Atomistic details from molecular modeling

pH‐responsive molecular tweezers have been proposed as an approach for targeting drug‐delivery to tumors, which tend to have a lower pH than normal cells. We performed a computational study of a pH‐responsive molecular tweezer using ab initio quantum chemistry in the gas‐phase and molecular dynamics (MD) simulations in solution. The binding free energy in solution was calculated using steered MD. We observe, in atomistic detail, the pH‐induced conformational switch of the tweezer and the resulting release of the drug molecule. Even when the tweezer opens, the drug molecule remains near a hydrophobic arm of the molecular tweezer. Drug release cannot occur, it seems, unless the tweezer is in a hydrophobic environment with low pH. © 2014 Wiley Periodicals, Inc.     

Multidimensional De Novo Design Reveals 5‐HT2B Receptor‐Selective Ligands

We report a multi‐objective de novo design study driven by synthetic tractability and aimed at the prioritization of computer‐generated 5‐HT_2B receptor ligands with accurately predicted target‐binding affinities. Relying on quantitative bioactivity models we designed and synthesized structurally novel, selective, nanomolar, and ligand‐efficient 5‐HT_2B modulators with sustained cell‐based effects. Our results suggest that seamless amalgamation of computational activity prediction and molecular design with microfluidics‐assisted synthesis enables the swift generation of small molecules with the desired polypharmacology.     

Using the local elevation method to construct optimized umbrella sampling potentials: Calculation of the relative free energies and interconversion barriers of glucopyranose ring conformers in water

A method is proposed to combine the local elevation (LE) conformational searching and the umbrella sampling (US) conformational sampling approaches into a single local elevation umbrella sampling (LEUS) scheme for (explicit‐solvent) molecular dynamics (MD) simulations. In this approach, an initial (relatively short) LE build‐up (searching) phase is used to construct an optimized biasing potential within a subspace of conformationally relevant degrees of freedom, that is then used in a (comparatively longer) US sampling phase. This scheme dramatically enhances (in comparison with plain MD) the sampling power of MD simulations, taking advantage of the fact that the preoptimized biasing potential represents a reasonable approximation to the negative of the free energy surface in the considered conformational subspace. The method is applied to the calculation of the relative free energies of β‐D ‐glucopyranose ring conformers in water (within the GROMOS 45A4 force field). Different schemes to assign sampled conformational regions to distinct states are also compared. This approach, which bears some analogies with adaptive umbrella sampling and metadynamics (but within a very distinct implementation), is shown to be: (i) efficient (nearly all the computational effort is invested in the actual sampling phase rather than in searching and equilibration); (ii) robust (the method is only weakly sensitive to the details of the build‐up protocol, even for relatively short build‐up times); (iii) versatile (a LEUS biasing potential database could easily be preoptimized for small molecules and assembled on a fragment basis for larger ones). © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

X‐Ray Crystallography and Free Energy Calculations Reveal the Binding Mechanism of A2A Adenosine Receptor Antagonists

We present a robust protocol based on iterations of free energy perturbation (FEP) calculations, chemical synthesis, biophysical mapping and X‐ray crystallography to reveal the binding mode of an antagonist series to the A_2A adenosine receptor (AR). Eight A_2AAR binding site mutations from biophysical mapping experiments were initially analyzed with sidechain FEP simulations, performed on alternate binding modes. The results distinctively supported one binding mode, which was subsequently used to design new chromone derivatives. Their affinities for the A_2AAR were experimentally determined and investigated through a cycle of ligand‐FEP calculations, validating the binding orientation of the different chemical substituents proposed. Subsequent X‐ray crystallography of the A_2AAR with a low and a high affinity chromone derivative confirmed the predicted binding orientation. The new molecules and structures here reported were driven by free energy calculations, and provide new insights on antagonist binding to the A_2AAR, an emerging target in immuno‐oncology.     

Theoretical Insights into Chirality‐Controlled SWCNT Growth from a Cycloparaphenylene Template

A self‐assembly mechanism for low‐temperature SWCNT growth from a [6]cycloparaphenylene ([6]CPP) precursor via ethynyl (C₂H) radical addition is presented, based on non‐equilibrium quantum chemical molecular dynamics (QM/MD) simulations and density functional theory (DFT) calculations. This mechanism, which maintains the (6,6) armchair chirality of a SWCNT fragment throughout the growth process, is energetically more favorable than a previously proposed Diels–Alder‐based growth mechanisms [E. H. Fort, et al., J. Mater. Chem. 2011 , 21, 1373]. QM/MD simulations and DFT calculations show that C₂H radicals play dual roles during SWCNT growth, by abstracting hydrogen from the SWCNT fragment and providing the carbon source necessary for growth itself. Simulations demonstrate that chirality‐controlled SWCNT growth from macrocyclic hydrocarbon seed molecules with pre‐selected edge structure can be accomplished when the reaction conditions are carefully selected for hydrogen abstraction by radical species during the growth process.     

Conformational Dynamics of apo‐GlnBP Revealed by Experimental and Computational Analysis

The glutamine binding protein (GlnBP) binds l ‐glutamine and cooperates with its cognate transporters during glutamine uptake. Crystal structure analysis has revealed an open and a closed conformation for apo‐ and holo‐GlnBP, respectively. However, the detailed conformational dynamics have remained unclear. Herein, we combined NMR spectroscopy, MD simulations, and single‐molecule FRET techniques to decipher the conformational dynamics of apo‐GlnBP. The NMR residual dipolar couplings of apo‐GlnBP were in good agreement with a MD‐derived structure ensemble consisting of four metastable states. The open and closed conformations are the two major states. This four‐state model was further validated by smFRET experiments and suggests the conformational selection mechanism in ligand recognition of GlnBP.