Molecular modeling of organic and biomolecular systems using BOSS and MCPRO

Application of a BOSS-Gaussian interface for QM/MM simulations of Henry and methyl transfer reactions

Hybrid quantum mechanics and molecular mechanics (QM/MM) computer simulations have become an indispensable tool for studying chemical and biological phenomena for systems too large to treat with QM alone. For several decades, semiempirical QM methods have been used in QM/MM simulations. However, with increased computational resources, the introduction of ab initio and density function methods into on-the-fly QM/MM simulations is being increasingly preferred. This adaptation can be accomplished with a program interface that tethers independent QM and MM software packages. This report introduces such an interface for the BOSS and Gaussian programs, featuring modification of BOSS to request QM energies and partial atomic charges from Gaussian. A customizable C-shell linker script facilitates the interprogram communication. The BOSS-Gaussian interface also provides convenient access to Charge Model 5 (CM5) partial atomic charges for multiple purposes including QM/MM studies of reactions. In this report, the BOSS-Gaussian interface is applied to a nitroaldol (Henry) reaction and two methyl transfer reactions in aqueous solution. Improved agreement with experiment is found by determining free-energy surfaces with MP2/CM5 QM/MM simulations than previously reported investigations using semiempirical methods.

Importance of topology for glycocluster binding to Pseudomonas aeruginosa and Burkholderia ambifaria bacterial lectins

Pseudomonas aeruginosa (PA) and Burkholderia ambifaria (BA) are two opportunistic Gram negative bacteria and major infectious agents involved in lung infection of cystic fibrosis patients. Both bacteria can develop resistance to conventional antibiotherapies. An alternative strategy consists of targeting virulence factors in particular lectins with high affinity ligands such as multivalent glycoclusters. LecA (PA-IL) and LecB (PA-IIL) are two tetravalent lectins from PA that recognise galactose and fucose respectively. BambL lectin from BA is trimeric with 2 binding sites per monomer and is also specific for fucose. These three lectins are potential therapeutic targets in an anti-adhesive anti-bacterial approach. Herein, we report the synthesis of 18 oligonucleotide pentofuranose-centered or mannitol-centered glycoclusters leading to tri-, penta- or decavalent clusters with different topologies. The linker arm length between the core and the carbohydrate epitope was also varied leading to 9 galactoclusters targeting LecA and 9 fucoclusters targeting both LecB and BambL. Their dissociation constants (Kd) were determined using a DNA-based carbohydrate microarray technology. The trivalent xylo-centered galactocluster and the ribo-centered fucocluster exhibited the best affinity for LecA and LecB respectively while the mannitol-centered decafucocluster displayed the best affinity to BambL. These data demonstrated that the topology and nature of linkers were the predominant factors for achieving high affinity rather than valency.

Evaluation of CM5 Charges for Nonaqueous Condensed-Phase Modeling

Partial atomic charges for neutral molecules from quantum mechanical calculations are typically scaled for use in molecular modeling of liquid-phase systems. Optimal scale factors of 1.14 for CM1A and 1.27 for CM5 charges were previously determined for minimizing errors in free energies of hydration. The adequacy of the 1.14*CM1A and 1.27*CM5 models are evaluated here in pure liquid simulations in combination with the OPLS-AA force field. For 22 organic liquids, the 1.14*CM1A and 1.27*CM5 models yield mean unsigned errors (MUEs) of ca. 1.40 kcal/mol for heats of vaporization. Not surprisingly, this reflects overpolarization with the scale factors derived for aqueous media. Prediction of pure liquid properties using CM5 charges is optimized using a scale factor of 1.14, which reduces the MUE for heats of vaporization to 0.89 kcal/mol. However, due to the impracticality of using different scale factors in different explicit-solvent condensed-phase simulations, a universal scale factor of 1.20 emerged for CM5 charges. This provides a balance between errors in computed pure liquid properties and free energies of hydration. Computation of free energies of hydration by the GB/SA method further found that 1.20 is equally suited for use in explicit or implicit treatments of aqueous solvation. With 1.20*CM5 charges, a variety of condensed-phase simulations can be pursued while maintaining average errors of 1.0 kcal/mol in key thermodynamic properties.

Improved Peptide and Protein Torsional Energetics with the OPLSAA Force Field

