Molecular Dynamics Free Energy Simulations Research Articles

Thermodynamic coupling between activation and inactivation gating in potassium channels revealed by free energy molecular dynamics simulations

The amount of ionic current flowing through K+ channels is determined by the interplay between two separate time-dependent processes: activation and inactivation gating. Activation is concerned with the stimulus-dependent opening of the main intracellular gate, whereas inactivation is a spontaneous conformational transition of the selectivity filter toward a nonconductive state occurring on a variety of timescales. A recent analysis of multiple x-ray structures of open and partially open KcsA channels revealed the mechanism by which movements of the inner activation gate, formed by the inner helices from the four subunits of the pore domain, bias the conformational changes at the selectivity filter toward a nonconductive inactivated state. This analysis highlighted the important role of Phe103, a residue located along the inner helix, near the hinge position associated with the opening of the intracellular gate. In the present study, we use free energy perturbation molecular dynamics simulations (FEP/MD) to quantitatively elucidate the thermodynamic basis for the coupling between the intracellular gate and the selectivity filter. The results of the FEP/MD calculations are in good agreement with experiments, and further analysis of the repulsive, van der Waals dispersive, and electrostatic free energy contributions reveals that the energetic basis underlying the absence of inactivation in the F103A mutation in KcsA is the absence of the unfavorable steric interaction occurring with the large Ile100 side chain in a neighboring subunit when the intracellular gate is open and the selectivity filter is in a conductive conformation. Macroscopic current analysis shows that the I100A mutant indeed relieves inactivation in KcsA, but to a lesser extent than the F103A mutant.

Free energy simulations of a GTPase: GTP and GDP binding to archaeal initiation factor 2.

Archaeal initiation factor 2 (aIF2) is a protein involved in the initiation of protein biosynthesis. In its GTP-bound, “ON” conformation, aIF2 binds an initiator tRNA and carries it to the ribosome. In its GDP-bound, “OFF” conformation, it dissociates from tRNA. To understand the specific binding of GTP and GDP and its dependence on the ON or OFF conformational state of aIF2, molecular dynamics free energy simulations (MDFE) are a tool of choice. However, the validity of the computed free energies depends on the simulation model, including the force field and the boundary conditions, and on the extent of conformational sampling in the simulations. aIF2 and other GTPases present specific difficulties; in particular, the nucleotide ligand coordinates a divalent Mg2+ ion, which can polarize the electronic distribution of its environment. Thus, a force field with an explicit treatment of electronic polarizability could be necessary, rather than a simpler, fixed charge force field. Here, we begin by comparing a fixed charge force field to quantum chemical calculations and experiment for Mg2+:phosphate binding in solution, with the force field giving large errors. Next, we consider GTP and GDP bound to aIF2 and we compare two fixed charge force fields to the recent, polarizable, AMOEBA force field, extended here in a simple, approximate manner to include GTP. We focus on a quantity that approximates the free energy to change GTP into GDP. Despite the errors seen for Mg2+:phosphate binding in solution, we observe a substantial cancellation of errors when we compare the free energy change in the protein to that in solution, or when we compare the protein ON and OFF states. Finally, we have used the fixed charge force field to perform MDFE simulations and alchemically transform GTP into GDP in the protein and in solution. With a total of about 200 ns of molecular dynamics, we obtain good convergence and a reasonable statistical uncertainty, comparable to the force field uncertainty, and somewhat lower than the predicted GTP/GDP binding free energy differences. The sign and magnitudes of the differences can thus be interpreted at a semiquantitative level, and are found to be consistent with the experimental binding preferences of ON- and OFF-aIF2.

An atomistic model for simulations of nilotinib and nilotinib/kinase binding

Nilotinib is a novel anticancer drug, which specifically binds to the Abl kinase and blocks its signaling activity. In order to model the nilotinib/protein interactions, we have developed a molecular mechanics force field for nilotinib, consistent with the CHARMM force field for proteins and nucleic acids. Atomic charges were derived by utilizing a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule that interacts with nilotinib fragments at six different positions. We investigated both neutral and protonated states of nilotinib. The final rms deviation between the ab initio and the force field energies, averaged over both forms, was equal 0.2 kcal/mol. The model reproduces the ab initio geometry and flexibility of nilotinib. To apply the force field to nilotinib/Abl simulations, it is also necessary to determine the most likely protein and nilotinib protonation state when it binds to Abl. This task was carried out using molecular dynamics free energy simulations. The simulations indicate that nilotinib can interact with Abl in protonated and deprotonated forms, with the protonated form more favoured for the interaction. In the course of our calculations, we established that the His361, a titratable amino acid residue that mediates the interaction, prefers to be neutral. These insights and models should be of interest for drug design.

