Alchemical Free Energy Perturbation Research Articles

How MeCP2 and R.DpnI Proteins Recognize Methylated DNA

DNA methylation plays a major role in organismal development and the regulation of gene expression. Methylation of cytosine bases and its cellular roles in eukaryotes are well established, as well as methylation of adenine bases in bacterial genomes. Here, we present results from molecular dynamics simulations, alchemical free energy perturbation, and MM-PBSA calculations to explain the specificity of the R.DpnI enzyme for binding to adenine-methylated DNA in both its catalytic and winged-helix domains. We find that adenine-methylated DNA binds more favorably to the catalytic subunit of R.DpnI (−4 kcal/mol) and to the winged-helix domain (−1.6 kcal/mol) than unmethylated DNA.In particular, N6-adenine methylation is found to enthalpically stabilize binding to R.DpnI. In contrast, C5-cytosine methylation stabilizes binding to the MBD domain of the MeCP2 entropically with almost no difference in binding enthalpy.

Prediction of hydration free energies for the SAMPL4 diverse set of compounds using molecular dynamics simulations with the OPLS-AA force field

All-atom molecular dynamics computer simulations were used to blindly predict the hydration free energies of a range of small molecules as part of the SAMPL4 challenge. Compounds were parametrized on the basis of the OPLS-AA force field using three different protocols for deriving partial charges: (1) using existing OPLS-AA atom types and charges with minor adjustments of partial charges on equivalent connecting atoms and derivation of new parameters for a number of distinct chemical groups (N-alkyl imidazole, nitrate) that were not present in the published force field; (2) calculation of quantum mechanical charges via geometry optimization, followed by electrostatic potential (ESP) fitting, using Jaguar at the LMP2/cc-pVTZ(-F) level; and (3) via geometry optimization and CHelpG charges (Gaussian09at the HF/6-31G* level), followed by two-stage RESP fitting. The absolute hydration free energy was computed by an established protocol including alchemical free energy perturbation with thermodynamic integration. The use of standard OPLS-AA charges (protocol 1) with a number of newly parametrized charges and the use of histidine derived parameters for imidazole yielded an overall root mean square deviation of the prediction from the experimental data of 1.75kcal/mol. The precision of our results appears to be mainly limited by relatively poor reproducibility of the Lennard-Jones contribution towards the solvation free energy, for which we observed large variability that could be traced to a strong dependence on the initial system conditions.

Generalized Scalable Multiple Copy Algorithms for Biological Molecular Dynamics Simulations in NAMD

Computational methodologies that couple the dynamical evolution of a collection of copies, or replicas, of a system of interest offer powerful and flexible approaches to characterize complex molecular processes. Such multiple copy algorithms (MCAs) can be used to enhance sampling, compute reversible work and free energies, as well as refine transition pathways. Widely used examples of MCAs include temperature and Hamiltonian-tempering replica-exchange molecular dynamics (T-REMD and H-REMD), alchemical free energy perturbation with lambda replica exchange (FEP/-REMD), umbrella sampling with Hamiltonian replica exchange (US/H-REMD), and string method with swarms-of-trajectories to determine conformational transition pathways. Here, we report a robust and general implementation of MCAs for molecular dynamics (MD) simulations in the highly scalable program NAMD built upon the parallel programming system charm++. Multiple concurrent NAMD instances are launched with internal partitions of charm++ and located continuously within a single communication world. Messages between NAMD instances are passed by low-level point-to-point communication functions, which are accessible through NAMD's Tcl scripting interface. The communication-enabled Tcl scripting provides a sustainable application interface for end users to realize generalized MCAs without modifying source code. Illustrative applications of MCAs with fine-grained inter-copy communication structure, including global lambda exchange in FEP/-REMD, window swapping US/H-REMD in multidimensional order parameter space, and string method with swarms-of-trajectories were carried out on IBM Blue Gene/Q to demonstrate the versatility and massive scalability of the present implementation.

Generalized scalable multiple copy algorithms for molecular dynamics simulations in NAMD

