Toward quantitative estimates of binding affinities for protein–ligand systems involving large inhibitor compounds: A steered molecular dynamics simulation route

Understanding binding mechanisms between enzymes and potential inhibitors and quantifying protein – ligand affinities in terms of binding free energy is of primary importance in drug design studies. In this respect, several approaches based on molecular dynamics simulations, often combined with docking techniques, have been exploited to investigate the physicochemical properties of complexes of pharmaceutical interest. Even if the geometric properties of a modeled protein – ligand complex can be well predicted by computational methods, it is still challenging to rank with chemical accuracy a series of ligand analogues in a consistent way. In this article, we face this issue calculating relative binding free energies of a focal adhesion kinase, an important target for the development of anticancer drugs, with pyrrolopyrimidine‐based ligands having different inhibitory power. To this aim, we employ steered molecular dynamics simulations combined with nonequilibrium work theorems for free energy calculations. This technique proves very powerful when a series of ligand analogues is considered, allowing one to tackle estimation of protein – ligand relative binding free energies in a reasonable time. In our cases, the calculated binding affinities are comparable with those recovered from experiments by exploiting the Michaelis – Menten mechanism with a competitive inhibitor.     

Fast and accurate predictions of binding free energies using MM‐PBSA and MM‐GBSA

In the drug discovery process, accurate methods of computing the affinity of small molecules with a biological target are strongly needed. This is particularly true for molecular docking and virtual screening methods, which use approximated scoring functions and struggle in estimating binding energies in correlation with experimental values. Among the various methods, MM‐PBSA and MM‐GBSA are emerging as useful and effective approaches. Although these methods are typically applied to large collections of equilibrated structures of protein‐ligand complexes sampled during molecular dynamics in water, the possibility to reliably estimate ligand affinity using a single energy‐minimized structure and implicit solvation models has not been explored in sufficient detail. Herein, we thoroughly investigate this hypothesis by comparing different methods for the generation of protein‐ligand complexes and diverse methods for free energy prediction for their ability to correlate with experimental values. The methods were tested on a series of structurally diverse inhibitors of Plasmodium falciparum DHFR with known binding mode and measured affinities. The results showed that correlations between MM‐PBSA or MM‐GBSA binding free energies with experimental affinities were in most cases excellent. Importantly, we found that correlations obtained with the use of a single protein‐ligand minimized structure and with implicit solvation models were similar to those obtained after averaging over multiple MD snapshots with explicit water molecules, with consequent save of computing time without loss of accuracy. When applied to a virtual screening experiment, such an approach proved to discriminate between true binders and decoy molecules and yielded significantly better enrichment curves. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Computer simulations of ligand–protein binding with ensembles of protein conformations: A Monte Carlo study of HIV‐1 protease binding energy landscapes

We present the results of molecular docking simulations with HIV‐1 protease for the sb203386 and skf107457 inhibitors by Monte Carlo simulated annealing. A simplified piecewise linear energy function, the standard AMBER force field, and the AMBER force field with solvation and a soft‐core smoothing component are employed in simulations with a single‐protein conformation to determine the relationship between docking simulations with a simple energy function and more realistic force fields. The temperature‐dependent binding free energy profiles of the inhibitors interacting with a single protein conformation provide a detailed picture of relative thermodynamic stability and a distribution of ligand binding modes in agreement with experimental crystallographic data. Using the simplified piecewise linear energy function, we also performed Monte Carlo docking simulations with an ensemble of protein conformations employing preferential biased sampling of low‐energy protein conformations, and the results are analyzed in connection with the free energy profiles. ©1999 John Wiley & Sons, Inc. Int J Quant Chem 72: 73–84, 1999     

Accurate and Rigorous Prediction of the Changes in Protein Free Energies in a Large‐Scale Mutation Scan

The prediction of mutation‐induced free‐energy changes in protein thermostability or protein–protein binding is of particular interest in the fields of protein design, biotechnology, and bioengineering. Herein, we achieve remarkable accuracy in a scan of 762 mutations estimating changes in protein thermostability based on the first principles of statistical mechanics. The remaining error in the free‐energy estimates appears to be due to three sources in approximately equal parts, namely sampling, force‐field inaccuracies, and experimental uncertainty. We propose a consensus force‐field approach, which, together with an increased sampling time, leads to a free‐energy prediction accuracy that matches those reached in experiments. This versatile approach enables accurate free‐energy estimates for diverse proteins, including the prediction of changes in the melting temperature of the membrane protein neurotensin receptor 1.     