The hidden energetics of ligand binding and activation in a glutamate receptor

Ionotropic glutamate receptors (iGluRs) are ligand-gated ion channels that mediate the majority of excitatory synaptic transmission in the central nervous system. The free energy of neurotransmitter-binding to the ligand-binding domains (LBDs) of iGluRs is converted into useful work to drive receptor activation. Here, the principal thermodynamic contributions from ligand-docking and ligand-induced LBD closure are computed for nine ligands of GluA2 using all-atom molecular dynamics free energy simulations. The results are validated by a comparison with experimentally measured apparent affinities to the isolated LBD. Features in the free energy landscapes governing LBD closure are critical determinants of binding free energies. An analysis of accessible LBD conformations transposed into the context of an intact GluA2 receptor reveals that the relative displacement of specific diagonal subunits in the tetrameric structure may be key to the action of partial agonists.

Molecular Mechanisms of Monovalent Cation Selectivity in Na+-Coupled Secondary Transporters from Free Energy Simulations with Implicit and Explicit Treatment of Electronic Effects

The ability to efficiently discriminate different cations is essential for the proper physiological functioning of channels and transporters. Strict monovalent selectivity of Na+-coupled secondary transporters is essential to realize the strong coupling between the transport of the substrate and Na+ flux. Loss of selectivity would result in the substrate being pumped outward. The x-ray structures of Na+-dependent amino-acid transporters (LeuT and Glt)[1,2], provides a great opportunity to better understand the molecular basis of monovalent cation selectivity.Extensive all-atom free energy molecular dynamics simulations with hamiltonian construted with additive, polarizable force-fields and QM/MM descriptions of the binding pocket are performed at various occupancy states of the binding sites. In this work, we demonstrate that there is a collective effect of multiple binding sites on a total selectivity for Na+ over K+ (and Tl+) in LeuT and Glt transporters. The polarizable models of metal ions and protein (for sites Na1 and Na2 identified in structural studies) employed to provide better understanding for ion selectivity [3] as well as to mechanism of ion-substrate coupling in Glt transporter. The role of local connectivity, site rigidity and atomic polarization in monovalent cation selectivity is discussed. We will also present results from decomposition analysis.[1] Boudker, O., Ryan, R.M., Yernool, D., Shimamoto, K., Gouaux, E., 2007, Nature 445: 387-393.[2] Yamashita, A., Singh, S.K., Kawate, T., Jin, Y., Gouaux, E., 2005, Nature 437: 215-223.[3] Lev B., Salohub D., Noskov S., Computational Life Sciences,Vol. 2, No. 1, pp. 12-20.

The Hidden Energetics of Ligand-Binding and Activation in a Glutamate Receptor

Ionotropic glutamate receptors (iGluRs) are ligand-gated ion channels that mediate the majority of excitatory synaptic transmission in the central nervous system. The free energy of neurotransmitter-binding to the ligand-binding domains (LBDs) of iGluRs is converted into useful work to induce conformational changes that cause receptor activation. Here, the energetic contributions are dissected into the “hidden” components associated with ligand-docking and ligand-induced LBD closure. The relevant thermodynamic contributions have been computed for nine ligands of GluA2 using all-atom molecular dynamics free energy simulations with explicit solvent. The results are compared with experimentally measured apparent affinities to the isolated LBD. An analysis of accessible LBD conformations transposed into the context of an intact GluA2 receptor reveals that the relative displacement of specific diagonal subunits in the tetrameric structure may play a key role in the action of partial agonists. Our X-ray crystal structure of a crosslink-stabilized LBD assembly shows that these diagonal subunits can be covalently linked without disrupting the fold of the LBDs.

On How the Conformational Cycle of the AcrB Efflux Pump is Coupled to Proton Translocation: A Theoretical Study Based on High-Resolution Structural Data