The development and validation of new peptide dihedral parameters are reported for the OPLS-AA force field. High accuracy quantum chemical methods were used to scan φ, ψ, χ1, and χ2 potential energy surfaces for blocked dipeptides. New Fourier coefficients for the dihedral angle terms of the OPLS-AA force field were fit to these surfaces, utilizing a Boltzmann-weighted error function and systematically examining the effects of weighting temperature. To prevent overfitting to the available data, a minimal number of new residue-specific and peptide-specific torsion terms were developed. Extensive experimental solution-phase and quantum chemical gas-phase benchmarks were used to assess the quality of the new parameters, named OPLS-AA/M, demonstrating significant improvement over previous OPLS-AA force fields. A Boltzmann weighting temperature of 2000 K was determined to be optimal for fitting the new Fourier coefficients for dihedral angle parameters. Conclusions are drawn from the results for best practices for developing new torsion parameters for protein force fields.

Structure-based virtual screening for drug discovery: principles, applications and recent advances

Structure-based drug discovery (SBDD) is becoming an essential tool in assisting fast and cost-efficient lead discovery and optimization. The application of rational, structure-based drug design is proven to be more efficient than the traditional way of drug discovery since it aims to understand the molecular basis of a disease and utilizes the knowledge of the three-dimensional structure of the biological target in the process. In this review, we focus on the principles and applications of Virtual Screening (VS) within the context of SBDD and examine different procedures ranging from the initial stages of the process that include receptor and library pre-processing, to docking, scoring and post-processing of topscoring hits. Recent improvements in structure-based virtual screening (SBVS) efficiency through ensemble docking, induced fit and consensus docking are also discussed. The review highlights advances in the field within the framework of several success studies that have led to nM inhibition directly from VS and provides recent trends in library design as well as discusses limitations of the method. Applications of SBVS in the design of substrates for engineered proteins that enable the discovery of new metabolic and signal transduction pathways and the design of inhibitors of multifunctional proteins are also reviewed. Finally, we contribute two promising VS protocols recently developed by us that aim to increase inhibitor selectivity. In the first protocol, we describe the discovery of micromolar inhibitors through SBVS designed to inhibit the mutant H1047R PI3Kα kinase. Second, we discuss a strategy for the identification of selective binders for the RXRα nuclear receptor. In this protocol, a set of target structures is constructed for ensemble docking based on binding site shape characterization and clustering, aiming to enhance the hit rate of selective inhibitors for the desired protein target through the SBVS process.

Identifying and designing chemicals with minimal acute aquatic toxicity

Industrial ecology has revolutionized our understanding of material stocks and flows in our economy and society. For this important discipline to have even deeper impact, we must understand the inherent nature of these materials in terms of human health and the environment. This paper focuses on methods to design synthetic chemicals to reduce their intrinsic ability to cause adverse consequence to the biosphere. Advances in the fields of computational chemistry and molecular toxicology in recent decades allow the development of predictive models that inform the design of molecules with reduced potential to be toxic to humans or the environment. The approach presented herein builds on the important work in quantitative structure-activity relationships by linking toxicological and chemical mechanistic insights to the identification of critical physical-chemical properties needed to be modified. This in silico approach yields design guidelines using boundary values for physiochemical properties. Acute aquatic toxicity serves as a model endpoint in this study. Defining value ranges for properties related to bioavailability and reactivity eliminates 99% of the chemicals in the highest concern for acute aquatic toxicity category. This approach and its future implementations are expected to yield very powerful tools for life cycle assessment practitioners and molecular designers that allow rapid assessment of multiple environmental and human health endpoints and inform modifications to minimize hazard.

Equilibrium molecular thermodynamics from Kirkwood sampling

We present two methods for barrierless equilibrium sampling of molecular systems based on the recently proposed Kirkwood method (J. Chem. Phys.2009, 130, 134102). Kirkwood sampling employs low-order correlations among internal coordinates of a molecule for random (or non-Markovian) sampling of the high dimensional conformational space. This is a geometrical sampling method independent of the potential energy surface. The first method is a variant of biased Monte Carlo, where Kirkwood sampling is used for generating trial Monte Carlo moves. Using this method, equilibrium distributions corresponding to different temperatures and potential energy functions can be generated from a given set of low-order correlations. Since Kirkwood samples are generated independently, this method is ideally suited for massively parallel distributed computing. The second approach is a variant of reservoir replica exchange, where Kirkwood sampling is used to construct a reservoir of conformations, which exchanges conformations with the replicas performing equilibrium sampling corresponding to different thermodynamic states. Coupling with the Kirkwood reservoir enhances sampling by facilitating global jumps in the conformational space. The efficiency of both methods depends on the overlap of the Kirkwood distribution with the target equilibrium distribution. We present proof-of-concept results for a model nine-atom linear molecule and alanine dipeptide.

