Accounting for receptor flexibility, and enhanced sampling methods in computer aided drug designSinko, W., S. Lindert, J.A. McCammonChem. Biol. Drug Disc., Vol. 81, Issue 1, pp. 41−49 (2013) [PubMed 23253130]Protein flexibility plays a major role in biomolecular recognition. In
many cases, it is not obvious how molecular structure will change upon
association with other molecules. In proteins, these changes can be major,
with large deviations in overall backbone structure, or they can be more
subtle as in a side-chain rotation. Either way the algorithms that predict
the favorability of biomolecular association require relatively accurate
predictions of the bound structure to give an accurate assessment of the
energy involved in association. Here, we review a number of techniques that
have been proposed to accommodate receptor flexibility in the simulation of
small molecules binding to protein receptors. We investigate modifications
to standard rigid receptor docking algorithms and also explore enhanced
sampling techniques, and the combination of free energy calculations and
enhanced sampling techniques. The understanding and allowance for receptor
flexibility are helping to make computer simulations of ligand protein
binding more accurate. These developments may help improve the efficiency
of drug discovery and development. Efficiency will be essential as we begin
to see personalized medicine tailored to individual patients, which means
specific drugs are needed for each patient's genetic makeup.
Molecular Recognition and Ligand AssociationBaron, R., J. A.McCammonAnn. Revs. Phys. Chem., Vol. 64, pp. 151−176 (2013) [PubMed 23473376]We review recent developments in our understanding of molecular recognition
and ligand association, focusing on two major viewpoints: First, studies that
highlight new physical insight into the molecular recognition process and the
driving forces determining thermodynamic signatures of binding. Second, recent
methodological advances in applications to protein-ligand binding. In particular,
we highlight the challenges posed by compensating enthalpic and entropic terms,
competing solute and solvent contributions, and the relevance of complex
configurational ensembles composed of multiple protein, ligand, and solvent
intermediate states. As more complete physics is taken into account,
computational approaches increase their ability to complement experimental
measurements, by providing a microscopic, dynamic view of ensemble-averaged
experimental observables. Physics-based approaches are increasingly expanding
their power in pharmacology applications.
Accelerated Molecular Dynamics: An efficient enhanced sampling method for biomolecular simulationsPierce, L.C.T., W. Sinko, J.A. McCammonRamachandran Symposium Book (2013)
On the Role of Dewetting Transitions in Host-Guest Binding Free Energy CalculationsRogers, K., J.M. Ortiz Sánchez, R. Baron, M. Fajer, C.A.F. de Oliveira, J.A. McCammonJ. Chem. Theory Comput.,Vol. 9, Issue 1, pp. 46-53 (2013) [PubMed 23316123]We use thermodynamic integration (TI) and explicit solvent molecular dynamics (MD)
simulation to estimate the absolute free energy of host-guest binding. In the
unbound state, water molecules visit all of the internally accessible volume of the
host, which is fully hydrated on all sides. Upon binding of an apolar guest, the
toroidal host cavity is fully dehydrated; thus, during the intermediate λ stages
along the integration, the hydration of the host fluctuates between hydrated and
dehydrated states. Estimating free energies by TI can be especially challenging when
there is a considerable difference in hydration between the two states of interest.
We investigate these aspects using the popular TIP3P and TIP4P water models. TI free
energy estimates through MD largely depend on water-related interactions, and water
dynamics significantly affect the convergence of binding free energy calculations.
Our results indicate that wetting/dewetting transitions play a major role in slowing
the convergence of free energy estimation. We employ two alternative approaches-one
analytical and the other empirically based on actual MD sampling-to correct for the
standard state free energy. This correction is sizable (up to 4 kcal/mol), and the
two approaches provide corrections that differ by about 1 kcal/mol. For the system
considered here, the TIP4P water model combined with an analytical correction for
the standard state free energy provides higher overall accuracy. This observation
might be transferable to other systems in which water-related contributions dominate
the binding process.
Structural insight into the separate roles of IP4 and DAD in activation of histone deacetylase 3Arrar, M., R. Turnham, L. Pierce, C.A.F. de Oliveira, J.A. McCammonProtein Sci., Vol. 22, Issue 1, pp. 83-92 (2013) [PubMed 23139175]Histone deacetylases (HDACs) repress transcription by deacetylating acetyllysines
on specific histone tails. HDAC3 is implicated in neurodegenerative diseases, certain
leukemias, and even in disrupting HIV-1 latency. A recent crystal structure of HDAC3
in complex with the deacetylase-activating domain (DAD) of its corepressor complex
revealed an inositol tetraphosphate (IP4) molecule at the protein-protein interface.
IP4 was shown to play an important, yet mechanistically ambiguous, role in the activity
of HDAC3. The goal of this article is to explore the conformational ensemble of HDAC3
in its inactive apo state and in the presence of each or both of DAD and IP4. Using
triplicate, 100 ns molecular dynamic simulations, we study the apo, ternary, and
intermediate DAD-bound or IP4-bound HDAC3 states. We find that a population-shift
effect is induced by the presence of each corepressor, and is most notable in the
presence of both. Our results offer new insights into the change in dynamics necessary
for the activation of HDAC3 and highlight the roles of IP4 and DAD in this process.
Antibacterial drug leads targeting isoprenoid biosynthesisZhu, W., Y. Zhang, W. Sinko, M. Hensler, J. Olson, K.J. Molohon, S. Lindert, R. Cao, K. Li, K. Wang, Y. Wang, Y.L. Liu, A. Sankovsky, C.A.F. de Oliveira, D.A. Mitchell, V. Nizet, J.A. McCammon, E. Oldfield.Proc. Natl. Acad. Sci. USA, Vol. 110, Issue 1, pp. 123−128 (2013) [PubMed 23248302]With the rise in resistance to antibiotics such as methicillin, there is a need
for new drugs. We report here the discovery and X-ray crystallographic structures of
10 chemically diverse compounds (benzoic, diketo, and phosphonic acids, as well as
a bisamidine and a bisamine) that inhibit bacterial undecaprenyl diphosphate
synthase, an essential enzyme involved in cell wall biosynthesis. The inhibitors
bind to one or more of the four undecaprenyl diphosphate synthase inhibitor binding
sites identified previously, with the most active leads binding to site 4, outside
the catalytic center. The most potent leads are active against Staphylococcus aureus
[minimal inhibitory concentration (MIC)(90) ∼0.25 µg/mL], and one potently synergizes
with methicillin (fractional inhibitory concentration index = 0.25) and is protective
in a mouse infection model. These results provide numerous leads for antibacterial
development and open up the possibility of restoring sensitivity to drugs such as
methicillin, using combination therapies.
w-REXAMD: A Hamiltonian replica exchange approach to improve free energy calculations for systems with kinetically-trapped conformationsArrar, M., M. Fajer, W. Sinko, C. de Oliveira, J.A. McCammonJ Chem Theory Comput., Vol. 9, Issue 1, pp. 18−23 (2013) [PubMed 23316122]Free energy governs the equilibrium extent of many biological processes.
High barriers separating free energy minima often limit the sampling in
molecular dynamics (MD) simulations, leading to inaccurate free energies.
Here, we demonstrate enhanced sampling and improved free energy calculations,
relative to conventional MD, using windowed accelerated MD within a Hamiltonian
replica exchange framework (w-REXAMD). We show that for a case in which multiple
conformations are separated by large free energy barriers, w-REXAMD is a useful enhanced
sampling technique, without any necessary reweighting.
Solvent fluctuations in hydrophobic cavity-ligand binding kineticsSetny, P., R. Baron, P. Kekenes-Huskey, J.A. McCammon, J. DzubiellaProc. Natl. Acad. Sci. USA, Vol. 110, Issue 4, pp. 1197−1202 (2013) [PubMed 23297241]Water plays a crucial part in virtually all protein-ligand binding processes
in and out of equilibrium. Here, we investigate the role of water in the binding
kinetics of a ligand to a prototypical hydrophobic pocket by explicit-water
molecular dynamics (MD) simulations and implicit diffusional approaches. The
concave pocket in the unbound state exhibits wet/dry hydration oscillations whose
magnitude and time scale are significantly amplified by the approaching ligand.
In turn, the ligand's stochastic motion intimately couples to the slow hydration
fluctuations, leading to a sixfold-enhanced friction in the vicinity of the pocket
entrance. The increased friction considerably decelerates association in the otherwise
barrierless system, indicating the importance of molecular-scale hydrodynamic effects
in cavity-ligand binding arising due to capillary fluctuations. We derive and analyze
the diffusivity profile and show that the mean first passage time distribution from
the MD simulation can be accurately reproduced by a standard Brownian dynamics simulation
if the appropriate position-dependent friction profile is included. However, long-time
decays in the water-ligand (random) force autocorrelation demonstrate violation of the
Markovian assumption, challenging standard diffusive approaches for rate prediction.
Remarkably, the static friction profile derived from the force correlations strongly
resembles the profile derived on the Markovian assumption apart from a simple shift in
space, which can be rationalized by a time-space retardation in the ligand's downhill
dynamics toward the pocket. The observed spatiotemporal hydrodynamic coupling may be of
biological importance providing the time needed for conformational receptor-ligand
adjustments, typical of the induced-fit paradigm.
Fluoroketone Inhibition of Ca++-Independent Phospholipase A2 through Binding Pocket Association Defined by Hydrogen/Deuterium Exchange and Molecular Dynamics.Hsu, Y.H., D. Bucher, J. Cao, S. Li, S.W. Yang, G. Kokotos, V. Woods, J. McCammon, E. DennisJ. Amer. Chem. Soc., Vol. 135, Issue 4, pp. 1330−1337 (2013) [PubMed 23256506]The mechanism of inhibition of Group VIA Ca
2+-independent phospholipase A
2 (iPLA
2) by
fluoroketone (FK) ligands is examined using a combination of deuterium exchange mass
spectrometry (DXMS) and molecular dynamics (MD). Models for iPLA
2 were built by
homology with the known structure of patatin and equilibrated by extensive MD
simulations. Empty pockets were identified during the simulations and studied for
their ability to accommodate FK inhibitors. Ligand docking techniques showed that
the potent inhibitor 1,1,1,3-tetrafluoro-7-phenylheptan-2-one (PHFK) forms favorable
interactions inside an active site pocket where it blocks the entrance of phospholipid
substrates. The polar fluoroketone head group is stabilized by hydrogen bonds with
residues Gly486, Gly487, and Ser519. The nonpolar aliphatic chain and aromatic group
are stabilized by hydrophobic contacts with Met544, Val548, Phe549, Leu560, and Ala640.
The binding mode is supported by DXMS experiments showing an important decrease of
deuteration in the contact regions in the presence of the inhibitor. The discovery of
the precise binding mode of FK ligands to the iPLA
2sub> should greatly improve our ability
to design new inhibitors with higher potency and selectivity.
Evaluation of Hydration Free Energy by the Level-Set Variational Implicit-Solvent Model with the Coulomb-Field ApproximationGuo, Z., B. Li, J. Dzubiella, L.T. Cheng, J.A. McCammon, J. CheJ. Chem. Theory Comput., Vol. 9, Issue 3, pp. 1778−1787 (2013) [PubMed 23505345]In this article, we systematically apply a novel implicit-solvent model, the variational
implicit-solvent model (VISM) together with the Coulomb-Field Approximation (CFA), to
calculate the hydration free energy of a large set of small organic molecules. Because
these molecules have been studied in detail by molecular dynamics simulations and other
implicit-solvent models, they provide a good benchmark for evaluating the performance of
VISM-CFA. With all-atom Amber force field parameters, VISM-CFA is able to reproduce well
not only the experimental and MD simulated total hydration free energy but also the polar
and nonpolar contributions individually. The correlation between VISM-CFA and experiments
is R 2 = 0.763 for the total hydration free energy, with a root-mean-square deviation (RMSD)
of 1.83 kcal/mol, and the correlation to results from TIP3P explicit water MD simulations is
R 2 = 0.839 with a RMSD = 1.36 kcal/mol. In addition, we demonstrate that VISM captures
dewetting phenomena in the p53/MDM2 complex and hydrophobic characteristics in the system.
This work demonstrates that the level-set VISM-CFA can be used to study the energetic
behavior of realistic molecular systems with complicated geometries in solvation, protein-ligand
binding, protein-protein association, and protein folding processes.
Structure-based Network Analysis of an Evolved G-Protein Coupled Receptor Homodimer InterfaceNichols, S.E., C.X. Hernández, Y. Wang, J.A. McCammonProtein Science, Volume 22, Issue 6, 745-754 (2013) [PubMed 23553730]Crystallographic structures and experimental assays of human CXC chemokine receptor type 4
(CXCR4) provide strong evidence for the capacity to homodimerize, potentially as a means of
allosteric regulation. Even so, how this homodimer forms and its biological significance has
yet to be fully characterized. By applying principles from network analysis, sequence-based
approaches such as statistical coupling analysis to determine coevolutionary residues, can
be used in conjunction with molecular dynamics simulations to identify residues relevant to
dimerization. Here, the predominant coevolution sector lies along the observed dimer interface,
suggesting functional relevance. Furthermore, coevolution scoring provides a basis for
determining significant nodes, termed hubs, in the network formed by residues found along the
interface of the homodimer. These node residues coincide with hotspots indicating potential
druggability. Drug design efforts targeting such key residues could potentially result in
modulation of binding and therapeutic benefits for disease states, such as lung cancers,
lymphomas and latent HIV-1 infection. Furthermore, this method may be applied to any protein-protein
interaction.
Farnesyl diphosphate synthase inhibitors from in silico screeningLindert, S., W. Zhu, Y.L. Liu, R. Pang, E. Oldfield, J.A. McCammonChem. Biol. Drug Des. 81(6): 742-748 (2013) [PubMed PMC3671582]The relaxed complex scheme is an in silico drug screening method that accounts for receptor
flexibility by using molecular dynamics simulations. Here, we used this approach combined
with similarity searches and experimental inhibition assays to identify several low micro-molar,
non-bisphosphonate inhibitors, bisamidines, of farnesyl diphosphate synthase (FPPS), an enzyme
targeted by some anti-cancer and anti-microbial agents and for the treatment of bone resorption
diseases. This novel class of FPPS inhibitors have more drug-like properties than existing bisphosphonate
inhibitors, making them interesting pharmaceutical leads.
Effect of histidine protonation and rotameric states on virtual screening of M. tuberculosis RmlCKim, M.O., S.E. Nichols, Y. Wang, J.A. McCammonJ. Comp. Aid. Mol. Des., Vol. 27, Issue 3, pp. 235−246 (2013) [PubMed 23579613]While it is well established that protonation and tautomeric states of ligands can significantly
affect the results of virtual screening, such effects of ionizable residues of protein receptors
are less well understood. In this study, we focus on histidine protonation and rotameric states and
their impact on virtual screening of Mycobacterium tuberculosis enzyme RmlC. Depending on the net
charge and the location of proton(s), a histidine can adopt three states: HIP (+1 charged, both δ-
and ε-nitrogens protonated), HID (neutral, δ-nitrogen protonated), and HIE (neutral, ε-nitrogen protonated).
Due to common ambiguities in X-ray crystal structures, a histidine may also be resolved as three additional
states with its imidazole ring flipped. Here, we systematically investigate the predictive power of 36
receptor models with different protonation and rotameric states of two histidines in the RmlC active site
by using results from a previous high-throughput screening. By measuring enrichment factors and area under
the receiver operating characteristic curves, we show that virtual screening results vary depending on
hydrogen bonding networks provided by the histidines, even in the cases where the ligand does not obviously
interact with the side chain. Our results also suggest that, even with the help of widely used pKa
prediction software, assigning histidine protonation and rotameric states for virtual screening can
still be challenging and requires further examination and systematic characterization of the receptor-ligand
complex.
AutoGrow 3.0: An Improved Algorithm for Chemically Tractable, Semi-Automated Protein Inhibitor DesignJacob D. Durrant, Steffen Lindert, J. Andrew McCammonJournal of Molecular Graphics and Model., Volume 44, 104-112 (2013) [PubMed PMC3842281]We here present an improved version of AutoGrow (version 3.0), an evolutionary algorithm that works in conjunction with existing open-source software to automatically optimize candidate ligands for predicted binding affinity and other druglike properties. Though no substitute for the medicinal chemist, AutoGrow 3.0, unlike its predecessors, attempts to introduce some chemical intuition into the automated optimization process. AutoGrow 3.0 uses the rules of click chemistry to guide optimization, greatly enhancing synthesizability. Additionally, the program discards any growing ligand whose physical and chemical properties are not druglike. By carefully crafting chemically feasible druglike molecules, we hope that AutoGrow 3.0 will help supplement the chemist's efforts.
To demonstrate the utility of the program, we used AutoGrow 3.0 to generate predicted inhibitors of three important drug targets: T. brucei RNA editing ligase 1, peroxisome proliferator-activated receptor γ, and dihydrofolate reductase. In all cases, AutoGrow generated druglike molecules with high predicted binding affinities.
AutoGrow 3.0 is available free of charge (
autogrow.ucsd.edu) under the terms of the GNU General Public License and has been tested on Linux and Mac OS X
Correlated Motions and Residual Frustration in ThrombinFuglestad, B., P.M. Gasper, J.A. McCammon, P.R.L. Markwick, E.A. KomivesJ Phys Chem B. (Peter Wolynes Festschrift) 117(42):12857-63 (2013) [PubMed 23621631]Thrombin is the central protease in the cascade of blood coagulation proteases. The structure of
thrombin consists of a double β-barrel core surrounded by connecting loops and helices. Compared
to chymotrypsin, thrombin has more extended loops that are thought to have arisen from insertions
in the serine protease that evolved to impart greater specificity. Previous experiments showed
thermodynamic coupling between ligand binding at the active site and distal exosites. We present
a combined approach of molecular dynamics (MD), accelerated molecular dynamics (AMD), and analysis
of the residual local frustration of apo-thrombin and active site-bound (PPACK-thrombin). Community
analysis of the MD ensembles identified changes upon active site occupation in groups of residues
linked through correlated motions and physical contacts. AMD simulations, calibrated on measured
residual dipolar couplings, reveal that upon active site ligation correlated loop motions are
quenched, but new ones connecting the active site with distal sites where allosteric regulators
bind, emerge. Residual local frustration analysis reveals a striking correlation between frustrated
contacts and regions undergoing slow timescale dynamics. The results elucidate a motional network
that probably evolved through retention of frustrated contacts to provide facile conversion between
ensembles of states.
Discovery of Staphylococcus aureus Sortase A Inhibitors Using Virtual Screening and the Relaxed Complex Scheme.Chan, A.H., J. Wereszczynski, B.R. Amer, S.W. Yi, M.E. Jung, J.A. McCammon, R.T. ClubbChem. Biol. Drug Disc. 82, 418-428 (2013) [PubMed 23701677]Staphylococcus aureus is the leading cause of hospital-acquired infections in the United States. The emergence of multidrug-resistant strains of
S. aureus has created an urgent need for new antibiotics.
Staphylococcus aureus uses the sortase A enzyme to display surface virulence factors suggesting that compounds that inhibit its activity will function as potent anti-infective agents. Here, we report the identification of several inhibitors of sortase A using virtual screening methods that employ the relaxed complex scheme, an advanced computer-docking methodology that accounts for protein receptor flexibility. Experimental testing validates that several compounds identified in the screen inhibit the activity of sortase A. A lead compound based on the 2-phenyl-2,3-dihydro-1
H-perimidine scaffold is particularly promising, and its binding mechanism was further investigated using molecular dynamics simulations and conducting preliminary structure–activity relationship studies.
Assessing the Two-body Diffusion Tensor Calculated by Bead ModelsNuo Wang, Gary A. Huber, J. Andrew McCammonJ. Chem. Phys. 138(20), article 204117 (2013) [PubMed PMC3683057]The diffusion tensor of complex macromolecules in Stokes flow is often approximated by the bead models. The bead models are known to reproduce the experimental diffusion coefficients of a single macromolecule, but the accuracy of their calculation of the whole multi-body diffusion tensor, which is important for Brownian dynamics simulations, has not been closely investigated. As a first step, we assess the accuracy of the bead model calculated diffusion tensor of two spheres. Our results show that the bead models produce very accurate diffusion tensors for two spheres where a reasonable number of beads are used and there is no bead overlap.
Population Based Reweighting of Scaled Molecular DynamicsWilliam Sinko, Yinglong Miao, César Augusto F. de Oliveira, and J. Andrew McCammonJ. Phys. Chem. B, 117(42): 12759-12768 (2013) [PubMed PMC3808002]Molecular dynamics simulation using enhanced sampling methods is one of the powerful computational tools used to explore protein conformations and free energy landscapes. Enhanced sampling methods often employ either an increase in temperature or a flattening of the potential energy surface to rapidly sample phase space, and a corresponding reweighting algorithm is used to recover the Boltzmann statistics. However, potential energies of complex biomolecules usually involve large fluctuations on a magnitude of hundreds of kcal/mol despite minimal structural changes during simulation. This leads to noisy reweighting statistics and complicates the obtainment of accurate final results. To overcome this common issue in enhanced conformational sampling, we propose a scaled molecular dynamics method, which modifies the biomolecular potential energy surface and employs a reweighting scheme based on configurational populations. Statistical mechanical theory is applied to derive the reweighting formula, and the canonical ensemble of simulated structures is recovered accordingly. Test simulations on alanine dipeptide and the fast folding polypeptide Chignolin exhibit sufficiently enhanced conformational sampling and accurate recovery of free energy surfaces and thermodynamic properties. The results are comparable to long conventional molecular dynamics simulations and exhibit better recovery of canonical statistics over methods which employ a potential energy term in reweighting.
Phase-Field Approach to Implicit Solvation of Biomolecules with Coulomb-Field Approximation.Zhao, Y., Y.Y. Kwan, J. Che, B Li, and J.A. McCammonJ. Chem. Phys. 139, article 024111 (2013). [PubMed PMC3724799]A phase-field variational implicit-solvent approach is developed for the solvation of charged
molecules. The starting point of such an approach is the representation of a solute-solvent
interface by a phase field that takes one value in the solute region and another in the solvent
region, with a smooth transition from one to the other on a small transition layer. The
minimization of an effective free-energy functional of all possible phase fields determines
the equilibrium conformations and free energies of an underlying molecular system. All the
surface energy, the solute-solvent van der Waals interaction, and the electrostatic
interaction are coupled together self-consistently through a phase field. The surface energy
results from the minimization of a double-well potential and the gradient of a field. The
electrostatic interaction is described by the Coulomb-field approximation. Accurate and
efficient methods are designed and implemented to numerically relax an underlying charged
molecular system. Applications to single ions, a two-plate system, and a two-domain protein
reveal that the new theory and methods can capture capillary evaporation in hydrophobic
confinement and corresponding multiple equilibrium states as found in molecular dynamics
simulations. Comparisons of the phase-field and the original sharp-interface variational
approaches are discussed.
Comparing Neural-Network Scoring Functions and the State of the Art: Applications to Common Library Screening.Durrant, J., A.J. Friedman, K.E. Rogers, J.A. McCammonJ. Chem. Info. Model. 53, 1726–1735 (2013). [PubMed PMC3735370]We compare established docking programs, AutoDock Vina and Schrödinger’s Glide, to the recently published NNScore scoring functions. As expected, the best protocol to use in a virtual-screening project is highly dependent on the target receptor being studied. However, the mean screening performance obtained when candidate ligands are docked with Vina and rescored with NNScore 1.0 is not statistically different than the mean performance obtained when docking and scoring with Glide. We further demonstrate that the Vina and NNScore docking scores both correlate with chemical properties like small-molecule size and polarizability. Compensating for these potential biases leads to improvements in virtual screen performance. Composite NNScore-based scoring functions suited to a specific receptor further improve performance. We are hopeful that the current study will prove useful for those interested in computer-aided drug design.
Activation and dynamic network of the M2 muscarinic receptor.Miao, Y., S.E. Nichols, P.M. Gasper, V.T. Metzger, J.A. McCammonProc. Natl. Acad. Sci. USA 110, 10982–10987 (2013). [PubMed PMC3703993]G-protein-coupled receptors (GPCRs) mediate cellular responses to various hormones and neurotransmitters and are important targets for treating a wide spectrum of diseases. Although significant advances have been made in structural studies of GPCRs, details of their activation mechanism remain unclear. The X-ray crystal structure of the M2 muscarinic receptor, a key GPCR that regulates human heart rate and contractile forces of cardiomyocytes, was determined recently in an inactive antagonist-bound state. Here, activation of the M2 receptor is directly observed via accelerated molecular dynamics simulation, in contrast to previous microsecond-timescale conventional molecular dynamics simulations in which the receptor remained inactive. Receptor activation is characterized by formation of a Tyr206(5.58)-Tyr440(7.53) hydrogen bond and ∼6-Å outward tilting of the cytoplasmic end of transmembrane α-helix 6, preceded by relocation of Trp400(6.48) toward Phe195(5.47) and Val199(5.51) and flipping of Tyr430(7.43) away from the ligand-binding cavity. Network analysis reveals that communication in the intracellular domains is greatly weakened during activation of the receptor. Together with the finding that residue motions in the ligand-binding and G-protein-coupling sites of the apo receptor are correlated, this result highlights a dynamic network for allosteric regulation of the M2 receptor activation.
Iterative Molecular Dynamics – Rosetta Protein Structure Refinement Protocol to Improve Model Quality.Lindert, S., J. Meiler, J.A. McCammonJ. Chem. Theory Comp. 9, 3843−3847 (2013). [PubMed PMC3744128]Rosetta is one of the prime tools for high resolution protein structure refinement. While its scoring function can distinguish native-like from non-native-like conformations in many cases, the method is limited by conformational sampling for larger proteins, that is, leaving a local energy minimum in which the search algorithm may get stuck. Here, we test the hypothesis that iteration of Rosetta with an orthogonal sampling and scoring strategy might facilitate exploration of conformational space. Specifically, we run short molecular dynamics (MD) simulations on models created by de novo folding of large proteins into cryoEM density maps to enable sampling of conformational space not directly accessible to Rosetta and thus provide an escape route from the conformational traps. We present a combined MD–Rosetta protein structure refinement protocol that can overcome some of these sampling limitations. Two of four benchmark proteins showed incremental improvement through all three rounds of the iterative refinement protocol. Molecular dynamics is most efficient in applying subtle but important rearrangements within secondary structure elements and is thus highly complementary to the Rosetta refinement, which focuses on side chains and loop regions.
Insertion of the Ca2+-Independent Phospholipase A2 into a Phospholipid Bilayer via Coarse-Grained and Atomistic Molecular Dynamics Simulations.Bucher, D., Y.H. Hsu, V.D. Mouchlis, E.A. Dennis, J.A. McCammonPLoS Comp. Biol. 9, article e1003156 (2013). [PubMed PMC3723492]Group VI Ca
2+-independent phospholipase A
2 (iPLA
2) is a water-soluble enzyme that is active when associated with phospholipid membranes. Despite its clear pharmaceutical relevance, no X-ray or NMR structural information is currently available for the iPLA
2 or its membrane complex. In this paper, we combine homology modeling with coarse-grained (CG) and all-atom (AA) molecular dynamics (MD) simulations to build structural models of iPLA
2 in association with a phospholipid bilayer. CG-MD simulations of the membrane insertion process were employed to provide a starting point for an atomistic description. Six AA-MD simulations were then conducted for 60 ns, starting from different initial CG structures, to refine the membrane complex. The resulting structures are shown to be consistent with each other and with deuterium exchange mass spectrometry (DXMS) experiments, suggesting that our approach is suitable for the modeling of iPLA
2 at the membrane surface. The models show that an anchoring region (residues 710–724) forms an amphipathic helix that is stabilized by the membrane. In future studies, the proposed iPLA
2 models should provide a structural basis for understanding the mechanisms of lipid extraction and drug-inhibition. In addition, the dual-resolution approach discussed here should provide the means for the future exploration of the impact of lipid diversity and sequence mutations on the activity of iPLA
2 and related enzymes.
Influence of neighboring reactive particles on diffusion-limited reactions.Changsun Eun, Peter M. Kekenes-Huskey, and J. Andrew McCammonJ. Chem. Phys. 139, article number 044117 (2013). [PubMed PMC3745503]Competition between reactive species is commonplace in typical chemical reactions. Specifically the primary reaction between a substrate and its target enzyme may be altered when interactions with secondary species in the system are substantial. We explore this competition phenomenon for diffusion-limited reactions in the presence of neighboring particles through numerical solution of the diffusion equation. As a general model for globular proteins and small molecules, we consider spherical representations of the reactants and neighboring particles; these neighbors vary in local density, size, distribution, and relative distance from the primary target reaction, as well as their surface reactivity. Modulations of these model variables permit inquiry into the influence of excluded volume and competition on the primary reaction due to the presence of neighboring particles. We find that the surface reactivity effect is long-ranged and a strong determinant of reaction kinetics, whereas the excluded volume effect is relatively short-ranged and less influential in comparison. As a consequence, the effect of the excluded volume is only modestly dependent on the neighbor distribution and is approximately additive; this additivity permits a linear approximation to the many-body effect on the reaction kinetics. In contrast, the surface reactivity effect is non-additive, and thus it may require higher-order approximations to describe the reaction kinetics. Our model study has broad implications in the general understanding of competition and local crowding on diffusion-limited chemical reactions.
Multi-core CPU or GPU-accelerated Multiscale Modeling for Biomolecular Complexes.Liao, T., Y. Zhang, P. Kekenes-Huskey, Y. Cheng, A. Michailova, A.D. McCulloch, M. Holst, J.A. McCammonMolecular Based Mathematical Biology 1, 164-179 (2013). [PubMed PMC3858848]Multi-scale modeling plays an important role in understanding the structure and biological functionalities of large biomolecular complexes. In this paper, we present an efficient computational framework to construct multi-scale models from atomic resolution data in the Protein Data Bank (PDB), which is accelerated by multi-core CPU and programmable Graphics Processing Units (GPU). A multi-level summation of Gaus-sian kernel functions is employed to generate implicit models for biomolecules. The coefficients in the summation are designed as functions of the structure indices, which specify the structures at a certain level and enable a local resolution control on the biomolecular surface. A method called neighboring search is adopted to locate the grid points close to the expected biomolecular surface, and reduce the number of grids to be analyzed. For a specific grid point, a KD-tree or bounding volume hierarchy is applied to search for the atoms contributing to its density computation, and faraway atoms are ignored due to the decay of Gaussian kernel functions. In addition to density map construction, three modes are also employed and compared during mesh generation and quality improvement to generate high quality tetrahedral meshes: CPU sequential, multi-core CPU parallel and GPU parallel. We have applied our algorithm to several large proteins and obtained good results.
Inactivating mutation in Histone deacetylase 3 stabilizes its active conformation.Arrar, M., C.A.F. de Oliveira, J.A. McCammonProtein Science 22,1306-1312 (2013) [PubMed 23904210]Histone deacetylases (HDACs), together with histone acetyltransferases (HATs), regulate gene expression by modulating the acetylation level of chromatin. HDAC3 is implicated in many important cellular processes, particularly in cancer cell proliferation and metastasis, making inhibition of HDAC3 a promising epigenetic treatment for certain cancers. HDAC3 is activated upon complex formation with both inositol tetraphosphate (IP4) and the deacetylase-activating domain (DAD) of multi-protein nuclear receptor corepressor complexes. In previous studies, we have shown that binding of DAD and IP4 to HDAC3 significantly restricts its conformational space towards its stable ternary complex conformation, and suggest this to be the active conformation. Here, we report a single mutation of HDAC3 that is capable of mimicking the stabilizing effects of DAD and IP4, without the presence of either. This mutation, however, results in a total loss of deacetylase activity, prompting a closer evaluation of our understanding of the activation of HDAC3.
Simulations of Biased Agonists in the β2 Adrenergic Receptor with Accelerated Molecular Dynamics.Tikhonova, I., B. Selvam, A. Ivetac, J. Wereszczynski, J.A. McCammonBiochemistry 52 (33), pp 5593-5603 (2013) [PubMed 23879802]The biased agonism of the G protein-coupled receptors (GPCRs), where in addition to a traditional G protein-signaling pathway a GPCR promotes intracellular signals though β-arrestin, is a novel paradigm in pharmacology. Biochemical and biophysical studies have suggested that aGPCR forms a distinct ensemble of conformations signaling through the G protein and β-arrestin. Here we report on the dynamics of the β2 adrenergic receptor bound to the β-arrestin and G protein-biased agonists and the empty receptor to further characterize the receptor conformational changes caused by biased agonists. We use conventional and accelerated molecular dynamics (aMD) simulations to explore the conformational transitions of the GPCR from the active state to the inactive state. We found that aMD simulations enable monitoring of the transition within the nanosecond time scale while capturing the known microscopic characteristics of the inactive states, such as the ionic lock, the inward position of F6.44, and water clusters. Distinct conformational states are shown to be stabilized by each biased agonist. In particular, in simulations of the receptor with the β-arrestin-biased agonist N-cyclopentylbutanepherine, we observe a different pattern of motions in helix 7 when compared to simulations with the G protein-biased agonist salbutamol that involves perturbations of the network of interactions within the NPxxY motif. Understanding the network of interactions induced by biased ligands and the subsequent receptor conformational shifts will lead to development of more efficient drugs.
Variational Implicit-Solvent Modeling of Host-Guest Binding: A Case Study on Cucurbit[7]uril.Zhou, S., K.E. Rogers, C.A.F. de Oliveira, R. Baron, L.T. Cheng, J. Dzubiella, B. Li, J.A. McCammon.J. Chem. Theory Comput. 9(9): 4195–4204 (2013) [PubMed PMC3770055]The synthetic host cucurbit[7]uril (CB[7]) binds aromatic guests or metal complexes with ultrahigh affinity compared with that typically displayed in protein–ligand binding. Due to its small size, CB[7] serves as an ideal receptor–ligand system for developing computational methods for molecular recognition. Here, we apply the recently developed variational implicit-solvent model (VISM), numerically evaluated by the level-set method, to study hydration effects in the high-affinity binding of the B2 bicyclo[2.2.2]octane derivative to CB[7]. For the unbound host, we find that the host cavity favors the hydrated state over the dry state due to electrostatic effects. For the guest binding, we find reasonable agreement to experimental binding affinities. Dissection of the individual VISM free-energy contributions shows that the major driving forces are water-mediated hydrophobic interactions and the intrinsic (vacuum) host–guest van der Waals interactions. These findings are in line with recent experiments and molecular dynamics simulations with explicit solvent. It is expected that the level-set VISM, with further refinement on the electrostatic descriptions, can efficiently predict molecular binding and recognition in a wide range of future applications.
Substrate-dependent dynamics of UDP-galactopyranose mutase: implications for drug design.Boechi, L., C.A.F. de Oliveira, I. Da Fonseca, K. Kizjakina, P. Sobrado, J.J. Tanner, J.A. McCammonProtein Science Vol. 22, Issue 11, 14901-1501 (2013). [PubMed 23934860]Trypanosoma cruzi is the causative agent of Chagas disease, a neglected tropical disease that represents one of the major health challenges of the Latin American countries. Successful efforts were made during the last few decades to control the transmission of this disease, but there is still no treatment for the 10 million adults in the chronic phase of the disease. In
T. cruzi, as well as in other pathogens, the flavoenzyme UDP-galactopyranose mutase (UGM) catalyzes the conversion of UDP-galactopyranose to UDP-galactofuranose, a precursor of the cell surface β-galactofuranose that is involved in the virulence of the pathogen. The fact that UGM is not present in humans makes inhibition of this enzyme a good approach in the design of new Chagas therapeutics. By performing a series of computer simulations of
T. cruzi UGM in the presence or absence of an active site ligand, we address the molecular details of the mechanism that controls the uptake and retention of the substrate. The simulations suggest a modular mechanism in which each moiety of the substrate controls the flexibility of a different protein loop. Furthermore, the calculations indicate that interactions with the substrate diphosphate moiety are especially important for stabilizing the closed active site. This hypothesis is supported with kinetics measurements of site-directed mutants of
T. cruzi UGM. Our results extend our knowledge of UGM dynamics and offer new alternatives for the prospective design of drugs.
Improving the Efficiency of Free Energy Calculations in the Amber Molecular Dynamics Package.Kaus, J.W., L.T. Pierce, R.C. Walker, J.A. McCammon.J. Chem. Theory Comp. 9, 4131−4139 (2013).
Alchemical transformations are widely used methods to calculate free energies. Amber has
traditionally included support for alchemical transformations as part of the sander molecular
dynamics (MD) engine. Here, we describe the implementation of a more efficient approach to
alchemical transformations in the Amber MD package. Specifically, we have implemented this
new approach within the more computationally efficient and scalable pmemd MD engine that is
included with the Amber MD package. The majority of the gain in efficiency comes from the
improved design of the calculation, which includes better parallel scaling and reduction in
the calculation of redundant terms. This new implementation is able to reproduce results
from equivalent simulations run with the existing functionality but at 2.5 times greater
computational efficiency. This new implementation is also able to run softcore simulations
at the λ end states making direct calculation of free energies more accurate, compared to
the extrapolation required in the existing implementation. The updated alchemical
transformation functionality is planned to be included in the next major release of Amber
(scheduled for release in Q1 2014), available at http://ambermd.org, under the Amber license.
Mapping of Allosteric Druggable Sites in Activation-Associated Conformers of the M2 Muscarinic Receptor.Miao, Y., S.E. Nichols, J.A. McCammon.Chem. Biol. Drug Design 83, 237-246 (2014). [PubMed PMC2918726]G-protein coupled receptors (GPCRs) are key cellular signaling proteins and have
been targeted by ~30-40% of marketed drugs for treating many human diseases
including cancer and heart failure. Recently, we directly observed activation of
the M2 muscarinic receptor through long-timescale accelerated molecular dynamics
(aMD) simulation, which revealed distinct inactive, intermediate and active
conformers of the receptor. Here, FTMAP is applied to search for “hot spots” in
these activation-associated conformers using a library of 16 organic probe
molecules that represent fragments of potential drugs. Seven allosteric
(non-orthosteric) binding sites are identified in the M2 receptor through the
FTMAP analysis. These sites are distributed in the solvent-exposed extracellular
and intracellular mouth regions, as well as the lipid-exposed pockets formed by
the transmembrane α helices TM3-TM4, TM5-TM6 and TM7-TM1/TM2. They serve as
promising target sites for designing novel allosteric modulators as
receptor-selective drugs.
Molecular and Subcellular-Scale Modeling of Nucleotide Diffusion in the Cardiac Myofilament LatticePeter M. Kekenes-Huskey, Tao Liao, Andrew K. Gillette, Johan E. Hake, Yongjie Zhang, Anushka P. Michailova, Andrew D. McCulloch, and J. Andrew McCammonBiophys. J 105,2130-2140 (2013) [PubMed PMC3835335]Contractile function of cardiac cells is driven by the sliding displacement of myofilaments powered by the cycling myosin crossbridges. Critical to this process is the availability of ATP, which myosin hydrolyzes during the cross-bridge cycle. The diffusion of adenine nucleotides through the myofilament lattice has been shown to be anisotropic, with slower radial diffu- sion perpendicular to the filament axis relative to parallel, and is attributed to the periodic hexagonal arrangement of the thin (actin) and thick (myosin) filaments. We investigated whether atomistic-resolution details of myofilament proteins can refine coarse-grain estimates of diffusional anisotropy for adenine nucleotides in the cardiac myofibril, using homogenization theory and atomistic thin filament models from the Protein Data Bank. Our results demonstrate considerable anisotropy in ATP and ADP diffusion constants that is consistent with experimental measurements and dependent on lattice spacing and myofilament overlap. A reaction-diffusion model of the half-sarcomere further suggests that diffusional anisotropy may lead to modest adenine nucleotide gradients in the myoplasm under physiological conditions.
Structure, mechanism, and dynamics of UDP-galactopyranose mutase.Tanner, J.J., L Boechi, J.A. McCammon, P. Sobrado.Archives of Biochemistry and Biophysics, Vol. 544, 128-141 (2014) [PubMed 24096172]The flavoenzyme UDP-galactopyranose mutase (UGM) is a key enzyme in
galactofuranose biosynthesis. The enzyme catalyzes the 6-to-5 ring contraction
of UDP-galactopyranose to UDP-galactofuranose. Galactofuranose is absent in
humans yet is an essential component of bacterial and fungal cell walls and a
cell surface virulence factor in protozoan parasites. Thus, inhibition of
galactofuranose biosynthesis is a valid strategy for developing new
antimicrobials. UGM is an excellent target in this effort because the product of
the UGM reaction represents the first appearance of galactofuranose in the
biosynthetic pathway. The UGM reaction is redox neutral, which is atypical for
flavoenzymes, motivating intense examination of the chemical mechanism and
structural features that tune the flavin for its unique role in catalysis. These
studies show that the flavin functions as nucleophile, forming a flavin-sugar
adduct that facilitates galactose-ring opening and contraction. The
3-dimensional fold is novel and conserved among all UGMs, however the larger
eukaryotic enzymes have additional secondary structure elements that lead to
significant differences in quaternary structure, substrate conformation, and
conformational flexibility. Here we present a comprehensive review of UGM
three-dimensional structure, provide an update on recent developments in
understanding the mechanism of the enzyme, and summarize computational studies
of active site flexibility.
Accelerated Molecular Dynamics Simulations with the AMOEBA Polarizable Force Field on graphics processing units.Lindert, S., D. Bucher, P. Eastman, V. Pande, J.A. McCammon.J. Chem. Theory Comp. 9 (11), pp 4684-4691 (2013) [PubMed PMC3948463]The accelerated molecular dynamics (aMD) method has recently been shown to enhance the sampling of biomolecules in molecular dynamics (MD) simulations, often by several orders of magnitude. Here, we describe an implementation of the aMD method for the OpenMM application layer that takes full advantage of graphics processing units (GPUs) computing. The aMD method is shown to work in combination with the AMOEBA polarizable force field (AMOEBA-aMD), allowing the simulation of long timescale events with a polarizable force field. Benchmarks are provided to show that the AMOEBA-aMD method is efficiently implemented and produces accurate results in its standard parametrization. For the BPTI protein, we demonstrate that the protein structure described with AMOEBA remains stable even on the extended time scales accessed at high levels of accelerations. For the DNA repair metalloenzyme endonuclease IV, we show that the use of the AMOEBA force field is a significant improvement over fixed charged models for describing the enzyme active-site. The new AMOEBA-aMD method is publicly available (http://wiki.simtk. org/openmm/VirtualRepository) and promises to be interesting for studying complex systems that can benefit from both the use of a polarizable force field and enhanced sampling.
Utilizing a dynamical description of IspH to aid in the development of novel antimicrobial drugs.Blachly, P.G., C.A.F. de Oliveira, S.L. Williams, J.A. McCammon.PLoS Comp. Biol. v.9(12) (2013). [PubMed PMC3868525]The nonmevalonate pathway is responsible for isoprenoid production in microbes, including H. pylori, M. tuberculosis and P. falciparum, but is nonexistent in humans, thus providing a desirable route for antibacterial and antimalarial drug discovery. We coordinate a structural study of IspH, a [4Fe-4S] protein responsible for converting HMBPP to IPP and DMAPP in the ultimate step in the nonmevalonate pathway. By performing accelerated molecular dynamics simulations on both substrate-free and HMBPP-bound [Fe4S4]2+ IspH, we elucidate how substrate binding alters the dynamics of the protein. Using principal component analysis, we note that while substrate-free IspH samples various open and closed conformations, the closed conformation observed experimentally for HMBPP-bound IspH is inaccessible in the absence of HMBPP. In contrast, simulations with HMBPP bound are restricted from accessing the open states sampled by the substrate-free simulations. Further investigation of the substrate-free simulations reveals large fluctuations in the HMBPP binding pocket, as well as allosteric pocket openings – both of which are achieved through the hinge motions of the individual domains in IspH. Coupling these findings with solvent mapping and various structural analyses reveals alternative druggable sites that may be exploited in future drug design efforts.