The AcrA/AcrB/TolC multidrug efflux pump confers Escherichia coli with antibiotic resistance by sequestering toxic compounds found within the periplasm and inner membrane and extruding them into the extracellular space. The AcrB trimer is the central component of this efflux complex; anchored in the inner membrane, it forms an asymmetric assembly that undergoes a conformational cycle in which each protomer adopts three different structures. As a result, substrates bound in the periplasmic domain of AcrB are projected into the TolC channel, which reaches beyond the outer membrane. Crucially, the conformational cycle within AcrB is driven by the translocation of protons down the gradient sustained by the inner membrane, through a mechanism that has not been characterized so far.Here, we investigate this microscopic mechanism through atomistic free-energy molecular dynamics simulations and electrostatic calculations, based upon novel high-resolution structural data for wild-type and mutagenized AcrB. Specifically, we assess the events associated with binding and release of protons within the membrane domain, and determine the mechanism by which these events are coupled to the reorganization of key transmembrane helices within each protomer. This investigation reveals how proton translocation influences both local and remote interactions within the protein, thereby modulating its structure.

Membrane/Toxin Interaction Energetics via Serial Multiscale Molecular Dynamics Simulations.

Computing free energies of complex biomolecular systems via atomistic (AT) molecular dynamics (MD) simulations remains a challenge due to the need for adequate sampling and convergence. Recent coarse-grained (CG) methodology allows simulations of significantly larger systems (∼10(6) to 10(8) atoms) over longer (μs/ms) time scales. Such CG models appear to be capable of making semiquantitative predictions. However, their ability to reproduce accurate thermodynamic quantities remains uncertain. We have recently used CG MD simulations to compute the potential of mean force (PMF) or free energy profile of a small peptide toxin interacting with a lipid bilayer along a 1D reaction coordinate. The toxin studied was VSTx1 (Voltage Sensor Toxin 1) from spider venom which inhibits the archeabacterial voltage-gated potassium (Kv) channel KvAP by binding to the voltage-sensor (VS) domains. Here, we re-estimate this PMF profile using (i) AT MD simulations with explicit membrane and solvent and (ii) an implicit membrane and solvent (generalized Born; GBIM) model where only the peptide was explicit. We used the CG MD free energy simulations to guide the setup of the corresponding AT MD simulations. The aim was to avoid local minima in the AT simulations which would be difficult over shorter AT time scales. A cross-comparison of the PMF profiles revealed a conserved topology, although there were differences in the magnitude of the free energies. The CG and AT simulations predicted a membrane/water interface free energy well of -27 and -23 kcal/mol, respectively (with respect to water). The GBIM model, however, gave a reduced interfacial free energy well (-12 kcal/mol). In addition, the CG and GBIM models predicted a free energy barrier of +61 and +96 kcal/mol, respectively, for positioning the toxin at the center of the bilayer, which was considerably smaller in the AT simulations (+26 kcal/mol). Thus, we present a framework for serially combining CG and AT simulations to estimate the free energy of peptide/membrane interactions. Such approaches for combining simulations at different levels of granularity will become increasingly important in future studies of complex membrane/protein systems.

A molecular mechanics model for imatinib and imatinib:kinase binding

Imatinib is an important anticancer drug, which binds specifically to the Abl kinase and blocks its signalling activity. To model imatinib:protein interactions, we have developed a molecular mechanics force field for imatinib and four close analogues, which is consistent with the CHARMM force field for proteins and nucleic acids. Atomic charges and Lennard-Jones parameters were derived from a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule interacting with imatinib fragments at 32 different positions. We considered both a neutral and a protonated imatinib. The final RMS deviation between the ab initio and force field energies, averaged over both forms, was 0.2 kcal/mol. The model also reproduces the ab initio geometry and flexibility of imatinib. To apply the force field to imatinib:Abl simulations, it is also necessary to determine the most likely imatinib protonation state when it binds to Abl. This was done using molecular dynamics free energy simulations, where imatinib is reversibly protonated during a series of MD simulations, both in solution and in complex with Abl. The simulations indicate that imatinib binds to Abl in its protonated, positively-charged form. To help test the force field and the protonation prediction, we did MD free energy simulations that compare the Abl binding affinities of two imatinib analogs, obtaining good agreement with experiment. Finally, two new imatinib variants were considered, one of which is predicted to have improved Abl binding. This variant could be of interest as a potential drug.

Tetracycline-Tet Repressor Binding Specificity: Insights from Experiments and Simulations