Determination of partial molar volumes from free energy perturbation theory

Partial molar volume is an important thermodynamic property that gives insights into molecular size and intermolecular interactions in solution. Theoretical frameworks for determining the partial molar volume (V°) of a solvated molecule generally apply Scaled Particle Theory or Kirkwood-Buff theory. With the current abilities to perform long molecular dynamics and Monte Carlo simulations, more direct methods are gaining popularity, such as computing V° directly as the difference in computed volume from two simulations, one with a solute present and another without. Thermodynamically, V° can also be determined as the pressure derivative of the free energy of solvation in the limit of infinite dilution. Both approaches are considered herein with the use of free energy perturbation (FEP) calculations to compute the necessary free energies of solvation at elevated pressures. Absolute and relative partial molar volumes are computed for benzene and benzene derivatives using the OPLS-AA force field. The mean unsigned error for all molecules is 2.8 cm(3) mol(-1). The present methodology should find use in many contexts such as the development and testing of force fields for use in computer simulations of organic and biomolecular systems, as a complement to related experimental studies, and to develop a deeper understanding of solute-solvent interactions.

Design, synthesis, and protein crystallography of biaryltriazoles as potent tautomerase inhibitors of macrophage migration inhibitory factor

Optimization is reported for biaryltriazoles as inhibitors of the tautomerase activity of human macrophage migration inhibitory factor (MIF), a proinflammatory cytokine associated with numerous inflammatory diseases and cancer. A combined approach was taken featuring organic synthesis, enzymatic assaying, crystallography, and modeling including free-energy perturbation (FEP) calculations. X-ray crystal structures for 3a and 3b bound to MIF are reported and provided a basis for the modeling efforts. The accommodation of the inhibitors in the binding site is striking with multiple hydrogen bonds and aryl-aryl interactions. Additional modeling encouraged pursuit of 5-phenoxyquinolinyl analogues, which led to the very potent compound 3s. Activity was further enhanced by addition of a fluorine atom adjacent to the phenolic hydroxyl group as in 3w, 3z, 3aa, and 3bb to strengthen a key hydrogen bond. It is also shown that physical properties of the compounds can be modulated by variation of solvent-exposed substituents. Several of the compounds are likely the most potent known MIF tautomerase inhibitors; the most active ones are more than 1000-fold more active than the well-studied (R)-ISO-1 and more than 200-fold more active than the chromen-4-one Orita-13.

Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field

Designing tight-binding ligands is a primary objective of small-molecule drug discovery. Over the past few decades, free-energy calculations have benefited from improved force fields and sampling algorithms, as well as the advent of low-cost parallel computing. However, it has proven to be challenging to reliably achieve the level of accuracy that would be needed to guide lead optimization (∼5× in binding affinity) for a wide range of ligands and protein targets. Not surprisingly, widespread commercial application of free-energy simulations has been limited due to the lack of large-scale validation coupled with the technical challenges traditionally associated with running these types of calculations. Here, we report an approach that achieves an unprecedented level of accuracy across a broad range of target classes and ligands, with retrospective results encompassing 200 ligands and a wide variety of chemical perturbations, many of which involve significant changes in ligand chemical structures. In addition, we have applied the method in prospective drug discovery projects and found a significant improvement in the quality of the compounds synthesized that have been predicted to be potent. Compounds predicted to be potent by this approach have a substantial reduction in false positives relative to compounds synthesized on the basis of other computational or medicinal chemistry approaches. Furthermore, the results are consistent with those obtained from our retrospective studies, demonstrating the robustness and broad range of applicability of this approach, which can be used to drive decisions in lead optimization.

Free energy simulations of active-site mutants of dihydrofolate reductase