Dielectric relaxation in proteins: a continuum electrostatics model incorporating dielectric heterogeneity of the protein and time‐dependent charges

A boundary element formulation of continuum electrostatics is used to examine time‐independent dielectric relaxation and screening in two proteins, and time‐dependent relaxation in two simpler solutes. Cytochrome c oxidation is modeled by inserting partial charges on the heme, using one to three dielectric regions in the protein. It was suggested recently that for charge insertion on a protein‐bound ligand, all or part of the ligand should be treated as a cavity within the protein medium. Here, the effect of an internal cavity surrounding the central heme atoms is examined, considering separately the static and relaxation (or reorganization) free energies. The former is the free energy to remove the redox electron while maintaining the rest of the structure and charge distribution fixed; the latter is the free energy associated with the relaxation into the product state after the corresponding constraints are released. The effect of the cavity is found to be small for the static free energy, while for the relaxation free energy it is large, as polarization of groups immediately around the heme dominates the relaxation. If the protein surface groups are treated as a distinct medium with a dielectric of 25 (as suggested by recent molecular dynamics simulations), the relaxation free energy decreases significantly (from −37.0 to −43.9 kcal/mol), compared to a model where the whole protein has a dielectric constant of two. Therefore, with this model, although polarization of groups immediately around the heme still dominates the relaxation, polar groups near the protein surface also contribute significantly, and solvent negligibly. The screening of an applied field within myoglobin is calculated, with the protein surrounded by either a low‐dielectric or a high‐dielectric glass. In the vicinity of the CO ligand, the screening is approximately isotropic with a low‐dielectric glass. It is anisotropic with a high‐dielectric glass, but the applied and local fields are still approximately parallel. This has implications for experiments that probe dielectric screening in proteins with the newly developed technique of vibrational Stark spectroscopy: with a high‐dielectric glass, a single, rotationally averaged screening factor can be used, the local field being about 1.65 times the applied field. Finally, we calculate the time‐dependent relaxation in response to instantaneous charge insertion within a spherical cavity in a Debye solvent, and to photoexcitation of a tryptophan solute, illustrating the extension of the boundary element formulation to time‐dependent problems. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 290–305, 2001     

WATsite: Hydration site prediction program with PyMOL interface

Water molecules that mediate protein–ligand interactions or are released from the binding site on ligand binding can contribute both enthalpically and entropically to the free energy of ligand binding. To elucidate the thermodynamic profile of individual water molecules and their potential contribution to ligand binding, a hydration site analysis program WATsite was developed together with an easy‐to‐use graphical user interface based on PyMOL. WATsite identifies hydration sites from a molecular dynamics simulation trajectory with explicit water molecules. The free energy profile of each hydration site is estimated by computing the enthalpy and entropy of the water molecule occupying a hydration site throughout the simulation. The results of the hydration site analysis can be displayed in PyMOL. A key feature of WATsite is that it is able to estimate the protein desolvation free energy for any user specified ligand. The WATsite program and its PyMOL plugin are available free of charge from http://people.pnhs.purdue.edu/~mlill/software . © 2014 Wiley Periodicals, Inc.     

Protein:Ligand binding free energies: A stringent test for computational protein design

A computational protein design method is extended to allow Monte Carlo simulations where two ligands are titrated into a protein binding pocket, yielding binding free energy differences. These provide a stringent test of the physical model, including the energy surface and sidechain rotamer definition. As a test, we consider tyrosyl‐tRNA synthetase (TyrRS), which has been extensively redesigned experimentally. We consider its specificity for its substrate l ‐tyrosine (l ‐Tyr), compared to the analogs d ‐Tyr, p‐acetyl‐, and p‐azido‐phenylalanine (ac‐Phe, az‐Phe). We simulate l ‐ and d ‐Tyr binding to TyrRS and six mutants, and compare the structures and binding free energies to a more rigorous “MD/GBSA” procedure: molecular dynamics with explicit solvent for structures and a Generalized Born + Surface Area model for binding free energies. Next, we consider l ‐Tyr, ac‐ and az‐Phe binding to six other TyrRS variants. The titration results are sensitive to the precise rotamer definition, which involves a short energy minimization for each sidechain pair to help relax bad contacts induced by the discrete rotamer set. However, when designed mutant structures are rescored with a standard GBSA energy model, results agree well with the more rigorous MD/GBSA. As a third test, we redesign three amino acid positions in the substrate coordination sphere, with either l ‐Tyr or d ‐Tyr as the ligand. For two, we obtain good agreement with experiment, recovering the wildtype residue when l ‐Tyr is the ligand and a d ‐Tyr specific mutant when d ‐Tyr is the ligand. For the third, we recover His with either ligand, instead of wildtype Gln. © 2015 Wiley Periodicals, Inc.     