Tetracycline (Tc) antibiotics have been put to new uses in the construction of artificial gene regulation systems, where they bind to the Tet repressor protein (TetR) and modulate its affinity for DNA. Many Tc variants have been produced, both to overcome bacterial resistance and to achieve a broad range of binding strengths. To better understand TetR-Tc binding, we investigate a library of 16 tetracyclines, using fluorescence experiments and molecular dynamics free energy simulations (MDFE). The relative TetR binding free energies are computed by reversibly transforming one Tc variant into another during the simulation, with no adjustable parameters. The chemical variations involve polar and nonpolar substitutions along one entire edge of the elongated Tc structure, which provides many of the protein-ligand contacts. The binding constants span five orders of magnitude. The simulations reproduce the experimental binding free energies, when available, within the uncertainty of either method (±0.5 kcal/mol), and reveal many additional details. Contributions of individual Tc substituents are evaluated, along with their additivity and transferability among different positions on the Tc scaffold; differences between D- and B-class repressors are quantified. With increasing computer power, the MDFE approach provides an attractive complement to experiment and should play an increasing role in the understanding and engineering of protein-ligand recognition.

Elbow Flexibility of the kt38 RNA Kink-Turn Motif Investigated by Free-Energy Molecular Dynamics Simulations

Kink-turns (K-turns) are common structural motifs that can introduce sharp kinks into double-stranded RNA, and have been proposed to mediate large-scale motions in the ribosome. K-turns consist of a bulge loop region flanked by trans sugar-Hoogsteen G:A pairs, and the sharp kink conformation is stabilized by A-minor interactions (adenine contacting a G:C basepair in the minor groove). Umbrella-sampling molecular dynamics simulations were used to disrupt an A-minor interaction in the ribosomal kt38 turn and to calculate the associated free-energy change. Coupling of umbrella sampling with replica exchanges between neighboring umbrella-sampling intervals could further improve the convergence of the free-energy calculations. The simulations revealed a coupled A-minor disruption and global opening of the K-turn motif, and allowed us to characterize several intermediate A-minor conformations. The calculated free-energy profile indicated a meta-stable, semi-open structure of slightly higher free energy (∼1 kcal mol −1), separated by a small free-energy barrier (∼1.5 kcal mol −1) from the closed (highly kinked) form. Both K-turn states are stabilized by distinct variants of the A-minor interaction. Further opening of the K-turn structure required significantly larger free-energy changes. The semi-open form had a reduced kink angle compatible with experimental data on K-turn solution structures, and opening was coupled to a continuous global unwinding of the K-turn motif. The range of free-energy changes associated with kt38 opening and unwinding are compatible with the idea that K-turns may facilitate biologically relevant motions during large-scale ribosome dynamics.

Mechanism Of Ion-ion And Ion-substrate Coupling In Secondary Amino-acid Transporters

The x-ray structures of LeuT and Glt, bacterial homologues of Na+/Cl--dependent amino-acid transporters, provides a great opportunity to better understand the molecular basis of monovalent cation. Both proteins possess ion-binding sites selective for Na+ over K+ and Li+. Extensive QM/MM minimization combined with all-atom free energy molecular dynamics simulations of the LeuT and Glt transporters embedded in an explicit membrane are performed at different temperatures and various occupancy states of the binding sites to dissect the molecular mechanism of ion selectivity, coupling between co-transporter substrate and ions occupying binding sites. In this work, we demonstrate that there is a collective effect of multiple binding sites on a total selectivity for Na+ over Li+ both in LeuT and Glt. We also will discuss functional roles of different ion binding sites in the transport cycle. The role of local connectivity, site rigidity, atomic polarization and partial charge transfer in monovalent cation selectivity is discussed.

Small molecule hydration free energies in explicit solvent: An extensive test of fixed-charge atomistic simulations.

Using molecular dynamics free energy simulations with TIP3P explicit solvent, we compute the hydration free energies of 504 neutral small organic molecules and compare them to experiments. We find, first, good general agreement between the simulations and the experiments, with an RMS error of 1.24 kcal/mol over the whole set (i.e., about 2 kT) and a correlation coefficient of 0.89. Second, we use an automated procedure to identify systematic errors for some classes of compounds, and suggest some improvements to the force field. We find that alkyne hydration free energies are particularly poorly predicted due to problems with a Lennard-Jones well depth, and find that an alternate choice for this well depth largely rectifies the situation. Third, we study the non-polar component of hydration free energies - that is, the part that is not due to electrostatics. While we find that repulsive and attractive components of the non-polar part both scale roughly with surface area (or volume) of the solute, the total non-polar free energy does not scale with the solute surface area or volume, because it is a small difference between large components and is dominated by the deviations from the trend. While the methods used here are not new, this is a more extensive test than previous explicit solvent studies, and the size of the test set allows identification of systematic problems with force field parameters for particular classes of compounds. We believe that the computed free energies and components will be valuable to others in future development of force fields and solvation models.

