Experimental Binding Free Energies Research Articles

Statistical Analysis on the Performance of Molecular Mechanics Poisson–Boltzmann Surface Area versus Absolute Binding Free Energy Calculations: Bromodomains as a Case Study

Binding free energy calculations that make use of alchemical pathways are becoming increasingly feasible thanks to advances in hardware and algorithms. Although relative binding free energy (RBFE) calculations are starting to find widespread use, absolute binding free energy (ABFE) calculations are still being explored mainly in academic settings due to the high computational requirements and still uncertain predictive value. However, in some drug design scenarios, RBFE calculations are not applicable and ABFE calculations could provide an alternative. Computationally cheaper end-point calculations in implicit solvent, such as molecular mechanics Poisson–Boltzmann surface area (MMPBSA) calculations, could too be used if one is primarily interested in a relative ranking of affinities. Here, we compare MMPBSA calculations to previously performed absolute alchemical free energy calculations in their ability to correlate with experimental binding free energies for three sets of bromodomain–inhibitor pairs. Different MMPBSA approaches have been considered, including a standard single-trajectory protocol, a protocol that includes a binding entropy estimate, and protocols that take into account the ligand hydration shell. Despite the improvements observed with the latter two MMPBSA approaches, ABFE calculations were found to be overall superior in obtaining correlation with experimental affinities for the test cases considered. A difference in weighted average Pearson () and Spearman () correlations of 0.25 and 0.31 was observed when using a standard single-trajectory MMPBSA setup ( = 0.64 and = 0.66 for ABFE; = 0.39 and = 0.35 for MMPBSA). The best performing MMPBSA protocols returned weighted average Pearson and Spearman correlations that were about 0.1 inferior to ABFE calculations: = 0.55 and = 0.56 when including an entropy estimate, and = 0.53 and = 0.55 when including explicit water molecules. Overall, the study suggests that ABFE calculations are indeed the more accurate approach, yet there is also value in MMPBSA calculations considering the lower compute requirements, and if agreement to experimental affinities in absolute terms is not of interest. Moreover, for the specific protein–ligand systems considered in this study, we find that including an explicit ligand hydration shell or a binding entropy estimate in the MMPBSA calculations resulted in significant performance improvements at a negligible computational cost.

Attach-Pull-Release Calculations of Ligand Binding and Conformational Changes on the First BRD4 Bromodomain.

Bromodomains, protein domains involved in epigenetic regulation, are able to bind small molecules with high affinity. In the present study, we report free energy calculations for the binding of seven ligands to the first BRD4 bromodomain, using the attach-pull-release (APR) method to compute the reversible work of removing the ligands from the binding site and then allowing the protein to relax conformationally. We test three different water models, TIP3P, TIP4PEw, and SPC/E, as well as the GAFF and GAFF2 parameter sets for the ligands. Our simulations show that the apo crystal structure of BRD4 is only metastable, with a structural transition happening in the absence of the ligand typically after 20 ns of simulation. We compute the free energy change for this transition with a separate APR calculation on the free protein and include its contribution to the ligand binding free energies, which generally causes an underestimation of the affinities. By testing different water models and ligand parameters, we are also able to assess their influence in our results and determine which one produces the best agreement with the experimental data. Both free energies associated with the conformational change and ligand binding are affected by the choice of water model, with the two sets of ligand parameters affecting their binding free energies to a lesser degree. Across all six combinations of water model and ligand potential function, the Pearson correlation coefficients between calculated and experimental binding free energies range from 0.55 to 0.83, and the root-mean-square errors range from 1.4–3.2 kcal/mol. The current protocol also yields encouraging preliminary results when used to assess the relative stability of ligand poses generated by docking or other methods, as illustrated for two different ligands. Our method takes advantage of the high performance provided by graphics processing units and can readily be applied to other ligands as well as other protein systems.

A Comprehensive Docking and MM/GBSA Rescoring Study of Ligand Recognition upon Binding Antithrombin.