Computational methodologies that couple the dynamical evolution of a set of replicated copies of a system of interest offer powerful and flexible approaches to characterize complex molecular processes. Such multiple copy algorithms (MCAs) can be used to enhance sampling, compute reversible work and free energies, as well as refine transition pathways. Widely used examples of MCAs include temperature and Hamiltonian-tempering replica-exchange molecular dynamics (T-REMD and H-REMD), alchemical free energy perturbation with lambda replica-exchange (FEP/λ-REMD), umbrella sampling with Hamiltonian replica exchange (US/H-REMD), and string method with swarms-of-trajectories conformational transition pathways. Here, we report a robust and general implementation of MCAs for molecular dynamics (MD) simulations in the highly scalable program NAMD built upon the parallel programming system Charm++. Multiple concurrent NAMD instances are launched with internal partitions of Charm++ and located continuously within a single communication world. Messages between NAMD instances are passed by low-level point-to-point communication functions, which are accessible through NAMD’s Tcl scripting interface. The communication-enabled Tcl scripting provides a sustainable application interface for end users to realize generalized MCAs without modifying the source code. Illustrative applications of MCAs with fine-grained inter-copy communication structure, including global lambda exchange in FEP/λ-REMD, window swapping US/H-REMD in multidimensional order parameter space, and string method with swarms-of-trajectories were carried out on IBM Blue Gene/Q to demonstrate the versatility and massive scalability of the present implementation.

Analysis of the Interactions between GABA(A) Receptors and T3 using Electrophysiology and Molecular Dynamics Simulations

Disorders of the thyroid cause a multitude of neurological dysfunctions including depression, anxiety, and psychosis. Thyroid hormones have been primarily thought to act via genomic mechanisms throughout the organism; however, recently another mechanism has been proposed for the adult brain. We hypothesize that the thyroid hormone triiodothyronine (T3) acts directly on GABA(A) receptors via a mechanism similar to that of neurosteroids. Previous electrophysiology experiments on expressed recombinant GABA(A) receptors demonstrated an inhibition of GABA responses in the presence of T3, and, at much higher concentrations of T3 alone, a direct stimulation of receptor activity. Current experiments are being conducted to investigate stereoselective effects and interactions with other compounds known to bind directly to the GABA(A) receptor. using a GABA(A) receptor homology model based on the crystal structure of the related Glutamate-gated chloride channel, atomic resolution molecular dynamics simulations were conducted to analyze the possible interaction of T3 and the GABA(A) receptor. In simulations, T3 is stable in binding sites in the transmembrane domain of the GABA(A)receptor in a region that is associated with the activation by neurosteroids. Alchemical free energy perturbation calculations are underway to test the affinity of T3 for the binding site in the GABA(A) receptor. Our results provide strong evidence supporting earlier experimental findings indicating a role of T3 in regulating the activity of GABA(A) receptors in brain, while bringing additional insight into both the molecular binding mode and the mechanism of modulation.

Theoretical studies of HIV-1 reverse transcriptase inhibition

Computational methods for accurately calculating the binding affinity of a ligand for a protein play a pivotal role in rational drug design. We herein present a theoretical study of the binding of five different ligands to one of the proteins responsible for the human immunodeficiency virus type 1 (HIV-1) cycle replication; the HIV-1 reverse transcriptase (RT). Two types of approaches are used based on molecular dynamics (MD) simulations within hybrid QM/MM potentials: the alchemical free energy perturbation method, FEP, and the pathway method, in which the ligand is physically pulled away from the binding site, thus rendering a potential of mean force (PMF) for the binding process. Our comparative analysis stresses their advantages and disadvantages and, although the results are not in quantitative agreement, both methods are capable of distinguishing the most and the less potent inhibitors of HIV-1 RT activity on an RNase H site. The methods can then be used to select the proper scaffold to design new drugs. A deeper analysis of these inhibitors through molecular electrostatic potential (MEP) and calculation of the binding contribution of the individual residues shows that, in a rational design, apart from the strong interactions established with the two magnesium cations present in the RNase H site, it is important to take into account interactions with His539 and with those residues that are anchoring the metals; Asp443, Glu478, Asp498 and Asp549. The MEPs of the active site of the protein and the different ligands show a better complementarity in those inhibitors that present higher binding energies, but there are still possibilities of improving the favourable interactions and decreasing those that are repulsive in order to design compounds with higher inhibitory activity.