A Bio‐inspired Approach for Chromophore Communication: Ligand‐to‐Ligand and Host‐to‐Guest Energy Transfer in Hybrid Crystalline Scaffolds

Efficient multiple‐chromophore coupling in a crystalline metal–organic scaffold was achieved by mimicking a protein system possessing 100 % energy‐transfer (ET) efficiency between a green fluorescent protein variant and cytochrome b₅₆₂. The two approaches developed for ET relied on the construction of coordination assemblies and host–guest coupling. Based on time‐resolved photoluminescence measurements in combination with calculations of the spectral overlap function and Förster radius, we demonstrated that both approaches resulted in a very high ET efficiency. In particular, the observed ligand‐to‐ligand ET efficiency value was the highest reported so far for two distinct ligands in a metal–organic framework. These studies provide important insights for the rational design of crystalline hybrid scaffolds consisting of a large ensemble of chromophore molecules with the capability of directional ET.     

Computational design of protein–ligand binding: Modifying the specificity of asparaginyl‐tRNA synthetase

A method for computational design of protein–ligand interactions is implemented and tested on the asparaginyl‐ and aspartyl‐tRNA synthetase enzymes (AsnRS, AspRS). The substrate specificity of these enzymes is crucial for the accurate translation of the genetic code. The method relies on a molecular mechanics energy function and a simple, continuum electrostatic, implicit solvent model. As test calculations, we first compute AspRS‐substrate binding free energy changes due to nine point mutations, for which experimental data are available; we also perform large‐scale redesign of the entire active site of each enzyme (40 amino acids) and compare to experimental sequences. We then apply the method to engineer an increased binding of aspartyl‐adenylate (AspAMP) into AsnRS. Mutants are obtained using several directed evolution protocols, where four or five amino acid positions in the active site are randomized. Promising mutants are subjected to molecular dynamics simulations; Poisson‐Boltzmann calculations provide an estimate of the corresponding, AspAMP, binding free energy changes, relative to the native AsnRS. Several of the mutants are predicted to have an inverted binding specificity, preferring to bind AspAMP rather than the natural substrate, AsnAMP. The computed binding affinities are significantly weaker than the native, AsnRS:AsnAMP affinity, and in most cases, the active site structure is significantly changed, compared to the native complex. This almost certainly precludes catalytic activity. One of the designed sequences has a higher affinity and more native‐like structure and may represent a valid candidate for Asp activity. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Connectivity and binding‐site recognition: Applications relevant to drug design

Here, we describe a family of methods based on residue–residue connectivity for characterizing binding sites and apply variants of the method to various types of protein–ligand complexes including proteases, allosteric‐binding sites, correctly and incorrectly docked poses, and inhibitors of protein–protein interactions. Residues within ligand‐binding sites have about 25% more contact neighbors than surface residues in general; high‐connectivity residues are found in contact with the ligand in 84% of all complexes studied. In addition, a k‐means algorithm was developed that may be useful for identifying potential binding sites with no obvious geometric or connectivity features. The analysis was primarily carried out on 61 protein–ligand structures from the MEROPS protease database, 250 protein–ligand structures from the PDBSelect (25%), and 30 protein–protein complexes. Analysis of four proteases with crystal structures for multiple bound ligands has shown that residues with high connectivity tend to have less variable side‐chain conformation. The relevance to drug design is discussed in terms of identifying allosteric‐binding sites, distinguishing between alternative docked poses and designing protein interface inhibitors. Taken together, this data indicate that residue–residue connectivity is highly relevant to medicinal chemistry. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

An estimation method of binding free energy in terms of ABEEMσπ/MM and continuum electrostatics fused into LIE method