Background: A high-throughput virtual screening pipeline has been extended from single energetically minimized structure Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) rescoring to ensemble-average MM/GBSA rescoring. The correlation coefficient (R2) of calculated and experimental binding free energies for a series of antithrombin ligands has been improved from 0.36 to 0.69 when switching from the single-structure MM/GBSA rescoring to ensemble-average one. The electrostatic interactions in both solute and solvent are identified to play an important role in determining the binding free energy after the decomposition of the calculated binding free energy. The increasing negative charge of the compounds provides a more favorable electrostatic energy change but creates a higher penalty for the solvation free energy. Such a penalty is compensated by the electrostatic energy change, which results in a better binding affinity. A highly hydrophobic ligand is determined by the docking program to be a non-specific binder.Results: Our results have demonstrated that it is very important to keep a few top poses for rescoring, if the binding is non-specific or the binding mode is not well determined by the docking calculation.

A Hybrid Knowledge-Based and Empirical Scoring Function for Protein–Ligand Interaction: SMoG2016

We present the third generation of our scoring function for the prediction of protein-ligand binding free energy. This function is now a hybrid between a knowledge-based potential and an empirical function. We constructed a diversified set of ∼1000 complexes from the PDBBinding-CN database for the training of the function, and we show that this number of complexes generates enough data to build the potential. The occurrence of 420 different types of atomic pairwise interactions is computed in up to five different ranges of distances to derive the knowledge-based part. All of the parameters were optimized, and we were able to considerably improve the accuracy of the scoring function with a Pearson correlation coefficient against experimental binding free energies of up to 0.57, which ranks our new scoring function as one of the best currently available and the second-best in terms of standard deviation (SD = 1.68 kcal/mol). The function was then further improved by inclusion of different terms taking into account repulsion and loss of entropy upon binding, and we show that it is capable of recovering native binding poses up to 80% of the time. All of the programs, tools, and protein sets are released in the Supporting Information or as open-source programs.

Fast and accurate determination of the relative binding affinities of small compounds to HIV-1 protease using non-equilibrium work.

The fast pulling ligand (FPL) out of binding cavity using non-equilibrium molecular dynamics (MD) simulations was demonstrated to be a rapid, accurate and low CPU demand method for the determination of the relative binding affinities of a large number of HIV-1 protease (PR) inhibitors. In this approach, the ligand is pulled out of the binding cavity of the protein using external harmonic forces, and the work of pulling force corresponds to the relative binding affinity of HIV-1 PR inhibitor. The correlation coefficient between the pulling work and the experimental binding free energy of R=-0.95 shows that FPL results are in good agreement with experiment. It is thus easier to rank the binding affinities of HIV-1 PR inhibitors, that have similar binding affinities because the mean error bar of pulling work amounts to δW=7%. The nature of binding is discovered using the FPL approach. © 2016 Wiley Periodicals, Inc.

An Atomistic Statistically Effective Energy Function for Computational Protein Design.

Shortcomings in the definition of effective free-energy surfaces of proteins are recognized to be a major contributory factor responsible for the low success rates of existing automated methods for computational protein design (CPD). The formulation of an atomistic statistically effective energy function (SEEF) suitable for a wide range of CPD applications and its derivation from structural data extracted from protein domains and protein-ligand complexes are described here. The proposed energy function comprises nonlocal atom-based and local residue-based SEEFs, which are coupled using a novel atom connectivity number factor to scale short-range, pairwise, nonbonded atomic interaction energies and a surface-area-dependent cavity energy term. This energy function was used to derive additional SEEFs describing the unfolded-state ensemble of any given residue sequence based on computed average energies for partially or fully solvent-exposed fragments in regions of irregular structure in native proteins. Relative thermal stabilities of 97 T4 bacteriophage lysozyme mutants were predicted from calculated energy differences for folded and unfolded states with an average unsigned error (AUE) of 0.84 kcal mol(-1) when compared to experiment. To demonstrate the utility of the energy function for CPD, further validation was carried out in tests of its capacity to recover cognate protein sequences and to discriminate native and near-native protein folds, loop conformers, and small-molecule ligand binding poses from non-native benchmark decoys. Experimental ligand binding free energies for a diverse set of 80 protein complexes could be predicted with an AUE of 2.4 kcal mol(-1) using an additional energy term to account for the loss in ligand configurational entropy upon binding. The atomistic SEEF is expected to improve the accuracy of residue-based coarse-grained SEEFs currently used in CPD and to extend the range of applications of extant atom-based protein statistical potentials.

A Three‐Site Mechanism for Agonist/Antagonist Selective Binding to Vasopressin Receptors