This study employs hybrid quantum mechanics-molecular mechanics (QM/MM) simulations to investigate the effect of mutations of the active-site residue I14 of E. coli dihydrofolate reductase (DHFR) on the hydride transfer. Recent kinetic measurements of the I14X mutants (X = V, A, and G) indicated slower hydride transfer rates and increasingly temperature-dependent kinetic isotope effects (KIEs) with systematic reduction of the I14 side chain. The QM/MM simulations show that when the original isoleucine residue is substituted in silico by valine, alanine, or glycine (I14V, I14A, and I14G DHFR, respectively), the free energy barrier height of the hydride transfer reaction increases relative to the wild-type enzyme. These trends are in line with the single-turnover rate measurements reported for these systems. In addition, extended dynamics simulations of the reactive Michaelis complex reveal enhanced flexibility in the mutants, and in particular for the I14G mutant, including considerable fluctuations of the donor-acceptor distance (DAD) and the active-site hydrogen bonding network compared with those detected in the native enzyme. These observations suggest that the perturbations induced by the mutations partly impair the active-site environment in the reactant state. On the other hand, the average DADs at the transition state of all DHFR variants are similar. Crystal structures of I14 mutants (V, A, and G) confirmed the trend of increased flexibility of the M20 and other loops.

Computational Studies of the Effect of Shock Waves on the Binding of Model Complexes

We have simulated effects of a shock wave in water that would result from the collapse of a cavitation bubble on binding in model complexes. We have considered a benzene dimer, a pair of uracil molecules, a complex of fragments of the X-linked inhibitor of apoptosis and caspase-9, and a fragment of c-Myc oncoprotein in binding with its dimerization partner Max. The effect of the shock waves was simulated by adding a momentum to a slab of solvent water molecules and observing the system as the slab moved and caused changes. In the cases of the small molecular pairs, the passage of the shock waves lead to dissociation of the complexes. The behavior of the protein systems was more complex, yet significant disruption of the binding and geometry was also observed. In all the cases, the effects did not occur during the immediate impact of the high-momentum solvent molecules, but rather during the expansion of the compressed system that followed the passage of the waves. The rationale of the studies was in attempting to understand the strong effects that irradiation with a low-intensity ultrasound can have on biomolecular systems, because such ultrasound irradiation can cause cavitation bubbles to be produced and collapse, thus leading to local shock wave generation. The long-term objective is to contribute to future design of synergetic ultrasound and chemical drug strategy of protein inhibition.

Structure-guided development of deoxycytidine kinase inhibitors with nanomolar affinity and improved metabolic stability

Recently, we have shown that small molecule dCK inhibitors in combination with pharmacological perturbations of de novo dNTP biosynthetic pathways could eliminate acute lymphoblastic leukemia cells in animal models. However, our previous lead compound had a short half-life in vivo. Therefore, we set out to develop dCK inhibitors with favorable pharmacokinetic properties. We delineated the sites of the inhibitor for modification, guided by crystal structures of dCK in complex with the lead compound and with derivatives. Crystal structure of the complex between dCK and the racemic mixture of our new lead compound indicated that the R-isomer is responsible for kinase inhibition. This was corroborated by kinetic analysis of the purified enantiomers, which showed that the R-isomer has >60-fold higher affinity than the S-isomer for dCK. This new lead compound has significantly improved metabolic stability, making it a prime candidate for dCK-inhibitor based therapies against hematological malignancies and, potentially, other cancers.

Simulating chemical reactions in ionic liquids using QM/MM methodology

The use of ionic liquids as a reaction medium for chemical reactions has dramatically increased in recent years due in large part to the numerous reported advances in catalysis and organic synthesis. In some extreme cases, ionic liquids have been shown to induce mechanistic changes relative to conventional solvents. Despite the large interest in the solvents, a clear understanding of the molecular factors behind their chemical impact is largely unknown. This feature article reviews our efforts developing and applying mixed quantum and molecular mechanical (QM/MM) methodology to elucidate the microscopic details of how these solvents operate to enhance rates and alter mechanisms for industrially and academically important reactions, e.g., Diels-Alder, Kemp eliminations, nucleophilic aromatic substitutions, and β-eliminations. Explicit solvent representation provided the medium dependence of the activation barriers and atomic-level characterization of the solute-solvent interactions responsible for the experimentally observed "ionic liquid effects". Technical advances are also discussed, including a linear-scaling pairwise electrostatic interaction alternative to Ewald sums, an efficient polynomial fitting method for modeling proton transfers, and the development of a custom ionic liquid OPLS-AA force field.

Challenges in computational studies of enzyme structure, function and dynamics

