Molecular insight into the interaction between IFABP and PA by using MM-PBSA and alanine scanning methods

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method combined with alanine-scanning mutagenesis is a very important tool for rational drug design. In this study, molecular dynamics (MD) and MM-PBSA were applied to calculate the binding free energy between the rat intestinal fatty acid binding protein (IFABP) and palmitic acid (PA) to gain insight to the interaction details. Equally spaced snapshots along the trajectory were chosen to perform the binding free energy calculation, which yields a result highly consistent with experimental value with a deviation of 0.4 kcal/mol. Computational alanine scanning was performed on the same set of snapshots by mutating the residues in IFABP to alanine and recomputing the DeltaDeltaG(binding). By postprocessing a single trajectory of the wild-type complex, the average unsigned error of our calculated DeltaDeltaG(binding) is below 1.5 kcal/mol for most of the alanine mutations of the noncharged residues (67% in total). To further investigate some particular mutants, three additional dynamical simulations of IFABP Arg126Ala, Arg106Ala, and Arg106Gln mutants were conducted. Recalculated binding free energies are well consistent with the experimental data. Moreover, the ambiguous role of Arg106 caused by the free energy change of the opposite sign when it is mutated to alanine and glutamine respectively is clarified both structurally and energetically. Typically, this can be attributed to the partial electrostatic compensation mainly from Arg56 and the obvious entropy gain in Arg106Ala mutant while not in Arg106Gln mutant. The presented structural model of IFABP-PA complex could be used to guide future studies.     

Computational alanine scanning of the 1:1 human growth hormone-receptor complex

The MM-PBSA (Molecular Mechanics-Poisson-Boltzmann surface area) method was applied to the human Growth Hormone (hGH) complexed with its receptor to assess both the validity and the limitations of the computational alanine scanning approach. A 400-ps dynamical trajectory of the fully solvated complex was simulated at 300 K in a 101 A x 81 A x 107 A water box using periodic boundary conditions. Long-range electrostatic interactions were treated with the particle mesh Ewald (PME) summation method. Equally spaced snapshots along the trajectory were chosen to compute the binding free energy using a continuum solvation model to calculate the electrostatic desolvation free energy and a solvent-accessible surface area approach to treat the nonpolar solvation free energy. Computational alanine scanning was performed on the same set of snapshots by mutating the residues in the structural epitope of the hormone and the receptor to alanine and recomputing the deltaGbinding. To further investigate a particular structure, a 200-ps dynamical trajectory of an R43A hormone-receptor complex was simulated. By postprocessing a single trajectory of the wild-type complex, the average unsigned error of our calculated deltadeltaGbinding is approximately1 kcal/mol for the alanine mutations of hydrophobic residues and polar/charged residues without buried salt bridges. When residues involved in buried salt bridges are mutated to alanine, it is demonstrated that a separate trajectory of the alanine mutant complex can lead to reasonable agreement with experimental results. Our approach can be extended to rapid screening of a variety of possible modifications to binding sites.     

Unraveling the importance of protein-protein interaction: application of a computational alanine-scanning mutagenesis to the study of the IgG1 streptococcal protein G (C2 fragment) complex

Alanine-scanning mutagenesis of protein-protein interfacial residues is a very important process for rational drug design. In this study, we have used the improved MM-PBSA approach that combining molecular mechanics and continuum solvent permits one to calculate the free energy differences through alanine mutation. To identify the binding determinants of the complex formed between the IgG1 (immunoglobulin-binding protein G) and protein G, we have extended the experimental alanine scanning mutagenesis study to both proteins of this complex and, therefore, to all interfacial residues of this binding complex. As a result, we present new residues that can be characterized as warm spots and, therefore, are important for complex formation. We have further increased the understanding of the functionality of this improved computational alanine-scanning mutagenesis approach testing its sensitivity to a protein-protein complex with an interface made up of residues mainly polar. In this study, we also have improved the method for the detection of an important amino acid residue that frequently constitutes a hot spot--tryptophan.     

Unravelling Hot Spots: a comprehensive computational mutagenesis study

