John C. Shelley

Prediction of Absolute Solvation Free Energies using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field

The accurate prediction of protein-ligand binding free energies is a primary objective in computer-aided drug design. The solvation free energy of a small molecule provides a surrogate to the desolvation of the ligand in the thermodynamic process of protein-ligand binding. Here, we use explicit solvent molecular dynamics free energy perturbation to predict the absolute solvation free energies of a set of 239 small molecules, spanning diverse chemical functional groups commonly found in drugs and drug-like molecules. We also compare the performance of absolute solvation free energies obtained using the OPLS_2005 force field with two other commonly used small molecule force fields-general AMBER force field (GAFF) with AM1-BCC charges and CHARMm-MSI with CHelpG charges. Using the OPLS_2005 force field, we obtain high correlation with experimental solvation free energies (R(2) = 0.94) and low average unsigned errors for a majority of the functional groups compared to AM1-BCC/GAFF or CHelpG/CHARMm-MSI. However, OPLS_2005 has errors of over 1.3 kcal/mol for certain classes of polar compounds. We show that predictions on these compound classes can be improved by using a semiempirical charge assignment method with an implicit bond charge correction.

Towards the comprehensive, rapid, and accurate prediction of the favorable tautomeric states of drug-like molecules in aqueous solution

Generating the appropriate protonation states of drug-like molecules in solution is important for success in both ligand- and structure-based virtual screening. Screening collections of millions of compounds requires a method for determining tautomers and their energies that is sufficiently rapid, accurate, and comprehensive. To maximise enrichment, the lowest energy tautomers must be determined from heterogeneous input, without over-enumerating unfavourable states. While computationally expensive, the density functional theory (DFT) method M06-2X/aug-cc-pVTZ(-f) [PB-SCRF] provides accurate energies for enumerated model tautomeric systems. The empirical Hammett–Taft methodology can very rapidly extrapolate substituent effects from model systems to drug-like molecules via the relationship between pKT and pKa. Combining the two complementary approaches transforms the tautomer problem from a scientific challenge to one of engineering scale-up, and avoids issues that arise due to the very limited number of measured pKT values, especially for the complicated heterocycles often favoured by medicinal chemists for their novelty and versatility. Several hundreds of pre-calculated tautomer energies and substituent pKa effects are tabulated in databases for use in structural adjustment by the program Epik, which treats tautomers as a subset of the larger problem of the protonation states in aqueous ensembles and their energy penalties. Accuracy and coverage is continually improved and expanded by parameterizing new systems of interest using DFT and experimental data. Recommendations are made for how to best incorporate tautomers in molecular design and virtual screening workflows.

Contribution of explicit solvent effects to the binding affinity of small-molecule inhibitors in blood coagulation factor serine proteases.

The prevention of blood coagulation is important in treating thromboembolic disorders, and several serine proteases involved in the coagulation cascade have been classified as pharmaceutically relevant. Whereas structure‐based drug design has contributed to the development of some serine protease inhibitors, traditional computational methods have not been able to fully describe structure–activity relationships (SAR). Here, we study the SAR for a number of serine proteases by using a method that calculates the thermodynamic properties (enthalpy and entropy) of the water that solvates the active site. We show that the displacement of water from specific subpockets (such as S1–4 and the ester binding pocket) of the active site by the ligand can govern potency, especially for cases in which small chemical changes (i.e., a methyl group or halogen) result in a substantial increase in potency. Furthermore, we describe how relative binding free energies can be estimated by combining the water displacement energy with complementary terms from an implicit solvent molecular mechanics description binding.

Computer simulation studies of biomembranes using a coarse grain model

A computationally efficient coarse grain (CG) model designed to mimic the lipid molecule, dimyristoylphosphatidylcholine (DMPC) is used to study the self-assembly of a lamellar bilayer starting from a disordered configuration. The utility of the CG model is illustrated by examining structural and dynamical properties of a hydrated bilayer system containing 1024 lipid molecules. Comparisons with results for an all-atom model of DMPC suggest that the CG model is about four orders of magnitude less demanding of CPU time.

Modeling Local Structural Rearrangements Using FEP/REST: Application to Relative Binding Affinity Predictions of CDK2 Inhibitors.

Accurate and reliable calculation of protein-ligand binding affinities remains a hotbed of computer-aided drug design research. Despite the potentially large impact FEP (free energy perturbation) may have in drug design projects, practical applications of FEP in industrial contexts have been limited. In this work, we use a recently developed method, FEP/REST (free energy perturbation/replica exchange with solute tempering), to calculate the relative binding affinities for a set of congeneric ligands binding to the CDK2 receptor. We compare the FEP/REST results with traditional FEP/MD (molecular dynamics) results and MM/GBSA (molecular mechanics/Generalized Born Surface Area model) results and examine why FEP/REST performed notably better than these other methods, as well as why certain ligand mutations lead to large increases of the binding affinity while others do not. We also introduce a mathematical framework for assessing the consistency and reliability of the calculations using cycle closures in FEP mutation paths.

Ligand-Dependent Activation and Deactivation of the Human Adenosine A2A Receptor