Molecular-dynamics simulations with metadynamics enhanced sampling reveal three distinct binding sites for arginine vasopressin (AVP) within its V2 -receptor (V2 R). Two of these, the vestibule and intermediate sites, block (antagonize) the receptor, and the third is the orthosteric activation (agonist) site. The contacts found for the orthosteric site satisfy all the requirements deduced from mutagenesis experiments. Metadynamics simulations for V2 R and its V1a R-analog give an excellent correlation with experimental binding free energies by assuming that the most stable binding site in the simulations corresponds to the experimental binding free energy in each case. The resulting three-site mechanism separates agonists from antagonists and explains subtype selectivity.

A Three‐Site Mechanism for Agonist/Antagonist Selective Binding to Vasopressin Receptors

AbstractMolecular‐dynamics simulations with metadynamics enhanced sampling reveal three distinct binding sites for arginine vasopressin (AVP) within its V2‐receptor (V2R). Two of these, the vestibule and intermediate sites, block (antagonize) the receptor, and the third is the orthosteric activation (agonist) site. The contacts found for the orthosteric site satisfy all the requirements deduced from mutagenesis experiments. Metadynamics simulations for V2R and its V1aR‐analog give an excellent correlation with experimental binding free energies by assuming that the most stable binding site in the simulations corresponds to the experimental binding free energy in each case. The resulting three‐site mechanism separates agonists from antagonists and explains subtype selectivity.

Structural and energetic comparison of the complexes of aminoglycosides with the model of the ribosomal A-site

Understanding the interactions which govern aminoglycoside antibiotic binding to ribosomal RNA is essential to propose modifications of these antibiotics. We have investigated the hydrogen bond patterns, solvent accessibility, and stacking energies of the nucleobases in twelve different aminoglycoside-RNA crystal complexes. These analyses pointed to some antibiotic-induced RNA structural differences that depend on the type of bound aminoglycoside. The largest differences were in the hydrogen bonding pattern in the vicinity of the U1406 and U1495 base pairs, especially in the complexes with geneticin and modified paromomycin. The complexes that stand out were the ones with neamine and kanamycin. We found that the solvent-accessible surface area buried upon aminoglycoside binding to RNA increases with the number of aminoglycoside rings but its correlation with the net charge of the antibiotic or experimental binding free energies was weak. We also investigated the dependence of other aminoglycoside characterizing descriptors, such as the number of rings and total charge, on the experimentally determined Gibbs energies. The correlation of ΔG with the total charge had the coefficient of determination R 2 over 0.8–0.9, depending on experimental data set, and was the highest of all descriptors.

Obtaining Binding Free Energy from a Path Sampling without Force Bias

Binding/dissociation free energy is among the most important quantities for computer-aided drug design. During generating dissociation pathways, the use of force bias can lead to the production of unnatural pathways, which may not sample more plausible dissociation configurations of protein/ligand complexes. Therefore, we apply the Parallel Cascade Selection Molecular Dynamics (PaCS-MD) [J. Chem. Phys. 139, 035103 (2013)] to generate multiple dissociation pathways without applying force biases and specific targeting scheme. We perform dissociation simulation for Hen-egg Lysozyme and a ligand Tri-N-Acetylchitotiose as a test case. Applying Umbrella Sampling (US) for the selected PaCS-MD snapshots, the obtained free energy profile is in good agreement with experimental binding free energy. In addition, from the trajectories of PaCS-MD, we build the Markov State Model and calculate the free energy profile. Without any further simulations as in US, the difference of free energy of the bound state and un-bound state is also well-reproduced.

Dynamics of Seeded Aβ40-Fibril Growth from Atomistic Molecular Dynamics Simulations: Kinetic Trapping and Reduced Water Mobility in the Locking Step.