Applying efficient implicit nongeometric constraints in alchemical free energy simulations.

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λ(i) ≤ 1 and Σ(i = 1)(N) λ(i) = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form: λ(i) = e(c sinθ(i))/Σ(j = 1)(N) e(c sinθ(j)) that oscillates between λ(i) = 0 and λ(i) = 1 and has relatively steep transitions between these endpoints.

Probing Carbohydrate Product Expulsion from a Processive Cellulase with Multiple Absolute Binding Free Energy Methods

Understanding the enzymatic mechanism that cellulases employ to degrade cellulose is critical to efforts to efficiently utilize plant biomass as a sustainable energy resource. A key component of cellulase action on cellulose is product inhibition from monosaccharide and disaccharides in the product site of cellulase tunnel. The absolute binding free energy of cellobiose and glucose to the product site of the catalytic tunnel of the Family 7 cellobiohydrolase (Cel7A) of Trichoderma reesei (Hypocrea jecorina) was calculated using two different approaches: steered molecular dynamics (SMD) simulations and alchemical free energy perturbation molecular dynamics (FEP/MD) simulations. For the SMD approach, three methods based on Jarzynski's equality were used to construct the potential of mean force from multiple pulling trajectories. The calculated binding free energies, −14.4 kcal/mol using SMD and −11.2 kcal/mol using FEP/MD, are in good qualitative agreement. Analysis of the SMD pulling trajectories suggests that several protein residues (Arg-251, Asp-259, Asp-262, Trp-376, and Tyr-381) play key roles in cellobiose and glucose binding to the catalytic tunnel. Five mutations (R251A, D259A, D262A, W376A, and Y381A) were made computationally to measure the changes in free energy during the product expulsion process. The absolute binding free energies of cellobiose to the catalytic tunnel of these five mutants are −13.1, −6.0, −11.5, −7.5, and −8.8 kcal/mol, respectively. The results demonstrated that all of the mutants tested can lower the binding free energy of cellobiose, which provides potential applications in engineering the enzyme to accelerate the product expulsion process and improve the efficiency of biomass conversion.

MSCALE: A General Utility for Multiscale Modeling

The combination of theoretical models of macromolecules that exist at different spatial and temporal scales has become increasingly important for addressing complex biochemical problems. This work describes the extension of concurrent multiscale approaches, introduces a general framework for carrying out calculations, and describes its implementation into the CHARMM macromolecular modeling package. This functionality, termed MSCALE, generalizes both the additive and subtractive multiscale scheme (e.g. QM/MM ONIOM-type), and extends its support to classical force fields, coarse grained modeling (e.g. ENM, GNM, etc.), and a mixture of them all. The MSCALE scheme is completely parallelized with each subsystem running as an independent, but connected calculation. One of the most attractive features of MSCALE is the relative ease of implementation using the standard MPI communication protocol. This allows external access to the framework and facilitates the combination of functionality previously isolated in separate programs. This new facility is fully integrated with free energy perturbation methods, Hessian based methods, and the use of periodicity and symmetry, which allows the calculation of accurate pressures. We demonstrate the utility of this new technique with four examples; (1) subtractive QM/MM and QM/QM calculations; (2) multi-force field alchemical free energy perturbation; (3) integration with the SANDER module of AMBER and the TINKER package to gain access to potentials not available in CHARMM; and (4) mixed resolution (i.e. coarse grain / all-atom) normal mode analysis. The potential of this new tool is clearly established and in conclusion an interesting mathematical problem is highlighted and future improvements are proposed.

Mechanistic and Thermodynamic Insights into the Transport Cycle of Lactose Permease