A method is proposed for the estimation of absolute binding free energy of interaction between proteins and ligands. The linear interaction energy method is combined with atom‐bond electronegativity equalization method at σπ level Force field (fused into molecular mechanics) and generalized Born continuum model calculation of electrostatic solvation for the estimation of the absolute free energy of binding. The parameters of this method are calibrated by using a training set of 24 HIV‐1 protease–inhibitor complexes (PDB entry 1AAQ). A correlation coefficient of 0.93 was obtained with a root mean square deviation of 0.70 kcal mol^?1. This approach is further tested on seven inhibitor and protease complexes, and it provides small root mean square deviation between the calculated binding free energy and experimental binding free energy without reparametrization. By comparing the radii of gyration and the hydrogen bond distances between ligand and protein of three training model molecules, the consistent comparison result of binding free energy is obtained. It proves that this method of calculating the binding free energy with appropriate structural analysis can be applied to quickly assess new inhibitors of HIV‐1 proteases. To test whether the parameters of this method can apply to other drug targets, we have validated this method for the drug target cyclooxygenase‐2. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Computational prediction of monosaccharide binding free energies to lectins with linear interaction energy models

The linear interaction energy (LIE) method to compute binding free energies is applied to lectin‐monosaccharide complexes. Here, we calculate the binding free energies of monosaccharides to the Ralstonia solanacearum lectin (RSL) and the Pseudomonas aeruginosa lectin‐II (PA‐IIL). The standard LIE model performs very well for RSL, whereas the PA‐IIL system, where ligand binding involves two calcium ions, presents a major challenge. To overcome this, we explore a new variant of the LIE model, where ligand–metal ion interactions are scaled separately. This model also predicts the saccharide binding preference of PA‐IIL on mutation of the receptor, which may be useful for protein engineering of lectins. © 2012 Wiley Periodicals, Inc.     

Site‐directed analysis on protein hydrophobicity

Hydrophobicity of a protein is considered to be one of the major intrinsic factors dictating the protein aggregation propensity. Understanding how protein hydrophobicity is determined is, therefore, of central importance in preventing protein aggregation diseases and in the biotechnological production of human therapeutics. Traditionally, protein hydrophobicity is estimated based on hydrophobicity scales determined for individual free amino acids, assuming that those scales are unaltered when amino acids are embedded in a protein. Here, we investigate how the hydrophobicity of constituent amino acid residues depends on the protein context. To this end, we analyze the hydration free energy—free energy change on hydration quantifying the hydrophobicity—of the wild‐type and 21 mutants of amyloid‐beta protein associated with Alzheimer's disease by performing molecular dynamics simulations and integral‐equation calculations. From detailed analysis of mutation effects on the protein hydrophobicity, we elucidate how the protein global factor such as the total charge as well as underlying protein conformations influence the hydrophobicity of amino acid residues. Our results provide a unique insight into the protein hydrophobicity for rationalizing and predicting the protein aggregation propensity on mutation, and open a new avenue to design aggregation‐resistant proteins as biotherapeutics. © 2014 Wiley Periodicals, Inc.     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Nonfitting protein–ligand interaction scoring function based on first‐principles theoretical chemistry methods: Development and application on kinase inhibitors

Targeted therapy is currently a hot topic in the fields of cancer research and drug design. An important requirement for this approach is the development of potent and selective inhibitors for the identified target protein. However, current ways to estimate inhibitor efficacy rely on empirical protein–ligand interaction scoring functions which, suffering from their heavy parameterizations, often lead to a low accuracy. In this work, we develop a nonfitting scoring function, which consists of three terms: (1) gas‐phase protein‐ligand binding enthalpy obtained by the eXtended ONIOM hybrid method based on an integration of density functional theory (DFT) methods (XYG3 and ωB97X‐D) and the semiempirical PM6 method, (2) solvation free energy based on DFT‐SMD solvation model, and (3) entropy effect estimated by using DFT frequency analysis. The new scoring function is tested on a cyclin‐dependent kinase 2 (CDK2) inhibitor database including 76 CDK2 protein inhibitors and a p21‐activated kinase 1 (PAK1) inhibitor database including 20 organometallic PAK1 protein inhibitors. From the results, good correlations are found between the calculated scores and the experimental inhibitor efficacies with the square of correlation coefficient R² of 0.76–0.88. This suggests a good predictive power of this scoring function. To the best of our knowledge, this is the first high level theory‐based nonfitting scoring function with such a good level of performance. This scoring function is recommended to be used in the final screening of lead structure derivatives. © 2013 Wiley Periodicals, Inc.     