Filamentous β-amyloid aggregates are crucial for the pathology of Alzheimer's disease. Despite the tremendous biomedical importance, the molecular pathway of growth propagation is not completely understood and remains challenging to investigate by simulations due to the long time scales involved. Here, we apply extensive all-atom molecular dynamics simulations in explicit water to obtain free energy profiles and kinetic information from position-dependent diffusion profiles for three different Aβ9-40-growth processes: fibril elongation by single monomers at the structurally unequal filament tips and association of larger filament fragments. Our approach provides insight into the molecular steps of the kinetic pathway and allows close agreement with experimental binding free energies and macroscopic growth rates. Water plays a decisive role, and solvent entropy is identified as the main driving force for assembly. Fibril growth is disfavored energetically due to cancellation of direct peptide-peptide interactions and solvation effects. The kinetics of growth is consistent with the characteristic dock/lock mechanism, and docking is at least 2 orders of magnitude faster. During initial docking, interactions are mediated by transient non-native hydrogen bonds, which efficiently catch the incoming monomer or fragment already at separations of about 3 nm. In subsequent locking, the dynamics is much slower due to formation of kinetically trapped conformations caused by long-lived non-native hydrogen bonds. Fibril growth additionally requires collective motion of water molecules to create a dry binding interface. Fibril growth is further retarded due to reduced mobility of the involved hydration water, evident from a 2-fold reduction of the diffusion coefficient.

Force Field for Mg(2+), Mn(2+), Zn(2+), and Cd(2+) Ions That Have Balanced Interactions with Nucleic Acids.

Divalent metal ions are of fundamental importance to the function and folding of nucleic acids. Divalent metal ion-nucleic acid interactions are complex in nature and include both territorial and site specific binding. Commonly employed nonbonded divalent ion models, however, are often parametrized against bulk ion properties and are subsequently utilized in biomolecular simulations without considering any data related to interactions at specific nucleic acid sites. Previously, we assessed the ability of 17 different nonbonded Mg(2+) ion models to reproduce different properties of Mg(2+) in aqueous solution including radial distribution functions, solvation free energies, water exchange rates, and translational diffusion coefficients. In the present work, we depart from the recently developed 12-6-4 potential models for divalent metal ions developed by Li and Merz and tune the pairwise parameters for Mg(2+), Mn(2+), Zn(2+), and Cd(2+) binding dimethyl phosphate, adenosine, and guanosine in order to reproduce experimental site specific binding free energies derived from potentiometric pH titration data. We further apply these parameters to investigate a metal ion migration previously proposed to occur during the catalytic reaction of the hammerhead ribozyme. The new parameters are shown to be accurate and balanced for nucleic acid binding in comparison with available experimental data and provide an important tool for molecular dynamics and free energy simulations of nucleic acids where these ions may exhibit different binding modes.

Examining the assumptions underlying continuum-solvent models.

Continuum-solvent models (CSMs) have successfully predicted many quantities, including the solvation-free energies (ΔG) of small molecules, but they have not consistently succeeded at reproducing experimental binding free energies (ΔΔG), especially for protein-protein complexes. Several CSMs break ΔG into the free energy (ΔGvdw) of inserting an uncharged molecule into solution and the free energy (ΔGel) gained from charging. Some further divide ΔGvdw into the free energy (ΔGrep) of inserting a nearly hard cavity into solution and the free energy (ΔGatt) gained from turning on dispersive interactions between the solute and solvent. We show that for 9 protein-protein complexes neither ΔGrep nor ΔGvdw was linear in the solvent-accessible area A, as assumed in many CSMs, and the corresponding components of ΔΔG were not linear in changes in A. We show that linear response theory (LRT) yielded good estimates of ΔGatt and ΔΔGatt, but estimates of ΔΔGatt obtained from either the initial or final configurations of the solvent were not consistent with those from LRT. The LRT estimates of ΔGel differed by more than 100 kcal/mol from the explicit solvent model's (ESM's) predictions, and its estimates of the corresponding component (ΔΔGel) of ΔΔG differed by more than 10 kcal/mol. Finally, the Poisson-Boltzmann equation produced estimates of ΔGel that were correlated with those from the ESM, but its estimates of ΔΔGel were much less so. These findings may help explain why many CSMs have not been consistently successful at predicting ΔΔG for many complexes, including protein-protein complexes.

Glycogen phosphorylase as a target for type 2 diabetes: synthetic, biochemical, structural and computational evaluation of novel N-acyl-N´-(β-D-glucopyranosyl) urea inhibitors.