The Major Facilitator Superfamily (MFS) is an important class of membrane transporters whose members are found in almost all types of organisms and are very diverse in terms of substrate transport. Lactose permease of E. coli (LacY), which transports various sugar molecules across the plasma membrane, is studied as a model for the MFS proteins. LacY undergoes widespread conformational changes between a periplasmic-open and a cytoplasmic-open state during the transport cycle. The focus of this work is to determine the periplasmic-open structure starting with the known cytoplasmic-open crystal structure with an emphasis on the underlying mechanism. A two-step hybrid simulation approach that involves self-guided Langevin dynamics (SGLD) simulations to enhance conformational sampling is used. Our model for the periplasmic-open structure agrees well with various experimental findings. The protonation of Glu269 and binding of a disaccharide to LacY is found to trigger the conformational transition to the periplasmic-open state. The residues and interhelical interactions that are important during the conformational transition are determined. This two-step simulation approach can be extended to determine the conformational changes in other membrane proteins. In addition, sugar binding during structural changes follows the well-established ‘alternating access model’. A more extensive study on sugar binding is carried out to understand the ‘anomeric binding phenomenon’ of LacY (different binding affinities for different anomers of a disaccharide). The alchemical free energy perturbation (FEP) method is used to quantify binding affinities. The calculated binding free energy of α-D-galactopyranoside (NPG) to LacY agrees well with the experimental value (Nie et al., JBC, 2006). The FEP calculations for αβ-(Galp)2 and ββ-(Galp)2 reveal a detailed description of LacY's anomeric binding phenomenon. This method is extended to different anomers of other disaccharides to elucidate anomeric binding.

Dominant Site Involved in Inhibition of Targets by Anesthetics Suggested Through Novel Free Energy Perturbation Calculations

It is well-established that general anesthetics bind directly to many of their targets, but the location of binding sites in a well-investigated target family, the cys-loop receptors, remains unresolved. Although the structure of eukaryotic cys-loop receptors has proven challenging to determine crystographically, a recent prokaryotic anesthetic-sensitive member of the family (GLIC) has been crystallized in complex with both a volatile and injected anesthetic. Detergent molecules occupy the pore of both the apo and liganded GLIC receptor, preventing any direct observation of pore-block by anesthetics from the crystal structure. In recent microsecond atomistic molecular dynamics (MD) simulations of flooding by isoflurane of both the eukaryotic nicotinic acetylcholine receptor (nAChR) and GLIC, we observed two isoflurane molecules bind to the pore of each receptor. Circumventing the experimental obstacles, computational measurements of binding affinity to the pore site provide an alternate method for determining the relevance of pore-block to GLIC inhibition by anesthetics. In addition to a calculated binding affinity for isoflurane in the pore that corresponds well to the micromolar concentrations at which inhibition is observed, we also present novel measurements using the alchemical free energy perturbation method (FEP) to quantify negative cooperativity between the two isoflurane molecules in the pore. The low Hill coefficient (l0.5) calculated in this manner corresponds to that measured experimentally within 10%. The results suggest that the pore is probably occupied in GLIC exposed to experimental concentrations of isoflurane, and pore-block may provide the dominant source of inhibition.

Evidence for Substrate Preorganization in the Peptidylglycine α-Amidating Monooxygenase Reaction Describing the Contribution of Ground State Structure to Hydrogen Tunneling

Peptidylglycine α-amidating monooxygenase (PAM) is a bifunctional enzyme which catalyzes the post-translational modification of inactive C-terminal glycine-extended peptide precursors to the corresponding bioactive α-amidated peptide hormone. This conversion involves two sequential reactions both of which are catalyzed by the separate catalytic domains of PAM. The first step, the copper-, ascorbate-, and O(2)-dependent stereospecific hydroxylation at the α-carbon of the C-terminal glycine, is catalyzed by peptidylglycine α-hydroxylating monooxygenase (PHM). The second step, the zinc-dependent dealkylation of the carbinolamide intermediate, is catalyzed by peptidylglycine amidoglycolate lyase. Quantum mechanical tunneling dominates PHM-dependent C(α)-H bond activation. This study probes the substrate structure dependence of this chemistry using a set of N-acylglycine substrates of varying hydrophobicity. Primary deuterium kinetic isotope effects (KIEs), molecular mechanical docking, alchemical free energy perturbation, and equilibrium molecular dynamics were used to study the role played by ground-state substrate structure on PHM catalysis. Our data show that all Ν-acylglycines bind sequentially to PHM in an equilibrium-ordered fashion. The primary deuterium KIE displays a linear decrease with respect to acyl chain length for straight-chain N-acylglycine substrates. Docking orientation of these substrates displayed increased dissociation energy proportional to hydrophobic pocket interaction. The decrease in KIE with hydrophobicity was attributed to a preorganization event which decreased reorganization energy by decreasing the conformational sampling associated with ground state substrate binding. This is the first example of preorganization in the family of noncoupled copper monooxygenases.