As protein–protein interactions are critical for all biological functions, representing a large and important class of targets for human therapeutics, identification of protein–protein interaction sites and detection of specific amino acid residues that contribute to the specificity and strength of protein interactions is very important in the biochemistry field. Alanine scanning mutagenesis has allowed the discovery of energetically crucial determinants for protein association that have been defined as hot spots. Systematic experimental mutagenesis is very laborious and time-consuming to perform, and thus it is important to achieve an accurate, predictive computational methodology for alanine scanning mutagenesis, capable of reproducing the experimental mutagenesis values. Having as a basis the MM–PBSA approach first developed by Massova et al., we performed a complete study of the influence of the variation of different parameters, such as the internal dielectric constant, the solvent representation, and the number of trajectories, in the accuracy of the free energy binding differences. As a result, we present here a very simple and fast methodological approach that achieved an overall success rate of 82% in reproducing the experimental mutagenesis data.Electronic Supplementary Material  Supplementary material is available for this article at and is accessible for authorized users.     

Relationship between Hot Spot Residues and Ligand Binding Hot Spots in Protein-Protein Interfaces

In the context of protein-protein interactions, the term "hot spot" refers to a residue or cluster of residues that makes a major contribution to the binding free energy, as determined by alanine scanning mutagenesis. In contrast, in pharmaceutical research, a hot spot is a site on a target protein that has high propensity for ligand binding and hence is potentially important for drug discovery. Here we examine the relationship between these two hot spot concepts by comparing alanine scanning data for a set of 15 proteins with results from mapping the protein surfaces for sites that can bind fragment-sized small molecules. We find the two types of hot spots are largely complementary; the residues protruding into hot spot regions identified by computational mapping or experimental fragment screening are almost always themselves hot spot residues as defined by alanine scanning experiments. Conversely, a residue that is found by alanine scanning to contribute little to binding rarely interacts with hot spot regions on the partner protein identified by fragment mapping. In spite of the strong correlation between the two hot spot concepts, they fundamentally differ, however. In particular, while identification of a hot spot by alanine scanning establishes the potential to generate substantial interaction energy with a binding partner, there are additional topological requirements to be a hot spot for small molecule binding. Hence, only a minority of hot spots identified by alanine scanning represent sites that are potentially useful for small inhibitor binding, and it is this subset that is identified by experimental or computational fragment screening.     

Protein-Protein Interactions: Insight from Molecular Dynamics Simulations and Nanoparticle Tracking Analysis

Protein-protein interaction plays an essential role in almost all cellular processes and biological functions. Coupling molecular dynamics (MD) simulations and nanoparticle tracking analysis (NTA) assay offered a simple, rapid, and direct approach in monitoring the protein-protein binding process and predicting the binding affinity. Our case study of designed ankyrin repeats proteins (DARPins)—AnkGAG1D4 and the single point mutated AnkGAG1D4-Y56A for HIV-1 capsid protein (CA) were investigated. As reported, AnkGAG1D4 bound with CA for inhibitory activity; however, it lost its inhibitory strength when tyrosine at residue 56 AnkGAG1D4, the most key residue was replaced by alanine (AnkGAG1D4-Y56A). Through NTA, the binding of DARPins and CA was measured by monitoring the increment of the hydrodynamic radius of the AnkGAG1D4-gold conjugated nanoparticles (AnkGAG1D4-GNP) and AnkGAG1D4-Y56A-GNP upon interaction with CA in buffer solution. The size of the AnkGAG1D4-GNP increased when it interacted with CA but not AnkGAG1D4-Y56A-GNP. In addition, a much higher binding free energy (∆GB) of AnkGAG1D4-Y56A (−31 kcal/mol) obtained from MD further suggested affinity for CA completely reduced compared to AnkGAG1D4 (−60 kcal/mol). The possible mechanism of the protein-protein binding was explored in detail by decomposing the binding free energy for crucial residues identification and hydrogen bond analysis.     

Molecular dynamics simulation study on stabilities and reactivities of NADH cytochrome B5 reductase