Glycogen phosphorylase (GP), a validated target for the development of anti-hyperglycaemic agents, has been targeted for the design of novel glycopyranosylamine inhibitors. Exploiting the two most potent inhibitors from our previous study of N-acyl-β-D-glucopyranosylamines (Parmenopoulou et al., Bioorg. Med. Chem. 2014, 22, 4810), we have extended the linking group to -NHCONHCO- between the glucose moiety and the aliphatic/aromatic substituent in the GP catalytic site β-cavity. The N-acyl-N´-(β-D-glucopyranosyl) urea inhibitors were synthesized and their efficiency assessed by biochemical methods, revealing inhibition constant values of 4.95 µM and 2.53 µM. Crystal structures of GP in complex with these inhibitors were determined and analyzed, providing data for further structure based design efforts. A novel Linear Response - Molecular Mechanics Coulomb Surface Area (LR-MM-CBSA) method has been developed which relates predicted and experimental binding free energies for a training set of N-acyl-N´-(β-D-glucopyranosyl) urea ligands with a correlation coefficient R(2) of 0.89 and leave-one-out cross-validation (LOO-cv) Q(2) statistic of 0.79. The method has significant applications to direct future lead optimization studies, where ligand entropy loss on binding is revealed as a key factor to be considered. ADMET property predictions revealed that apart from potential permeability issues, the synthesized N-acyl-N´-(β-D-glucopyranosyl) urea inhibitors have drug-like potential without any toxicity warnings.

Structural and energetic insight into the interactions between the benzolactam inhibitors and tumor marker HSP90α

The heat shock protein 90α (HSP90α) provides a promising molecular target for cancer therapy. A series of novel benzolactam inhibitors exhibited distinct inhibitory activity for HSP90α. However, the structural basis for the impact of distinct R1 substituent groups of nine benzolactam inhibitors on HSP90α binding affinities remains unknown. In this study, we carried out molecular docking, molecular dynamics (MD) simulations, and molecular mechanics and generalized Born/surface area (MM–GBSA) binding free energy calculations to address the differences. Molecular docking studies indicated that all nine compounds presented one conformation in the ATP-binding site of HSP90α N-terminal domain. MD simulations and subsequent MM–GBSA calculations revealed that the hydrophobic interactions between all compounds and HSP90α contributed the most to the binding affinity and a good linear correlation was obtained between the calculated and the experimental binding free energies (R=0.88). The per residue decomposition revealed that the most remarkable differences of residue contributions were found in the residues Ala55, Ile96, and Leu107 defining a hydrophobic pocket for the R1 group, consistent with the analysis of binding modes. This study may be helpful for the future design of novel HSP90α inhibitors.

Dispersion- and Exchange-Corrected Density Functional Theory for Sodium Ion Hydration.

A challenge in density functional theory is developing functionals that simultaneously describe intermolecular electron correlation and electron delocalization. Recent exchange-correlation functionals address those two issues by adding corrections important at long ranges: an atom-centered pairwise dispersion term to account for correlation and a modified long-range component of the electron exchange term to correct for delocalization. Here we investigate how those corrections influence the accuracy of binding free energy predictions for sodium-water clusters. We find that the dual-corrected ωB97X-D functional gives cluster binding energies closest to high-level ab initio methods (CCSD(T)). Binding energy decomposition shows that the ωB97X-D functional predicts the smallest ion-water (pairwise) interaction energy and larger multibody contributions for a four-water cluster than most other functionals - a trend consistent with CCSD(T) results. Also, ωB97X-D produces the smallest amounts of charge transfer and the least polarizable waters of the density functionals studied, which mimics the lower polarizability of CCSD. When compared with experimental binding free energies, however, the exchange-corrected CAM-B3LYP functional performs best (error <1 kcal/mol), possibly because of its parametrization to experimental formation enthalpies. For clusters containing more than four waters, "split-shell" coordination must be considered to obtain accurate free energies in comparison with experiment.

Ametoctradin is a potent Qo site inhibitor of the mitochondrial respiration complex III.

Ametoctradin is a new Oomycete-specific fungicide under development by BASF. It is a potent inhibitor of the bc1 complex in mitochondrial respiration. However, its detailed action mechanism remains unknown. In the present work, the binding mode of ametoctradin was first uncovered by integrating molecular docking, MD simulations, and MM/PBSA calculations, which showed that ametoctradin should be a Q(o) site inhibitor of bc1 complex. Subsequently, a series of new 1,2,4-triazolo[1,5-a]pyrimidine derivatives were designed and synthesized to further understand the substituent effects on the 5- and 6-position of 1,2,4-triazolo[1,5-a]pyrimidine. The calculated binding free energies (ΔG(cal)) of newly synthesized analogues as Qo site inhibitors correlated very well (R(2) = 0.96) with their experimental binding free energies (ΔG(exp)). Two compounds (4a and 4c) with higher inhibitory activity against porcine SQR than ametoctradin were successfully identified. The structural and mechanistic insights obtained from the present study will provide a valuable clue for future designing of a new promising bc1 inhibitor.