In this review we give an overview of the field of Computational enzymology. We start by describing the birth of the field, with emphasis on the work of the 2013 chemistry Nobel Laureates. We then present key features of the state-of-the-art in the field, showing what theory, accompanied by experiments, has taught us so far about enzymes. We also briefly describe computational methods, such as quantum mechanics-molecular mechanics approaches, reaction coordinate treatment, and free energy simulation approaches. We finalize by discussing open questions and challenges.

Molecular dynamics and Monte Carlo simulations for protein-ligand binding and inhibitor design

BACKGROUND

Non-nucleoside inhibitors of HIV reverse transcriptase are an important component of treatment against HIV infection. Novel inhibitors are sought that increase potency against variants that contain the Tyr181Cys mutation.

METHODS

Molecular dynamics based free energy perturbation simulations have been run to study factors that contribute to protein-ligand binding, and the results are compared with those from previous Monte Carlo based simulations and activity data.

RESULTS

Predictions of protein-ligand binding modes are very consistent for the two simulation methods; the accord is attributed to the use of an enhanced sampling protocol. The Tyr181Cys binding pocket supports large, hydrophobic substituents, which is in good agreement with experiment.

CONCLUSIONS

Although some discrepancies exist between the results of the two simulation methods and experiment, free energy perturbation simulations can be used to rapidly test small molecules for gains in binding affinity.

GENERAL SIGNIFICANCE

Free energy perturbation methods show promise in providing fast, reliable and accurate data that can be used to complement experiment in lead optimization projects. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".

Parameter choice matters: validating probe parameters for use in mixed-solvent simulations

Probe mapping is a common approach for identifying potential binding sites in structure-based drug design; however, it typically relies on energy minimizations of probes in the gas phase and a static protein structure. The mixed-solvent molecular dynamics (MixMD) approach was recently developed to account for full protein flexibility and solvation effects in hot-spot mapping. Our first study used only acetonitrile as a probe, and here, we have augmented the set of functional group probes through careful testing and parameter validation. A diverse range of probes are needed in order to map complex binding interactions. A small variation in probe parameters can adversely effect mixed-solvent behavior, which we highlight with isopropanol. We tested 11 solvents to identify six with appropriate behavior in TIP3P water to use as organic probes in the MixMD method. In addition to acetonitrile and isopropanol, we have identified acetone, N-methylacetamide, imidazole, and pyrimidine. These probe solvents will enable MixMD studies to recover hydrogen-bonding sites, hydrophobic pockets, protein–protein interactions, and aromatic hotspots. Also, we show that ternary-solvent systems can be incorporated within a single simulation. Importantly, these binary and ternary solvents do not require artificial repulsion terms like other methods. Within merely 5 ns, layered solvent boxes become evenly mixed for soluble probes. We used radial distribution functions to evaluate solvent behavior, determine adequate mixing, and confirm the absence of phase separation. We recommend that radial distribution functions should be used to assess adequate sampling in all mixed-solvent techniques rather than the current practice of examining the solvent ratios at the edges of the solvent box.

Evaluation of CM5 Charges for Condensed-Phase Modeling

The recently developed Charge Model 5 (CM5) is tested for its utility in condensed-phase simulations. The CM5 approach, which derives partial atomic charges from Hirshfeld population analyses, provides excellent results for gas-phase dipole moments and is applicable to all elements of the periodic table. Herein, the adequacy of scaled CM5 charges for use in modeling aqueous solutions has been evaluated by computing free energies of hydration (ΔGhyd) for 42 neutral organic molecules via Monte Carlo statistical mechanics. An optimal scaling factor for the CM5 charges was determined to be 1.27, resulting in a mean unsigned error (MUE) of 1.1 kcal/mol for the free energies of hydration. Testing for an additional 20 molecules gave an MUE of 1.3 kcal/mol. The high precision of the results is confirmed by free energy calculations using both sequential perturbations and complete molecular annihilation. Performance for specific functional groups is discussed; sulfur-containing molecules yield the largest errors. In addition, the scaling factor of 1.27 is shown to be appropriate for CM5 charges derived from a variety of density functional methods and basis sets. Though the average errors from the 1.27*CM5 results are only slightly lower than those using 1.14*CM1A charges, the broader applicability and easier access to CM5 charges via the Gaussian program are additional attractive features. The 1.27*CM5 charge model can be used for an enormous variety of applications in conjunction with many fixed-charge force fields and molecular modeling programs.

A theoretical view of protein dynamics