Binding free energies between coenzyme (FAD and NADH) and the apoenzyme of NADH-cytochrome b5 reductase (b5R) were estimated by applying the continuum Poisson-Boltzmann (PB) model to structures sampled from molecular dynamics simulations in explicit water molecules. Important residues for the enzymatic catalysis were clarified using a computational alanine scanning method. The binding free energies calculated by applying an alanine scanning method can successfully reproduce the trends of the measured steady-state enzymatic activities kcatNADH/KmNADH. Significant decreases in the binding free energy are expected when one of the four residues Arg91, Lys110, Ser127, and Thr181 is mutated into Ala. According to the results of the molecular dynamics simulation, Thr181 is considered to be one of the key residues that helps NADH to approach the isoalloxazine in FAD. Finally, we have constructed very simplified model systems and carried out density functional theory calculations using B3LYP/LANL2DZ//ROHF(or RHF)/LANL2DZ level of theory in order to elucidate a realistic and feasible mechanism of the hydride-ion transfer from NADH to FAD affected by HEME(Fe3+) as an electron acceptor. Our calculated results suggest that the electron and/or hydride-ion transfer reaction from NADH to FAD can be accelerated in the presence of HEME(Fe3+).     

Hot spot occlusion from bulk water: a comprehensive study of the complex between the lysozyme HEL and the antibody FVD1.3

Alanine scanning of protein-protein interfaces has shown that there are some residues in the protein-protein interfaces, responsible for most of the binding free energy, which are called hot spots. Hot spots tend to exist in densely packed central clusters, and a hypothesis has been proposed that considers that inaccessibility to the solvent must be a necessary condition to define a residue as a binding hot spot. This O-ring hypothesis is mainly based on the analysis of the accessible surface area (ASA) of 23 static, crystallographic structures of protein complexes. It is known, however, that protein flexibility allows for temporary exposures of buried interfacial groups, and even though the ASA provides a general trend of the propensity for hydration, protein/solvent-specific interactions or hydrogen bonding cannot be considered here. Therefore, a microscopic level, atomistic picture of hot spot solvation is needed to support the O-ring hypothesis. In this study, we began by applying a computational alanine-scanning mutagenesis technique, which reproduces the experimental results and allows for decomposing the binding free energy difference in its different energetic factors. Subsequently, we calculated the radial distribution function and residence times of the water molecules near the hot/warm spots to study the importance of the water environment around those energetically important amino acid residues. This study shows that within a flexible, dynamic protein framework, the warm/hot spot residues are, indeed, kept sheltered from the bulk solvent during the whole simulation, which allows a better interacting microenvironment.     

Role of protein structure and the role of individual fingers in zinc finger protein–DNA recognition: a molecular dynamics simulation study and free energy calculations

Molecular dynamics and MM_GBSA energy calculations on various zinc finger proteins containing three and four fingers bound to their target DNA gave insights into the role of each finger in the DNA binding process as part of the protein structure. The wild type Zif 268 (PDB code: 1AAY) gave a ΔG value of ??76.1 (14) kcal/mol. Zinc fingers ZF1, ZF2 and ZF3 were mutated in one experiment and in another experiment one finger was cut and the rest of the protein was studied for binding. The ΔΔG values for the Zinc Finger protein with both ZF1 and ZF2 mutated was +?80 kcal/mol, while mutating only ZF1 the ΔΔG value was +?52 kcal/mol (relative to the wild type). Cutting ZF3 and studying the protein consisting only of ZF1 linked to ZF2 gave a ΔΔG value of +?68 kcal/mol. Upon cutting ZF1, the resulting ZF2 linked to ZF3 protein gave a ΔΔG value of +?41 kcal/mol. The above results shed light on the importance of each finger in the binding process, especially the role of ZF1 as the anchoring finger followed in importance by ZF2 and ZF3. The energy difference between the binding of the wild type protein Zif268 (1AAY) and that for individual finger binding to DNA according to the formula: ΔΔG_{linkers, otherstructuralfactors}?=?ΔG_zif268???(ΔG_F1+F2+F3) gave a value?=???44.5 kcal/mol. This stabilization can be attributed to the contribution of linkers and other structural factors in the intact protein in the DNA binding process. DNA binding energies of variant proteins of the wild type Zif268 which differ in their ZF1 amino acid sequence gave evidence of a good relationship between binding energy and recognition and specificity, this finding confirms the reported vital role of ZF1 in the ZF protein scanning and anchoring to the target DNA sequence. The role of hydrogen bonds in both specific and nonspecific amino acid-DNA contacts is discussed in relation to mutations. The binding energies of variant Zinc Finger proteins confirmed the role of ZF1 in the recognition, specificity and anchoring of the zinc finger protein to DNA.     