New parameterization approaches of the LIE method to improve free energy calculations of PlmII‐Inhibitors complexes

The standard parameterization of the Linear Interaction Energy (LIE) method has been applied with quite good results to reproduce the experimental absolute binding free energies for several protein–ligand systems. However, we found that this parameterization failed to reproduce the experimental binding free energy of Plasmepsin II (PlmII) in complexes with inhibitors belonging to four dissimilar scaffolds. To overcome this fact, we developed three approaches of LIE, which combine systematic approaches to predict the inhibitor‐specific values of α, β, and γ parameters, to gauge their ability to calculate the absolute binding free energies for these PlmII‐Inhibitor complexes. Specifically: (i) we modified the linear relationship between the weighted nonpolar desolvation ratio (WNDR) and the α parameter, by introducing two models of the β parameter determined by the free energy perturbation (FEP) method in the absence of the constant term γ, and (ii) we developed a new parameterization model to investigate the linear correlation between WNDR and the correction term γ. Using these parameterizations, we were able to reproduce the experimental binding free energy from these systems with mean absolute errors lower than 1.5 kcal/mol. © 2010 Wiley Periodicals, Inc. J Comput Chem 2010     

Expanding molecular modeling and design tools to non‐natural sidechains

Protein–protein interactions encode the wiring diagram of cellular signaling pathways and their deregulations underlie a variety of diseases, such as cancer. Inhibiting protein–protein interactions with peptide derivatives is a promising way to develop new biological and therapeutic tools. Here, we develop a general framework to computationally handle hundreds of non‐natural amino acid sidechains and predict the effect of inserting them into peptides or proteins. We first generate all structural files (pdb and mol2), as well as parameters and topologies for standard molecular mechanics software (CHARMM and Gromacs). Accurate predictions of rotamer probabilities are provided using a novel combined knowledge and physics based strategy. Non‐natural sidechains are useful to increase peptide ligand binding affinity. Our results obtained on non‐natural mutants of a BCL9 peptide targeting beta‐catenin show very good correlation between predicted and experimental binding free‐energies, indicating that such predictions can be used to design new inhibitors. Data generated in this work, as well as PyMOL and UCSF Chimera plug‐ins for user‐friendly visualization of non‐natural sidechains, are all available at http://www.swisssidechain.ch . Our results enable researchers to rapidly and efficiently work with hundreds of non‐natural sidechains. © 2012 Wiley Periodicals, Inc.     

Including ligand‐induced protein flexibility into protein tunnel prediction

In proteins with buried active sites, understanding how ligands migrate through the tunnels that connect the exterior of the protein to the active site can shed light on substrate specificity and enzyme function. A growing body of evidence highlights the importance of protein flexibility in the binding site on ligand binding; however, the influence of protein flexibility throughout the body of the protein during ligand entry and egress is much less characterized. We have developed a novel tunnel prediction and evaluation method named IterTunnel, which includes the influence of ligand‐induced protein flexibility, guarantees ligand egress, and provides detailed free energy information as the ligand proceeds along the egress route. IterTunnel combines geometric tunnel prediction with steered molecular dynamics in an iterative process to identify tunnels that open as a result of ligand migration and calculates the potential of mean force of ligand egress through a given tunnel. Applying this new method to cytochrome P450 2B6, we demonstrate the influence of protein flexibility on the shape and accessibility of tunnels. More importantly, we demonstrate that the ligand itself, while traversing through a tunnel, can reshape tunnels due to its interaction with the protein. This process results in the exposure of new tunnels and the closure of preexisting tunnels as the ligand migrates from the active site. © 2014 Wiley Periodicals, Inc.     

Ligand‐Promoted Rhodium(III)‐Catalyzed ortho‐C−H Amination with Free Amines

Ligand development for rhodium(III)‐catalyzed C−H activation reactions has largely been limited to cyclopentadienyl (Cp) based scaffolds. 2‐Methylquinoline has now been identified as a feasible ligand that can coordinate to the metal center of Cp*RhCl to accelerate the cleavage of the C−H bond of N ‐pentafluorophenylbenzamides, providing a new structural lead for ligand design. The compatibility of this reaction with secondary free amines and anilines also overcomes the limitations of palladium(II)‐catalyzed C−H amination reactions.