Enantiomeric discrimination of chiral organic salts by chiral aza-15-crown-5 ether with C 1 symmetry: experimental and theoretical approaches.

The work involves an experimental ((1)H NMR) and theoretical (MD, MM-PBSA and DFT) investigation of the molecular recognition and discrimination properties of a chiral aza-15-crown-5 against methyl esters of alanine, phenylalanine and valine hydrochloride salts. The results indicate that the receptor binds enantiomers with moderate binding constants (88-1,389 M(-1)), with phenylalanine being more discriminated. The difference in experimental binding free energies (ΔG(R) - ΔG(S)) for alanine, phenylalanine and valine enantiomers were calculated as -0.36, -1.58 and 0.80 kcal mol(-1), respectively. The differences in theoretical binding energies were calculated by MM-PBSA (ΔE(R)PB - ΔE(S)PB=) as -0.30, -1.45 and 0.88, by B3LYP/6-31+G(d) as -1.17, -0.84 and 0.74 and by M06-2X/6-31+G(d) as -1.40, -3.26 and 1.66 kcal mol(-1). The data obtained give valuable information regarding the molecular recognition mode of the organoammonium complexes of chiral aza-crown ether with C 1 symmetry, which may be relevant to biological systems.

Data-driven prediction and design of bZIP coiled-coil interactions.

Selective dimerization of the basic-region leucine-zipper (bZIP) transcription factors presents a vivid example of how a high degree of interaction specificity can be achieved within a family of structurally similar proteins. The coiled-coil motif that mediates homo- or hetero-dimerization of the bZIP proteins has been intensively studied, and a variety of methods have been proposed to predict these interactions from sequence data. In this work, we used a large quantitative set of 4,549 bZIP coiled-coil interactions to develop a predictive model that exploits knowledge of structurally conserved residue-residue interactions in the coiled-coil motif. Our model, which expresses interaction energies as a sum of interpretable residue-pair and triplet terms, achieves a correlation with experimental binding free energies of R = 0.68 and significantly out-performs other scoring functions. To use our model in protein design applications, we devised a strategy in which synthetic peptides are built by assembling 7-residue native-protein heptad modules into new combinations. An integer linear program was used to find the optimal combination of heptads to bind selectively to a target human bZIP coiled coil, but not to target paralogs. Using this approach, we designed peptides to interact with the bZIP domains from human JUN, XBP1, ATF4 and ATF5. Testing more than 132 candidate protein complexes using a fluorescence resonance energy transfer assay confirmed the formation of tight and selective heterodimers between the designed peptides and their targets. This approach can be used to make inhibitors of native proteins, or to develop novel peptides for applications in synthetic biology or nanotechnology.

Protocols utilizing constant pH molecular dynamics to compute pH-dependent binding free energies.

In protein–ligand binding, the electrostatic environments of the two binding partners may vary significantly in bound and unbound states, which may lead to protonation changes upon binding. In cases where ligand binding results in a net uptake or release of protons, the free energy of binding is pH-dependent. Nevertheless, conventional free energy calculations and molecular docking protocols typically do not rigorously account for changes in protonation that may occur upon ligand binding. To address these shortcomings, we present a simple methodology based on Wyman’s binding polynomial formalism to account for the pH dependence of binding free energies and demonstrate its use on cucurbit[7]uril (CB[7]) host–guest systems. Using constant pH molecular dynamics and a reference binding free energy that is taken either from experiment or from thermodynamic integration computations, the pH-dependent binding free energy is determined. This computational protocol accurately captures the large pKa shifts observed experimentally upon CB[7]:guest association and reproduces experimental binding free energies at different levels of pH. We show that incorrect assignment of fixed protonation states in free energy computations can give errors of >2 kcal/mol in these host–guest systems. Use of the methods presented here avoids such errors, thus suggesting their utility in computing proton-linked binding free energies for protein–ligand complexes.