Hot spot computational identification: Application to the complex formed between the hen egg white lysozyme (HEL) and the antibody HyHEL‐10

The definition and comprehension of the hot spots in an interface is a subject of primary interest for a variety of fields, including structure‐based drug design. Therefore, to achieve an alanine mutagenesis computational approach that is at the same time accurate and predictive, capable of reproducing the experimental mutagenesis values is a major challenge in the computational biochemistry field. Antibody/protein antigen complexes provide one of the greatest models to study protein–protein recognition process because they have three fundamentally features: specificity, high complementary association and a small epitope restricted to the diminutive complementary determining regions (CDR) region, while the remainder of the antibody is largely invariant. Thus, we apply a computational mutational methodological approach to the study of the antigen–antibody complex formed between the hen egg white lysozyme (HEL) and the antibody HyHEL‐10. A critical evaluation that focuses essentially on the limitations and advantages between different computational methods for hot spot determination, as well as between experimental and computational methodological approaches, is presented. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2007     

Homology modeling and molecular dynamic simulation of UDP-N-acetylmuramoyl-l-alanine-d-glutamate ligase (MurD) from Mycobacterium tuberculosis H37Rv using in silico approach

The present study aimed to identify the prospective inhibitors of MurD, a cytoplasmic enzyme that catalyzes the addition of d-glutamate to the UDP-N-acetylmuramoyl-l-alanine nucleotide precursor in Mycobacterium tuberculosis (MTB), using virtual screening, docking studies, pharmacokinetic analysis, Molecular Dynamic (MD) simulation, and Molecular Mechanics Generalized Born and Surface Area (MM-GBSA) analyses. The three dimensional (3D) structure was determined based on the homology technique using a template from Streptococcus agalactiae. The modeled structure had three binding sites, namely; substrate binding site (Val18, Thr19, Asp39, Asp40, Gly75, Asn147, Gln171 and His192), the ATP binding site (Gly123, Lys124, Thr125, Thr126, Glu166, Asp283, and Arg314) and the glutamic acid binding site (Arg382, Ser463, and Tyr470). These residues mentioned above play a critical role in the catalytic activity of the enzyme, and their inhibition could serve as a stumbling block to the normal function of the enzyme. A total of 10,344 obtained from virtual screened of Zinc and PubChem databases. These compounds further screened for Lipinski rule of five, docking studies and pharmacokinetic analysis. Four compounds with good binding energies (ZINC11881196 = −10.33 kcal/mol, ZINC12247644 = −8.90 kcal/mol, ZINC14995379 =−8.42 kcal/mol, and PubChem6185 = −8.20 kcal/mol), better than the binding energies of the ATP (−2.31 kcal/mol) and the ligand with known IC₅₀, Aminothiazole (−7.11 kcal/mol) were selected for the MD simulation and MM-GBSA analyses. The result of the analyses showed that all the four ligands formed a stable complex and had the binding free energies better than the binding energy of ATP. Therefore, these ligands considered as suitable prospective inhibitors of the MurD after experimental validation.     

Engineering protein therapeutics: predictive performances of a structure-based virtual affinity maturation protocol

The implementation of a structure based virtual affinity maturation protocol and evaluation of its predictivity are presented. The in silico protocol is based on conformational sampling of the interface residues (using the Dead End Elimination/A* algorithm), followed by the estimation of the change of free energy of binding due to a point mutation, applying MM/PBSA calculations. Several implementations of the protocol have been evaluated for 173 mutations in 7 different protein complexes for which experimental data were available: the use of the Boltzamnn averaged predictor based on the free energy of binding (ΔΔG(*)) combined with the one based on its polar component only (ΔΔE(pol*)) led to the proposal of a subset of mutations out of which 45% would have successfully enhanced the binding. When focusing on those mutations that are less likely to be introduced by natural in vivo maturation methods (99 mutations with at least two base changes in the codon), the success rate is increased to 63%. In another evaluation, focusing on 56 alanine scanning mutations, the in silico protocol was able to detect 89% of the hot-spots.     