Calculating Chemically Accurate Redox Potentials for Engineered Flavoproteins from Classical Molecular Dynamics Free Energy Simulations

The tricyclic isoalloxazine nucleus of the redox cofactors flavin mononucleotide (FMN) and flavin adenine dinucleotide (FAD) acts as an electron sink in life-sustaining biological electron transfer (eT). The functional diversity of flavin-containing proteins (flavoproteins) transcends that of free flavins. A large body of experimental evidence attributes natural control of flavoprotein-mediated eT to tuning of the thermodynamic driving force by the protein environment. Understanding and engineering such modulation by the protein environment of the flavin redox potential (DeltaE(o)) is valuable in biotechnology and device design. In this study we employed classical molecular dynamics free energy simulations (MDFES), within a thermodynamic integration (TI) formalism, to calculate the change in FMN first reduction potential (DeltaDeltaE(o)(ox/sq)) imparted by 6 flavoprotein active site mutations. The combined performance of the AMBER ff03 (protein) and GAFF (cofactor) force fields was benchmarked against experimental data for mutations close to the isoalloxazine re- and si-faces that perturb the wild-type DeltaE(o)(ox/sq) value in Anabaena flavodoxin. The classical alchemical approach used in this study overestimates the magnitude of DeltaE(o) values, in common with other studies. Nevertheless, chemically accurate DeltaDeltaE(o) values--calculated to within 1 kcal mol(-1) of the experimental value--were obtained for five of the six mutations studied. We have shown that this approach is practical for quantitative in silico screening of the effect of mutations on the first reduction potential where experimental values and structural data are available for the wild-type flavoprotein. This approach promises to be useful as an integral part of future interdisciplinary strategies to engineer desired thermodynamic properties in flavoproteins of biotechnological interest.

Control of Ion Selectivity in LeuT: Two Na+ Binding Sites with Two Different Mechanisms

The x-ray structure of LeuT, a bacterial homologue of Na+/Cl−-dependent neurotransmitter transporters, provides a great opportunity to better understand the molecular basis of monovalent cation selectivity in ion-coupled transporters. LeuT possesses two ion binding sites, NA1 and NA2, which are highly selective for Na+. Extensive all-atom free-energy molecular dynamics simulations of LeuT embedded in an explicit membrane are performed at different temperatures and various occupancy states of the binding sites to dissect the molecular mechanism of ion selectivity. The results show that the two binding sites display robust selectivity for Na+ over K+ or Li+, the competing ions of most similar radii. Of particular interest, the mechanism primarily responsible for selectivity for each of the two binding sites appears to be different. In NA1, selectivity for Na+ over K+ arises predominantly from the strong electrostatic field arising from the negatively charged carboxylate group of the leucine substrate coordinating the ion directly. In NA2, which comprises only neutral ligands, selectivity for Na+ is enforced by the local structural restraints arising from the hydrogen-bonding network and the covalent connectivity of the polypeptide chain surrounding the ion according to a “snug-fit” mechanism.

The Reaction Mechanism of Paraoxon Hydrolysis by Phosphotriesterase from Combined QM/MM Simulations

Molecular dynamics simulations employing combined quantum mechanical and molecular mechanical (QM/MM) potentials have been carried out to investigate the reaction mechanism of the hydrolysis of paraoxon by phosphotriesterase (PTE). We used a dual-level QM/MM approach that synthesizes accurate results from high-level electronic structure calculations with computational efficiency of semiempirical QM/MM potentials for free energy simulations. In particular, the intrinsic (gas-phase) energies of the active site in the QM region are determined by using density functional theory (B3LYP) and second-order Møller-Plesset perturbation theory (MP2) and the molecular dynamics free energy simulations are performed by using the mixed AM1:CHARMM potential. The simulation results suggest a revised mechanism for the phosphotriester hydrolysis mechanism by PTE. The reaction free energy profile is mirrored by structural motions of the binuclear metal center in the active site. The two zinc ions occupy a compact conformation with an average zinc-zinc distance of 3.5 +/- 0.1 A in the Michaelis complex, whereas it is elongated to 5.3 +/- 0.3 A at the transition state and product state. The substrate is loosely bound to the more exposed zinc ion (Znbeta2+) at an average distance of 3.8 A +/- 0.3 A. The P=O bond of the substrate paraoxon is activated by adopting a tight coordination to the Znbeta2+, releasing the coordinate to the bridging hydroxide ion and increasing its nucleophilicity. It was also found that a water molecule enters into the binding pocket of the loosely bound binuclear center, originally occupied by the nucleophilic hydroxide ion. We suggest that the proton of this water molecule is taken up by His254 at low pH or released to the solvent at high pH, resulting in a hydroxide ion that pulls the Znbeta2+ ion closer to form the compact configuration and restores the resting state of the enzyme.