Proteins are fascinating supramolecular structures, which are able to recognize ligands transforming binding information into chemical signals. They can transfer information across the cell, can catalyse complex chemical reactions, and are able to transform energy into work with much more efficiency than any human engine. The unique abilities of proteins are tightly coupled with their dynamic properties, which are coded in a complex way in the sequence and carefully refined by evolution. Despite its importance, our experimental knowledge of protein dynamics is still rather limited, and mostly derived from theoretical calculations. I will review here, in a systematic way, the current state-of-the-art theoretical approaches to the study of protein dynamics, emphasizing the most recent advances, examples of use and the expected lines of development in the near future.

GEM*: A Molecular Electronic Density-Based Force Field for Molecular Dynamics Simulations

GEM*, a force field that combines Coulomb and Exchange terms calculated with Hermite Gaussians with the polarization, bonded, and modified van der Waals terms from AMOEBA is presented. GEM* is tested on an initial water model fitted at the same level as AMOEBA. The integrals required for the evaluation of the intermolecular Coulomb interactions are efficiently evaluated by means of reciprocal space methods. The GEM* water model is tested by comparing energies and forces for a series of water oligomers and MD simulations. Timings for GEM* compared to AMOEBA are presented and discussed.

Cooperative effects and optimal halogen bonding motifs for self-assembling systems

Halogen bonding, due to its directionality and tunable strength, is being increasingly utilized in self-assembling materials and crystal engineering. Using density functional theory (DFT) and molecular mechanics (OPLS/CM1Ax) calculations, multiply halogen bonded complexes of brominated imidazole and pyridine are investigated along with their potential in construction of self-assembling architectures. Dimers with 1–10 halogen bonds are considered and reveal maximal binding energies of 3–36 kcal/mol. Cooperative (nonadditive) effects are found in complexes that extend both along and perpendicular to the halogen bonding axes, with interaction energies depending on polarization, secondary interactions, and ring spacers. Four structural motifs were identified to yield optimal halogen bonding. For the largest systems, the excellent agreement found between the DFT and OPLS/CM1Ax results supports the utility of the latter approach for analysis and design of self-assembling supramolecular structures.

Quantum and Molecular Mechanical (QM/MM) Monte Carlo Techniques for Modeling Condensed-Phase Reactions

A recent review (Acc. Chem. Res. 2010, 43:142-151) examined our use and development of a combined quantum and molecular mechanical (QM/MM) technique for modelling organic and enzymatic reactions. Advances included the PDDG/PM3 semiempirical QM (SQM) method, computation of multi-dimensional potentials of mean force (PMF), incorporation of on-the-fly QM in Monte Carlo simulations, and a polynomial quadrature method for rapidly treating proton-transfer reactions. The current article serves as a follow up on our progress. Highlights include new reactions, alternative SQM methods, a polarizable OPLS force field, and novel solvent environments, e.g., "on water" and room temperature ionic liquids. The methodology is strikingly accurate across a wide range of condensed-phase and antibody-catalyzed reactions including substitution, decarboxylation, elimination, isomerization, and pericyclic classes. Comparisons are made to systems treated with continuum-based solvents and ab initio or density functional theory (DFT) methods. Overall, the QM/MM methodology provides detailed characterization of reaction paths, proper configurational sampling, several advantages over implicit solvent models, and a reasonable computational cost.

Phenol-quinone tautomerism in (arylazo)naphthols and the analogous Schiff bases: benchmark calculations

Tautomerization energies of a series of isomeric [(4-R-phenyl)azo]naphthols and the analogous Schiff bases (R = N(CH3)2, OCH3, H, CN, NO2) are calculated by LPNO-CEPA/1-CBS using the def2-TZVPP and def2-QZVPP basis sets for extrapolation. The performance of various density functionals (B3LYP, M06-2X, PW6B95, B2PLYP, mPW2PLYP, PWPB95) as well as MP2 and SCS-MP2 is evaluated against these results. M06-2X and SCS-MP2 yield results close to the LPNO-CEPA/1-CBS values. Solvent effects (CCl4, CHCl3, CH3CN, and CH3OH) are treated by a variety of bulk solvation models (SM8, IEFPCM, COSMO, PBF, and SMD) as well as explicit solvation (Monte Carlo free energy perturbation using the OPLSAA force field).

Enhanced Monte Carlo Sampling through Replica Exchange with Solute Tempering