Combining ab initio quantum mechanics with a dipole-field model to describe acid dissociation reactions in water: First-principles free energy and entropy calculations

We introduce a novel approach to compute dissociation free energy and entropy values in simulations that employ a density functional theory description of the acidic moiety and of the solvent. The approach consists of utilizing an alchemical transformation of a weak acid A-COOH into the strong acid B-COOH, which makes it practical to employ alchemical free energy perturbation methods in the context of ab initio molecular dynamics simulations. The present alchemical transformation circumvents the need to tackle changes in the total number of electrons and atoms by replacing the chemical residue responsible for the change in acidity with an easily tunable external effective potential. Our investigation demonstrates that (1) a simple but effective class of external potentials that control acidity changes in the acetic/trifluoroacetic acid series can be achieved by replacing the methyl and trifluoromethyl substituents by screened dipoles. Using this dipole-field/quantum-mechanics (DF/QM) approach one can predict gas-phase geometries, proton dissociation energies, total dipole moments, and water binding energies in good agreement with full-QM values. (2) The resulting alchemical perturbation calculations are stable and well converged and allow one to compute absolute pK(a) values whose accuracy is limited primarily by the exchange-correlation functional employed: H-COOH=2.5+/-0.6 (full-QM calculation), 3.7 (exp); F(3)C-COOH=0.4+/-0.6 (DF/QM calculation), 0.5 (exp); H(3)C-COOH=3.1+/-0.7 (DF/QM calculation), 4.7 (exp); 3) Our DF/QM model predicts that the difference in acidity between H-COOH and H(3)C-COOH is dominated by solvent entropy effects, in excellent agreement with experimental observations. The calculated difference between the dissociation energies of these acids is DeltaDelta(d)U=0.0+/-0.26 kcal/mol while the experimental value is 0.0+/-0.1 kcal/mol.

Atomistic Model For the Outward-Facing State of Lactose Permease and Quantification of Ligand Binding

The Major Facilitator Superfamily (MFS) is an important class of membrane transporters whose members are found in unicellular organism to complex systems like the human central nervous system. Lactose Permease (LacY) of E. coli, an MFS protein that transports various sugar molecules across plasma membrane, has been studied as a prototype of secondary transport proteins. Though the structure of LacY in the cytoplasmic open state (inward-facing) has been determined, the atomic-level details of the periplasmic open structure (outward-facing) are unknown. Using a two-step hybrid simulation approach that involves Self-Guided Langevin Dynamics (SGLD) simulations with an implicit membrane followed by molecular dynamics simulations with an explicit membrane, we have obtained an opening of LacY on the periplasmic side that is consistent with the Double Electron-Electron Resonance (DEER) experiments (Smirnova et al., PNAS, 2007). The inward-facing state is stabilized by several helix-helix hydrogen bonds involving side chains of the residues N245, S41, E374, K42 and Q242 and mutations in these residues may destabilize these interactions and help crystallize the outward-facing state. We believe that this hybrid simulation approach can be extended to determine the unknown structures and mechanisms of other secondary active transporters. In addition, sugar binding to LacY has been investigated and our simulations support the alternating access model of sugar transport (same binding site accessible from either side of the membrane). More extensive studies on binding involved alchemical free energy perturbation calculations on αβ-(Galp)2 and ββ-(Galp)2 to better understand the anomeric binding effect. Our methods are tested by comparing the computed binding free energy values for NPG (p-nitrophenyl α-D-galactopyranoside) with the experimental values (Nie et al., JBC, 2006). This work will aid in understanding the effect of substrate structure and affinity to LacY.