Combined Simulation and Mutagenesis Analyses Reveal the Involvement of Key Residues for Peroxisome Proliferator-activated Receptorα Helix 12 Dynamic Behavior

The dynamic properties of helix 12 in the ligand binding domain of nuclear receptors are a major determinant of AF-2 domain activity. We investigated the molecular and structural basis of helix 12 mobility, as well as the involvement of individual residues with regard to peroxisome proliferator-activated receptor alpha (PPARalpha) constitutive and ligand-dependent transcriptional activity. Functional assays of the activity of PPARalpha helix 12 mutants were combined with free energy molecular dynamics simulations. The agreement between the results from these approaches allows us to make robust claims concerning the mechanisms that govern helix 12 functions. Our data support a model in which PPARalpha helix 12 transiently adopts a relatively stable active conformation even in the absence of a ligand. This conformation provides the interface for the recruitment of a coactivator and results in constitutive activity. The receptor agonists stabilize this conformation and increase PPARalpha transcription activation potential. Finally, we disclose important functions of residues in PPARalpha AF-2, which determine the positioning of helix 12 in the active conformation in the absence of a ligand. Substitution of these residues suppresses PPARalpha constitutive activity, without changing PPARalpha ligand-dependent activation potential.

Importance of Hydration and Dynamics on the Selectivity of the KcsA and NaK Channels

Fundamental concepts governing ion selectivity in narrow pores are reviewed and the microscopic factors responsible for the lack of selectivity of the NaK channel, which is structurally similar to the K+-selective KcsA channel, are elucidated on the basis of all-atom molecular dynamics free energy simulations. The results on NaK are contrasted and compared with previous studies of the KcsA channel. Analysis indicates that differences in hydration of the cation in the pore of NaK is at the origin of the lack of selectivity of NaK.

Structure-Selective Recognition by Voltammetry: Enantiomeric Determination of Amines Using Azophenolic Crowns in Aprotic Solvent

The enantiomeric recognition of amines by voltammetry using electroactive macrocyclic molecules, nitroazophenolic crown ethers, is reported. The oxidation potential of the nitroazophenol moiety in nitroazophenols with 18-crown-6 sensitively depends on the structure of alkyl amines. Based on this phenomenon, enantiomeric amines and even the quantitative assay of the R/S ratio in enantiomeric mixtures can be selectively recognized by using 18-crown-6 azophenol (3-H) with chiral centers. In the case of phenylglycinol, the association constants (K) of 3-H for the R and S forms have an R/S value of 3.5. The peak potential of the R form in square-wave voltammograms reproducibly differs from that of the S form by 32 mV, within which the peak potential linearly varies with the enantiomeric ratio. Free energy perturbation and molecular dynamics simulation provide deeper understanding of the enantiomeric recognition in this system. The theoretical analysis indicates that the free energy difference between diastereomeric complexes agrees well with the experimental results, and the pi-pi or charge-charge interaction plays a key role in enantiomeric recognition.

Binding Free Energy Contributions of Interfacial Waters in HIV-1 Protease/Inhibitor Complexes

Water molecules are commonly observed in crystal structures of protein-ligand complexes where they mediate protein-ligand binding. It is of considerable theoretical and practical importance to determine quantitatively the individual free energy contributions of these interfacial water molecules to protein-ligand binding and to elucidate factors that influence them. The double-decoupling free energy molecular dynamics simulation method has been used to calculate the binding free energy contribution for each of the four interfacial water molecules observed in the crystal structure of HIV-1 protease complexed with KNI-272, a potent inhibitor. While two of these water molecules contribute significantly to the binding free energy, the other two have close to zero contribution. It was further observed that the protonation states of two catalytic aspartate residues, Asp25 and Asp125, strongly influence the free energy contribution of a conserved water molecule Wat301 and that different inhibitors significantly influence the free energy contribution of Wat301. Our results have important implications on our understanding of the role of interfacial water molecules in protein-ligand binding and to structure-based drug design aimed at incorporating these interfacial water molecules into ligands.