Efficient calculation of relative binding free energies by umbrella sampling perturbation

An important task of biomolecular simulation is the calculation of relative binding free energies upon chemical modification of partner molecules in a biomolecular complex. The potential of mean force (PMF) along a reaction coordinate for association or dissociation of the complex can be used to estimate binding affinities. A free energy perturbation approach, termed umbrella sampling (US) perturbation, has been designed that allows an efficient calculation of the change of the PMF upon modification of a binding partner based on the trajectories obtained for the wild type reference complex. The approach was tested on the interaction of modified water molecules in aqueous solution and applied to in silico alanine scanning of a peptide‐protein complex. For the water interaction test case, excellent agreement with an explicit PMF calculation for each modification was obtained as long as no long range electrostatic perturbations were considered. For the alanine scanning, the experimentally determined ranking and binding affinity changes upon alanine substitutions could be reproduced within 0.1–2.0 kcal/mol. In addition, good agreement with explicitly calculated PMFs was obtained mostly within the sampling uncertainty. The combined US and perturbation approach yields, under the condition of sufficiently small system modifications, rigorously derived changes in free energy and is applicable to any PMF calculation. © 2014 Wiley Periodicals, Inc.     

Detailed potential of mean force studies on host–guest systems from the SAMPL6 challenge

Accurately predicting receptor–ligand binding free energies is one of the holy grails of computational chemistry with many applications in chemistry and biology. Many successes have been reported, but issues relating to sampling and force field accuracy remain significant issues affecting our ability to reliably calculate binding free energies. In order to explore these issues in more detail we have examined a series of small host–guest complexes from the SAMPL6 blind challenge, namely octa-acids (OAs)–guest complexes and Curcurbit[8]uril (CB8)–guest complexes. Specifically, potential of mean force studies using umbrella sampling combined with the weighted histogram method were carried out on both systems with both known and unknown binding affinities. We find that using standard force fields and straightforward simulation protocols we are able to obtain satisfactory results, but that simply scaling our results allows us to significantly improve our predictive ability for the unknown test sets: the overall RMSD of the binding free energy versus experiment is reduced from 5.59 to 2.36 kcal/mol; for the CB8 test system, the RMSD goes from 8.04 to 3.51 kcal/mol, while for the OAs test system, the RSMD goes from 2.89 to 0.95 kcal/mol. The scaling approach was inspired by studies on structurally related known benchmark sets: by simply scaling, the RMSD was reduced from 6.23 to 1.19 kcal/mol and from 2.96 to 0.62 kcal/mol for the CB8 benchmark system and the OA benchmark system, respectively. We find this scaling procedure to correct absolute binding affinities to be highly effective especially when working across a “congeneric” series with similar charge states. It is less successful when applied to mixed ligands with varied charges and chemical characteristics, but improvement is still realized in the present case. This approach suggests that there are large systematic errors in absolute binding free energy calculations that can be straightforwardly accounted for using a scaling procedure. Random errors are still an issue, but near chemical accuracy can be obtained using the present strategy in select cases.     

Improved interaction potentials for charged residues in proteins

Electrostatic interactions dominate the structure and free energy of biomolecules. To obtain accurate free energies involving charged groups from molecular simulations, OPLS-AA parameters have been reoptimized using Monte Carlo free energy perturbation. New parameters fit a self-consistent, experimental set of hydration free energies for acetate (Asp), propionate (Glu), 4-methylimidazolium (Hip), n-butylammonium (Lys), and n-propylguanidinium (Arg), all resembling charged residue side chains, including beta-carbons. It is shown that OPLS-AA free energies depend critically on the type of water model, TIP4P or TIP3P; i.e., each water model requires specific water-charged molecule interaction potentials. New models (models 1 and 3) are thus described for both water models. Uncertainties in relative free energies of charged residues are approximately 2 kcal/mol with the new parameters, due to variations in system setup (MAEs of ca. 1 kcal/mol) and noise from simulations (ca. 1 kcal/mol). The latter error of approximately 1 kcal/mol contrasts MAEs from standard OPLS-AA of up to 13 kcal/mol for the entire series of charged residues or up to 5 kcal/mol for the cationic series Lys, Arg, and Hip. The new parameters can be used directly in molecular simulations with no modification of neutral residues needed and are envisioned to be particular important in simulations where charged residues change environment.     

Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA.