G-protein-coupled receptors (GPCRs) are membrane proteins with critical functions in cellular signal transduction, representing a primary class of drug targets. Acting by direct binding, many drugs modulate GPCR activity and influence the signaling pathways associated with numerous diseases. However, complete details of ligand-dependent GPCR activation/deactivation are difficult to obtain from experiments. Therefore, it remains unclear how ligands modulate a GPCRs activity. To elucidate the ligand-dependent activation/deactivation mechanism of the human adenosine A2A receptor (AA2AR), a member of the class A GPCRs, we performed large-scale unbiased molecular dynamics and metadynamics simulations of the receptor embedded in a membrane. At the atomic level, we have observed distinct structural states that resemble the active and inactive states. In particular, we noted key structural elements changing in a highly concerted fashion during the conformational transitions, including six conformational states of a tryptophan (Trp246(6.48)). Our findings agree with a previously proposed view that, during activation, this tryptophan residue undergoes a rotameric transition that may be coupled to a series of coherent conformational changes, resulting in the opening of the G-protein binding site. Further, metadynamics simulations provide quantitative evidence for this mechanism, suggesting how ligand binding shifts the equilibrium between the active and inactive states. Our analysis also proposes that a few specific residues are associated with agonism/antagonism, affinity, and selectivity, and suggests that the ligand-binding pocket can be thought of as having three distinct regions, providing dynamic features for structure-based design. Additional simulations with AA2AR bound to a novel ligand are consistent with our proposed mechanism. Generally, our study provides insights into the ligand-dependent AA2AR activation/deactivation in addition to what has been found in crystal structures. These results should aid in the discovery of more effective and selective GPCR ligands.

Self-assembly of a phospholipid Langmuir monolayer using coarse-grained molecular dynamics simulations

Molecular dynamics simulations using a coarse-grained (CG) model for dimyristoyl-phosphatidyl-choline and water molecules have been carried out to follow the self-assembly process of a Langmuir monolayer. We expand on a previous study of the characteristics of the CG model where we compare the rotational and translational constants of the present model to those of an all-atom (AA) model, and find that the rotational and translational timescales are up to two orders of magnitude faster than in an AA model. We then apply the model to the self-assembly of a Langmuir monolayer. The initial randomly distributed system, which consists of 80 lipids and 5000 water sites, quickly self-assembles into two Langmuir monolayers and a micelle in the bulk water region. The micelle slowly diffuses towards and fuses with one of the interfacial monolayers, leaving the final equilibrated state with a Langmuir monolayer at each of the two air/water interfaces. The effective speed-up gained from the CG approach gives access to timescales and spatial scales that are much larger than those currently accessible with AA models.

Macrocycle Conformational Sampling with MacroModel

Sampling low energy conformations of macrocycles is challenging due to the large size of many of these molecules and the constraints imposed by the macrocycle. We present a new conformational search method (implemented in MacroModel) that uses brief MD simulations followed by minimization and normal-mode search steps. The method was parametrized using a set of 100 macrocycles from the PDB and CSD. It was then tested on a publicly available data set for which there are published results using alternative methods; we found that when the same force field is used (in this case MMFFs in vacuum), our method tended to identify conformations with lower energies than what the other methods identified. The performance on a new set of 50 macrocycles from the PDB and CSD was also quite good; the mean and median RMSD values for just the ring atoms were 0.60 and 0.33 Å, respectively. However, the RMSD values for macrocycles with more than 30 ring-atoms were quite a bit larger compared to the smaller macrocycles. Possible origins for this and ideas for improving the performance on very large macrocycles are discussed.

Melittin Aggregation in Aqueous Solutions: Insight from Molecular Dynamics Simulations

Melittin is a natural peptide that aggregates in aqueous solutions with paradigmatic monomer-to-tetramer and coil-to-helix transitions. Since little is known about the molecular mechanisms of melittin aggregation in solution, we simulated its self-aggregation process under various conditions. After confirming the stability of a melittin tetramer in solution, we observed—for the first time in atomistic detail—that four separated melittin monomers aggregate into a tetramer. Our simulated dependence of melittin aggregation on peptide concentration, temperature, and ionic strength is in good agreement with prior experiments. We propose that melittin mainly self-aggregates via a mechanism involving the sequential addition of monomers, which is supported by both qualitative and quantitative evidence obtained from unbiased and metadynamics simulations. Moreover, by combining computer simulations and a theory of the electrical double layer, we provide evidence to suggest why melittin aggregation in solution likely stops at the tetramer, rather than forming higher-order oligomers. Overall, our study not only explains prior experimental results at the molecular level but also provides quantitative mechanistic information that may guide the engineering of melittin for higher efficacy and safety.

A New Mixed All-Atom/Coarse-Grained Model: Application to Melittin Aggregation in Aqueous Solution

We introduce a new mixed resolution, all-atom/coarse-grained approach (AACG), for modeling peptides in aqueous solution and apply it to characterizing the aggregation of melittin. All of the atoms in peptidic components are represented, while a single site is used for each water molecule. With the full flexibility of the peptide retained, our AACG method achieves speedups by a factor of 3–4 for CPU time reduction and another factor of roughly 7 for diffusion. An Ewald treatment permits the inclusion of long-range electrostatic interactions. These characteristics fit well with the requirements for studying peptide association and aggregation, where the system sizes and time scales require considerable computational resources with all-atom models. In particular, AACG is well suited for biologics since changes in peptide shape and long-range electrostatics may play an important role. The application of AACG to melittin, a 26-residue peptide with a well-known propensity to aggregate in solution, serves as an initial demonstration of this technology for studying peptide aggregation. We observed the formation of melittin aggregates during our simulations and characterized the time-evolution of aggregate size distribution, buried surface areas, and residue contacts. Key interactions including π-cation and π-stacking involving TRP19 were also examined. Our AACG simulations demonstrated a clear salt effect and a moderate temperature effect on aggregation and support the molten globule model of melittin aggregates. As a showcase, this work illustrates the useful role for AACG in investigations of peptide aggregation and its potential to guide formulation and design of biologics.