With a view to improving the consistency of free energy perturbation calculations in Monte Carlo simulations of protein–ligand complexes, we have implemented the replica exchange with solute tempering (REST) method in the MCPRO software. By augmenting the standard REST approach with regular attempted jumps in selected dihedral angles, our combined method facilitates sampling of ligand binding modes that are separated by high free energy barriers and ensures that computed free energy changes are considerably less dependent on the starting conditions and the chosen mutation pathway than those calculated with standard Monte Carlo sampling. We have applied the enhanced sampling method to the calculation of the activities of seven non-nucleoside inhibitors of HIV-1 reverse transcriptase, and its Tyr181Cys variant, and have shown that a range of binding orientations is possible depending on the nature of the ligand and the presence of mutations at the binding site.

Exploring early stages of the chemical unfolding of proteins at the proteome scale

After decades of using urea as denaturant, the kinetic role of this molecule in the unfolding process is still undefined: does urea actively induce protein unfolding or passively stabilize the unfolded state? By analyzing a set of 30 proteins (representative of all native folds) through extensive molecular dynamics simulations in denaturant (using a range of force-fields), we derived robust rules for urea unfolding that are valid at the proteome level. Irrespective of the protein fold, presence or absence of disulphide bridges, and secondary structure composition, urea concentrates in the first solvation shell of quasi-native proteins, but with a density lower than that of the fully unfolded state. The presence of urea does not alter the spontaneous vibration pattern of proteins. In fact, it reduces the magnitude of such vibrations, leading to a counterintuitive slow down of the atomic-motions that opposes unfolding. Urea stickiness and slow diffusion is, however, crucial for unfolding. Long residence urea molecules placed around the hydrophobic core are crucial to stabilize partially open structures generated by thermal fluctuations. Our simulations indicate that although urea does not favor the formation of partially open microstates, it is not a mere spectator of unfolding that simply displaces to the right of the folded←→unfolded equilibrium. On the contrary, urea actively favors unfolding: it selects and stabilizes partially unfolded microstates, slowly driving the protein conformational ensemble far from the native one and also from the conformations sampled during thermal unfolding.

The delicate equilibrium between the folded and functional structure of a protein and its unfolded state is highly dependent on environmental variables such as the solvent. For example the co-solvent urea is a well-known protein denaturant that displaces the equilibrium towards unstructured and non-functional conformations of proteins. However the molecular mechanism behind its ability remains an enigma and the interpretation of the experimental data is still ambiguous. By analyzing a set of representative proteins through extensive molecular dynamics simulations in urea, we provide a robust and consensus picture of the first stages of urea-driven protein unfolding and elucidate the role of urea in accelerating protein unfolding. Our results suggest that urea, thanks to its stickiness and slow diffusion, benefits from the intrinsic flexibility of proteins and stabilizes partially open-states, slowly driving the protein toward unfolding.

Design, synthesis, and biological evaluation of 1-phenylpyrazolo[3,4-e]pyrrolo[3,4-g]indolizine-4,6(1H,5H)-diones as new glycogen synthase kinase-3β inhibitors

Compound 5 was selected from our in-house library as a suitable starting point for the rational design of new GSK-3β inhibitors. MC/FEP calculations of 5 led to the identification of a structural class of new GSK-3β inhibitors. Compound 18 inhibited GSK-3β with an IC50 of 0.24 μM and inhibited tau phosphorylation in a cell-based assay. It proved to be a selective inhibitor of GSK-3 against a panel of 17 kinases and showed >10-fold selectivity against CDK2. Calculated physicochemical properties and Volsurf predictions suggested that compound 18 has the potential to diffuse passively across the blood-brain barrier.

Picomolar inhibitors of HIV reverse transcriptase featuring bicyclic replacement of a cyanovinylphenyl group