In this work, a new ansatz is presented that combines molecular dynamics simulations with MM-PBSA (Molecular Mechanics Poisson-Boltzmann/surface area) to rank the binding affinities of 12 TIBO-like HIV-1 RT inhibitors. Encouraging results have been obtained not only for the relative binding free energies, but also for the absolute ones, which have a root-mean-square deviation of 1.0 kcal/mol (the maximum error is 1.89 kcal/mol). Since the root-mean-square error is rather small, this approach can be reliably applied in ranking the ligands from the databases for this important target. Encouraged by the results, we decided to apply MM-PBSA combined with molecular docking to determine the binding mode of efavirenz SUSTIVA(TM) another promising HIV-1 RT inhibitor for which no ligand-protein crystal structure had been published at the time of this work. To proceed, we define the following ansatz: Five hundred picosecond molecular dynamics simulations were first performed for the five binding modes suggested by DOCK 4.0, and then MM-PBSA was carried out for the collected snapshots. MM-PBSA successfully identified the correct binding mode, which has a binding free energy about 7 kcal/mol more favorable than the second best mode. Moreover, the calculated binding free energy (-13.2 kcal/mol) is in reasonable agreement with experiment (-11.6 kcal/mol). In addition, this procedure was also quite successful in modeling the complex and the structure of the last snapshot was quite close to that of the measured 2,3 A resolution crystal (structure the root-mean-square deviation of the 54 C(alpha) around the binding site and the inhibitor is 1.1 A). We want to point out that this result was achieved without prior knowledge of the structure of the efavirenz/RT complex. Therefore, molecular docking combined with MD simulations followed by MM-PBSA analysis is an attractive approach for modeling protein complexes a priori.     

Further development and validation of empirical scoring functions for structure-based binding affinity prediction

New empirical scoring functions have been developed to estimate the binding affinity of a given protein-ligand complex with known three-dimensional structure. These scoring functions include terms accounting for van der Waals interaction, hydrogen bonding, deformation penalty, and hydrophobic effect. A special feature is that three different algorithms have been implemented to calculate the hydrophobic effect term, which results in three parallel scoring functions. All three scoring functions are calibrated through multivariate regression analysis of a set of 200 protein-ligand complexes and they reproduce the binding free energies of the entire training set with standard deviations of 2.2 kcal/mol, 2.1 kcal/mol, and 2.0 kcal/mol, respectively. These three scoring functions are further combined into a consensus scoring function, X-CSCORE. When tested on an independent set of 30 protein-ligand complexes, X-CSCORE is able to predict their binding free energies with a standard deviation of 2.2 kcal/mol. The potential application of X-CSCORE to molecular docking is also investigated. Our results show that this consensus scoring function improves the docking accuracy considerably when compared to the conventional force field computation used for molecular docking.     

Predicting binding affinities of host-guest systems in the SAMPL3 blind challenge: the performance of relative free energy calculations

Relative free energy calculations based on molecular dynamics simulations are combined with available experimental binding free energies to predict unknown binding affinities of acyclic Cucurbituril complexes in the blind SAMPL3 competition. The predictions yield root mean square errors between 2.6 and 3.2 kcal/mol for seven host-guest systems. Those deviations are comparable to results for solvation free energies of small organic molecules. However, the standard deviations found in our simulations range from 0.4 to 2.4 kcal/mol, which indicates the need for better sampling. Three different approaches are compared. Bennett's Acceptance Ratio Method and thermodynamic integration based on the trapezoidal rule with 12 λ-points exhibit a root mean square error of 2.6 kcal/mol, while thermodynamic integration with Simpson's rule and 11 λ-points leads to a root mean square error of 3.2 kcal/mol. In terms of absolute median errors, Bennett's Acceptance Ratio Method performs better than thermodynamic integration with the trapezoidal rule (1.7 vs. 2.9 kcal/mol). Simulations of the deprotonated forms of the guest molecules exhibit a poorer correspondence to experimental results with a root mean square error of 5.2 kcal/mol. In addition, a decrease of the buffer concentration by approximately 20 mM in the simulations raises the root mean square error to 3.8 kcal/mol.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.