Members of the catechol diether class are highly potent non-nucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs). The most active compounds yield EC50 values below 0.5 nM in assays using human T-cells infected by wild-type HIV-1. However, these compounds such as rilpivirine, the most recently FDA-approved NNRTI, bear a cyanovinylphenyl (CVP) group. This is an uncommon substructure in drugs that gives reactivity concerns. In the present work, computer simulations were used to design bicyclic replacements for the CVP group. The predicted viability of a 2-cyanoindolizinyl alternative was confirmed experimentally and provided compounds with 0.4 nM activity against the wild-type virus. The compounds also performed well with EC50 values of 10 nM against the challenging HIV-1 variant that contains the Lys103Asn/Tyr181Cys double mutation in the RT enzyme. Indolyl and benzofuranyl analogues were also investigated; the most potent compounds in these cases have EC50 values toward wild-type HIV-1 near 10 nM and high-nanomolar activities toward the double-variant. The structural expectations from the modeling were much enhanced by obtaining an X-ray crystal structure at 2.88 Å resolution for the complex of the parent 2-cyanoindolizine 10b and HIV-1 RT. The aqueous solubilities of the most potent indolizine analogues were also measured to be ~40 μg/mL, which is similar to that for the approved drug efavirenz and ~1000-fold greater than for rilpivirine.

A β-peptide agonist of the GLP-1 receptor, a class B GPCR

Previous work has shown that certain β(3)-peptides can effectively mimic the side chain display of an α-helix and inhibit interactions between proteins, both in vitro and in cultured cells. Here we describe a β(3)-peptide analog of GLP-1, CC-3(Act), that interacts with the GLP-1R extracellular domain (nGLP-1R) in vitro in a manner that competes with exendin-4 and induces GLP-1R-dependent cAMP signaling in cultured CHO-K1 cells expressing GLP-1R.

QuanPol: a full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree-Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self-consistent field method, Møller-Plesset perturbation theory method, and time-dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self-consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented.

Optimization of diarylazines as anti-HIV agents with dramatically enhanced solubility

Non-nucleoside inhibitors of HIV-1 reverse transcriptase are reported that have ca. 100-fold greater solubility than the structurally related drugs etravirine and rilpivirine, while retaining high anti-viral activity. The solubility enhancements come from strategic placement of a morpholinylalkoxy substituent in the entrance channel of the NNRTI binding site. Compound 4d shows low-nanomolar activity similar to etravirine towards wild-type HIV-1 and key viral variants.

Development of new deoxycytidine kinase inhibitors and noninvasive in vivo evaluation using positron emission tomography

Combined inhibition of ribonucleotide reductase and deoxycytidine kinase (dCK) in multiple cancer cell lines depletes deoxycytidine triphosphate pools leading to DNA replication stress, cell cycle arrest, and apoptosis. Evidence implicating dCK in cancer cell proliferation and survival stimulated our interest in developing small molecule dCK inhibitors. Following a high throughput screen of a diverse chemical library, a structure-activity relationship study was carried out. Positron Emission Tomography (PET) using (18)F-L-1-(2'-deoxy-2'-FluoroArabinofuranosyl) Cytosine ((18)F-L-FAC), a dCK-specific substrate, was used to rapidly rank lead compounds based on their ability to inhibit dCK activity in vivo. Evaluation of a subset of the most potent compounds in cell culture (IC50 = ∼1-12 nM) using the (18)F-L-FAC PET pharmacodynamic assay identified compounds demonstrating superior in vivo efficacy.

QM/MM investigation on 1,3-dipolar cycloadditions of the phthalazinium dicyanomethanide with three different dipolarophiles on water and in solution

An "on water" environment, describing the reactions with insoluble reactants in water, has been reported to give high yields of products compared to organic solvents. The 1,3-dipolar cycloadditions of phthalazinium dicyanomethanide 1 with three different dipolarophiles, methyl vinyl ketone (MVK), methyl acrylate (MAC), and styrene (STY), have been investigated using QM/MM calculations in water, acetonitrile, and acetonitrile-water solvent mixtures, as well as at the vacuum-water interface. Monte Carlo statistical mechanics simulations utilizing the free-energy perturbation theory and PDDG/PM3 for the QM method have been used. The transition structures for all three reactions do not show large variations among different solvents. However, the calculated free energies of activation at the interface are found to be higher than those calculated in bulk water. Computed energy pair distributions and radial distribution functions reveal a uniform loss of hydrogen bonds for the reactants and transitions states in progressing from bulk water to the vacuum-water interface. The hydrophobic effects in the reactions of 1 with MVK and MAC are similar for both, and weaker than the effect in the reaction with STY. According to the results in water-acetonitrile mixtures at different molar ratios, it is clear that the special hydrogen bonding effects are the main reason which leads to the rapid rate enhancement in progressing from a water-acetonitrile molar ratio of 0.9 : 0.1 to pure water. New insights into solvent effects for 1,3-dipolar cycloadditions are presented herein.