Accuracy Comparison of Generalized Born Models ... - ACS Publications

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2018, 14, 1656−1670

Accuracy Comparison of Generalized Born Models in the Calculation of Electrostatic Binding Free Energies Saeed Izadi,† Robert C. Harris,‡ Marcia O. Fenley,*,¶ and Alexey V. Onufriev*,§,∥ †

Early Stage Pharmaceutical Development, Genentech Inc., 1 DNA Way, South San Francisco, California 94080, United States Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, Maryland 21201, United States ¶ Institute of Molecular Biophysics, Florida State University, Tallahassee, Florida 32306-3408, United States § Department of Computer Science and Physics and ∥Center for Soft Matter and Biological Physics, Virginia Tech, Blacksburg, Virginia 24061, United States ‡

S Supporting Information *

ABSTRACT: The need for accurate yet efficient representation of the aqueous environment in biomolecular modeling has led to the development of a variety of generalized Born (GB) implicit solvent models. While many studies have focused on the accuracy of available GB models in predicting solvation free energies, a systematic assessment of the quality of these models in binding free energy calculations, crucial for rational drug design, has not been undertaken. Here, we evaluate the accuracies of eight common GB flavors (GB-HCT, GB-OBC, GB-neck2, GBNSR6, GBSW, GBMV1, GBMV2, and GBMV3), available in major molecular dynamics packages, in predicting the electrostatic binding free energies (ΔΔGel) for a diverse set of 60 biomolecular complexes belonging to four main classes: protein−protein, protein-drug, RNA-peptide, and small complexes. The GB flavors are examined in terms of their ability to reproduce the results from the Poisson−Boltzmann (PB) model, commonly used as accuracy reference in this context. We show that the agreement with the PB of ΔΔGel estimates varies widely between different GB models and also across different types of biomolecular complexes, with R2 correlations ranging from 0.3772 to 0.9986. A surface-based “R6” GB model recently implemented in AMBER shows the closest overall agreement with reference PB (R2 = 0.9949, RMSD = 8.75 kcal/mol). The RNA-peptide and protein-drug complex sets appear to be most challenging for all but one model, as indicated by the large deviations from the PB in ΔΔGel. Small neutral complexes present the least challenge for most of the GB models tested. The quantitative demonstration of the strengths and weaknesses of the GB models across the diverse complex types provided here can be used as a guide for practical computations and future development efforts.

■

INTRODUCTION

environment can therefore have a profound effect on the outcomes of receptor−ligand binding calculations.18−20 Traditionally, the effects of aqueous solvation have been accounted for by placing a sufficiently large number of individual water molecules around the solute and simulating their motion on an equal footing with the molecule of interest.21−23 While arguably the most realistic of the current theoretical approaches, it suffers from considerable computational costs, which often become prohibitive, for example for the receptor−ligand docking problem in which one searches for the best position of the ligand molecule (drug candidate) on the target receptor.24,25 The cost of such calculations increases dramatically with the system size. Moreover, it is not a priori clear at this point which, if any, of the numerous available explict water models can be used as an accuracy standard for estimates of receptor−ligand binding affinities. For example, the popular

Accurate determination of the binding affinities of small molecules with a biological target is crucial in molecular design and drug discovery but remains challenging.1−9 Experimental approaches for assessing the free energy of binding are difficult and prohibitively expensive on a large scale.10 Therefore, considerable effort has been devoted to develop in silico approaches with the aim of reducing time and cost in modern drug design processes, namely free energy perturbation (FEP),11 thermodynamic integration (TI),12 linear response (LR),13 and molecular mechanics-Poisson−Boltzmann-surface area (MM-PBSA)14,15 methods. Among all, molecular simulations have gained significant popularity in estimating binding affinities due to substantial methodological improvements and advanced computational power.14,16,17 In essence, the free energy of binding is largely determined from the combined effect of the delicate balance between receptor−ligand, solvent−ligand, solvent-receptor, and solvent−solvent interactions. The quality of the underlying solvent © 2018 American Chemical Society

Received: August 20, 2017 Published: January 29, 2018 1656

DOI: 10.1021/acs.jctc.7b00886 J. Chem. Theory Comput. 2018, 14, 1656−1670

Article

Journal of Chemical Theory and Computation

solvation free energy. We adopt the same strategy here in terms of the PB reference and compare the predictions of ΔGel and ΔΔGel by different GB models to those given by the numerical PB. That is we focus on the electrostatic component of the binding free energy, since it is the electrostatics that the GB approximation is designed to estimate. We fully acknowledge the approximate nature of the PB;108 further justification for the use of the numerical PB reference for the GB models considered here is provided in “Methods”. Implications for practical calculations of binding affinity are discussed in “Conclusions”. We expressly avoid comparing computational efficiencies of different models from multiple packages− interested readers are referred to previous comprehensive studies64,100 and, for the newest models, to the original manuscripts.109,110 The remainder of the paper is organized as follows. First, we present the specifics of the receptor−ligand complexes used for the comparative studies, along with a brief overview of the GB flavors and the reference PB model used in this study. Second, we present the quality of the calculated ΔGel and ΔΔGel with respect to reproducing the values predicted by the reference PB models. Our findings, including implications of the GB model accuracy relative to the PB in the context of practical calculations are summarized and discussed in the Conclusions.

TIP3P model can lead to 40% (6 kcal/mol) error relative to experiment in binding enthalpy estimates of even small host− guest systems;26 differences between binding free energy estimates based on the widely used TIP3P and TIP4P-Ew models can be as high as 9 kcal/mol.27 The implicit solvation methodology28,29,31−50 provides significant computational advantages by representing the discrete solvent molecules by a uniform continuous medium with the average properties of the solvent. Currently, the practical “engine” of the implicit solvent models−responsible for the estimation of the key electrostatic interactions−is either the Poisson−Boltzmann (PB) model14,29,38,45,51−57 or, in the case where computational efficiency is paramount, the generalized Born (GB) approximation.36,58−63 The GB, of which a great variety of flavors and implementations are available,63−92 has become popular due to a reasonable balance between accuracy and speed and the availability of a diverse class of GB implementations in major molecular modeling packages.93−97 While arguably most widely used in molecular dynamics simulations, the relatively high efficiency of the GB model makes it also very promising in large-scale estimates of solvation free energies and, hence, in estimates of receptor− ligand binding interactions,16,30,34,61,98,99 especially in conjunction with approaches, such as MM-PB(GB)SA,14−16 where large numbers of snapshots need to be processed. In the past, systematic studies were performed to examine the accuracy and efficiency of computing the electrostatic solvation free energies (ΔGel) with various implicit solvent models, including the GB and PB.18,19,64,100−103 However, even a reasonably good accuracy of a GB (or the PB, for that matter) model in reproducing the electrostatic solvation free energy alone may not necessarily translate into the same level of accuracy for ligand binding calculations.56,63,104 One reason is that the values of binding free energies are relatively much smaller than the free energy terms corresponding to the bound and unbound states (e.g., tens of kcal/mol for binding affinities versus thousands of kcal/mol for solvation free energies), which leads to high sensitivity of the method to parameters used in the calculations of binding affinities.56,105−107 To the best of our knowledge, no systematic assessment of the GB accuracy in receptor−ligand binding was ever performed on a large and diverse set of complexes, simultaneously for most GB models available in major modeling packages. This study is intended to fill this critical gap. Here, we study the accuracies of the electrostatic binding free energies (ΔΔGel) computed for a set of 60 biomolecular complexes using eight commonly used GB models available in the major molecular modeling packages AMBER93 and CHARMM.94 Several of these GB models originally developed for these packages have since been adopted by other major modeling packages.95−97 The sets of biomolecular complexes tested here are diverse with respect to the biomolecule type (e.g., nucleic acids, proteins, peptides, and drugs), structure sizes, and net charges, as well as ΔΔGel values. The objective is to provide a quantitative assessment of the accuracies of common GB models across diverse complex types. All of the GB models (flavors) considered here have been developed in the context of simulating biological macromolecules; the PB model has been used as their natural accuracy reference, including previous large-scale estimates of the GB performance across multiple models64,100 focused on

■

METHODS Test Structures. Four test sets containing 60 biomolecular complexes were collected for this study: a set of small complexes (data set 1), a set of protein-drug complexes (data set 2), a set of protein−protein complexes (data set 3), and a set of RNA-peptide complexes (data set 4). The structures of these molecules were taken from the RCSB protein databank111 (PDB), and the PDB identification codes of these structures can be found in Table 1. These sets are diverse with respect to the biomolecule type (e.g., nucleic acids, proteins, and drugs), structure sizes, and net charges (see Table 2), as well as their calculated ΔΔGel (see the Results section). The structures in set 1 (small neutral complexes) were originally selected27,112 to contain no more than 2000 heavy atoms and no missing heavy atoms. They were protonated as described previously.27 Table 1. PDB ID of the Complexes in the Four Data Sets

1657

data set 1 (small complexes)

data set 2 (proteindrug)

1B11 1BKF 1F40 1FB7 1FKB 1FKF 1FKG 1FKH 1FKJ 1FKL 1PBK 1ZP8 1ZPA 2FKE 2HAH 3FKP

1JTX 1JTY 1JUM 1JUP 1QVT 1QVU 3BQZ 3BR3 3BTC 3BTI 3BTL 3PM1

data set 3 (protein−protein) 1B6C 1BEB 1BVN 1E96 1EMV 1FC2 1GRN 1HE1 1KXQ 1SBB 1FMQ 1UDI 2BTF 2PCC 2SIC 2SNI 7CEI

data set 4 (RNApeptide) 1A1T 1A4T 1BIV 1EXY 1HJI 1I9F 1MNB 1NYB 1QFQ 1ULL 1ZBN 2A9X 484D


Article


errors in ΔΔGel, which is the only quantity that the GB model predicts. Reference Poisson−Boltzmann Calculations. The finite-difference solution of the 3D linear PB equation was solved using the Adaptive Poisson−Boltzmann Solver (APBS).120 To obtain the ΔGel of a molecule the linear PB equation was solved twice, first in salt solution and then in vacuum, and then the difference between the two electrostatic free energies was taken. Each molecule was treated as a low dielectric region with a dielectric constant (ϵin) of 1. The solvent (here a salt solution) region was treated as a high dielectric region with a dielectric constant (ϵout) of 80 and NaCl concentration of 0.145 M, with the latter representing a physiological environment. The solvent temperature was set to 298.15 K. The dielectric boundary that separates the low dielectric molecule from the solvent region was defined using the solvent excluded molecular surface obtained with a 1.4 Å spherical water probe radius. No ion exclusion region was considered for consistency with the GB models here examined. We used an autofocusing protocol with a fine grid 20 Å larger than the molecule in each dimension, and a coarse grid 1.7 times the size of the molecule in each dimension, and a fine grid spacing of 0.3 or 0.5 Å (to check for convergence with respect to grid spacing105,107). The default values of all other APBS parameters were used. The Generalized Born Models. The “canonical” generalized Born model approximates the Green function of the Poisson equation via a simple closed form expression, first introduced by Still et al.,60 and later augmented to account for the presence of monovalent salt121

Table 2. Average Sizes (in Atoms) and Average Absolute Total Charges of the Complexes and Ligands for the Four Test Sets average size (atoms) small complexes protein-drug protein−protein RNA-peptide

average |charge| (e)

complex

ligand

complex

ligand

1769.9 6203.4 6147 1209.5

101.4 44.9 2475.1 393.5

0.00 3.615 6.00 16.615

0.00 1.154 4.667 7.769

For the proteins and RNA molecules (except for data set 1), the atomic charges and hydrogens were assigned using the PDB 2PQR utility.113,114 The ff99bsc0 parameters and the GAFF force field, both part of AMBER,115 were used for preparing the topology and coordinate files, which includes partial charges. The atomic charges for the drugs in data set (2) were assigned with the methods outlined elsewhere.63 For the sake of consistency, we used the Bondi radii set for all of the implicit solvent models studied here, as these radii are recommended for specific AMBER GB models.27 However, to test the sensitivities of the calculated ΔGel and ΔΔGel to the choice of atomic radii, some of the calculations were repeated with the atomic radii other than Bondi. The PQR files values for the test sets are available in the Supporting Information. The Thermodynamic Cycle for Electrostatic Binding Free Energies. The thermodynamic cycle shown in Figure 1 is

qiqj exp( −κfGB ) ⎞ 1⎛ 1 ⎟∑ ΔGel = − ⎜ − 2 ⎝ ϵ in ϵout ⎠ ij fGB (rij , R i , R j)

where rij is the distance between atoms i and j, qi and qj are the respective (partial) charges, the Ri and Rj are the so-called effective Born radii, and κ is the Debye−Hückel screening parameter used to incorporate121 the electrostatic screening effects of (monovalent) salt. While alternatives have been proposed,71,89,122,123 the most widely used functional form of f GB is

Figure 1. Thermodynamic cycle used to compute the electrostatic component of the binding free energy. The surrounding dielectric medium is shaded for water and is white for vacuum.

fGB (rij , R i , R j) = [rij2 + R iR j exp( −rij2/4R iR j)]1/2

(3)

Apart from a “dielectric” correction to eq 2 used by one of the models (see below), the GB Green function eq 2 is the same for all of the GB models considered in this work. That Green function can be rigorously derived from the Poisson equation under a certain set of assumptions,36 which further supports124 the use of the numerical PB for testing of the GB accuracy, at least in the context of the GB flavors tested here. The latter were developed primarily for modeling of biological macromolecules. At the same time, we acknowledge the existence of a view108 that the GB approximation in general is not a direct approximation to the PB, see e.g. ref 124 for a discussion. The way the effective Born radii (Ri) are calculated differ significantly between implementations of the GB models; in that aspect, the GB models (flavors) that are investigated here represent a variety of classes. Full details of these GB models are provided in the original papers; here, we briefly review the main differences between them.

used to calculate ΔΔGel: first, the receptor and the ligand are transferred from solvent into vacuum, second, the complex is formed in vacuum, and then the complex is transferred from vacuum into the solvent. The ΔGel of the receptor, the ligand, and the complex is designated as ΔGRel, ΔGLel, and ΔGRL el , respectively. The electrostatic component of the binding free energy in vacuum is ΔECel, which is calculated similarly for all models. Using this cycle ΔΔGel = (ΔGelRL − ΔGelR − ΔGelL) + ΔEelC

(2)

(1)

In all the calculations presented here we assume no conformational changes to the structures and do not consider the nonpolar component of the solvation free energy or the entropic contributions to binding. While estimating these quantities would be necessary for a direct comparison with experiment, the associated uncertainties in their theoretical estimates can be large9,14,15,104,116−119 and may obfuscate the 1658


Article


also implemented in the MM_PBSA137,138 module and the Free Energy Workflow (FEW)139 tool of AmberTools17.136 In the GBNSR6 model, the effective Born radii are computed via the “R6” numerical integration over the molecular surface of the solute, defined by the solute intrinsic radii and the solvent probe. GBNSR6 includes a single additional (to the PB) parameter: a constant offset B to the inverse effective Born radii. The value of the offset is fixed for a given probe radius; here we used B = 0.028 Å−1 as proposed by Mongan et al.135 for the probe radius of 1.4 Å. That is, as long as the solvent probe radius is set to the default 1.4 Å, and agreement with the PB is the goal, the GBNSR6 model becomes parameter-free in the same sense as the PB. Note that in addition to input parameters that are generally used by the PB, such as the intrinsic atomic radii, other GB models rely on multiple additional parameters to approximate the PB energies; the optimum values of these parameters generally depend on the atomic radii. By default, ΔGel in GBNSR6 is calculated by the ALPB model,134 which introduces a physically correct dependence on the dielectric constants into the original GB model (eq 2), while maintaining the efficiency of the original

While the GB does not employ the dielectric boundary in exactly the same manner as the PB model does, the solute/ solvent boundary plays the same conceptual role in the GB; technically, the dielectric boundary enters125 the GB model through the effective Born radii. In all of the models tested here the effective Born radii are computed as approximate or exact integrals over the suitably chosen solute volume or surface. The two-dielectric model considered here is defined by specifying the values of the internal (solute) and external (solvent) dielectrics and the postion of the sharp boundary between them. Thus, within this model it is natural to assume the same dielectric constants for the GB when testing it directly against the PB. Within this limited scope testing, the precise value of the dielectric constant used to describe the low dielectric interior of the macromolecule is not very relevant. However, we argue that using the lowest physically meaningful value of εin = 1 is best if the goal is to highlight the GB versus PB error−the error is maximized when the solute and solvent dielectrics differ the most and disappears when the two are equal.126 For all the GB models, we used interior and exterior dielectric constants of (εin = 1) and (εout = 80), respectively, a NaCl concentration of 0.145 M, the temperature set to 298.15 K, and a solvent probe radius of 1.4 Å where appropriate. In order to correspond with the reference PB calculations, we used the same intrinsic atomic radii for the GB and PB calculations; see below. Most of the reported results are based on Bondi radii, but we have explored the sensitivity of the conclusions to several other radii sets as described below. Appropriate radii were incorporated into the prmtop files using a variation of the AMBER2CHARMM tool. AMBER GB Models. Four variants of the GB model were used: GB-HCT,127−129 GB-OBC,130,131 GB-neck2,88,132 and GBNSR6.109,133 Except for GBNSR6, all of these GB models are based on the Coulomb field approximation and differ mainly in how the effective Born radii are computed. The analytical linearized PB (ALPB) correction134 was used with the GBNSR6 model only. GB-HCT uses the pairwise descreening approach of Hawkins, Cramer, and Truhlar, 127,128 with parameters described by Tsui and Case,129 which leads to a GB model implemented in Amber with igb = 1. The 3D integral used in the estimation of the effective radii is performed over the van der Waals spheres of solute atoms, which implies a definition of the solute volume in terms of a set of spheres, rather than the more complex molecular surface. In GB-OBC,131 the effective Born radii are rescaled to account for the interstitial spaces between atom spheres missed by the GB-HCT approximation. In that sense, GB-OBC is intended to be a closer approximation to the true molecular volume, albeit in an average sense. GB-OBC can be used in AMBER by setting igb = 5 or 2; here we used the igb = 5 version. GB-neck2 is the latest, most successful version of an earlier GBn model,135 with different additional parameters.88,132 GBneck2 uses a pairwise correction term to GB-OBC to approximate135 a molecular surface dielectric boundary; that is to eliminate interstitial regions of high dielectric smaller than a solvent molecule. This correction affects all atoms and is geometry-specific, going beyond the geometry-free, “average” rescaling approach of GB-OBC, which mostly affects buried atoms. GB-neck2 can be used in AMBER by setting igb = 8. GBNSR6109,133,136 (“Numerical Surface R6 GB”) is one of the latest additions to the AMBER family of GB models, now

1⎛ 1 1 ⎞ 1 ΔGel = − ⎜ − ⎟ 2 ⎝ ϵ in ϵout ⎠ 1 + βα

⎛ 1 αβ ⎞ ⎟ + GB A⎠ ⎝f

∑ qiqj⎜ ij

(4)

where β = ϵin/ϵout, α = 0.571412, and A is the electrostatic size of the molecule, which is essentially the overall size of the structure; A is estimated analytically.134 CHARMM GB Models. Four GB models available in the CHARMM software package94 were evaluated in this study: GBMV1 and GBMV2 are analytical methods in which the molecular volume is constructed from a superposition of atomic functions.73,140 GBMV1 is the oldest and least accurate of these two analytical GB models. GBMV373,100 is a grid-based method that uses a Cartesian grid representation of the molecular volume similar to that used in traditional Poisson−Boltzmann calculations, and GBSW87,141 uses a van der Waals-based surface with a simple polynomial switching function to smooth the dielectric boundary. A more detailed description of these models is available in a recent review.100 To compute these energies, CHARMM topology and coordinate files were created with a script. The notation used here for naming the different CHARMM GB parameters follows the CHARMM manual (see the documentation for the gbmv and gbsw modules). For GBSW, the molsurf option was used to mimic the molecular surface, conc = 0.145 M, sw = 0.2, and the parameters aa0 and aa1 were set to 1.2045 and 0.1866, respectively. For GBMV1 the following parameters were set: beta to −100, the shift parameter to −0.5, lambda1 to 0.1, p1 to 0.44, tt to 2.92, corr to 0, bufr to 0.5, mem to 20, cuta to 20, p3 to 0, wtyp to 0, and eshift to 0. For GBMV2 the following parameters were set: tol to 10−10, mem to 20, cuta to 20, dn to 1, bufr to 0.2, beta to −12, shift to −0.1, slope to 0.9, lambda1 to 0.5, p1 to 0.45, p2 to 1.25, p3 to 0.65, p6 to 8, onx to 1.9, offx to 2.1, corr to 1, son to 1.2, soff to 1.5, nphi to 5, sa to 0, sb to 0, and kappa to 0.1239950497. For GBMV3 the following parameters were set: dn to 0.2, p6 to 4.0, nphi to 38, shift to −0.007998, slope to 0.9026, wtyp to 2, and kappa to 0.1239950497. All other parameters used their default values. Choice of the Input Atomic Radii Set. Within the implicit solvation framework, the solvation energies are 1659


Article

Journal of Chemical Theory and Computation Table 3. Deviation of GB ΔGel Values (kcal/mol) from the PB Reference for the Four Data Setsa data set 1 (small complexes) RMSD GBNSR6 GB-HCT GB-OBC GB-neck2 GBSW GBMV1 GBMV2 GBMV3

GBNSR6 GB-HCT GB-OBC GB-neck2 GBSW GBMV1 GBMV2 GBMV3

21.8 33.6 24.7 6.7 11.5 135.4 14.0 25.0

R2

0.999987 18.6 0.999703 −25.3 0.999904 20.8 0.999919 −0.2 0.999799 5.4 0.999537 −109.2 0.999964 12.0 0.999945 20.7 data set 3 (protein−protein)

data set 2 (protein-drug)

RMS of worst 20%

RMSD

31.0 52.5 36.8 12.5 20.0 201.2 23.3 37.3

21.4 260.6 172.8 36.1 148.0 895.7 19.9 80.7

R2

0.999974 9.3 0.999938 −211.1 0.999764 134.7 0.999902 1.5 0.999937 −119.1 0.999854 −731.9 0.999971 2.5 0.999954 62.6 data set 4 (RNA-peptide)

RMS of worst 20% 42.7 358.3 271.8 70.7 216.8 1164.2 37.1 131.0

RMSD

R2

RMS of worst 20%

RMSD

R2

RMS of worst 20%

39.5 358.2 93.0 101.4 112.8 701.0 39.0 82.1

0.999944 0.996962 0.999094 0.999391 0.999709 0.995759 0.999786 0.999759

33.8 −275.7 −26.5 −62.2 −89.8 −610.7 13.0 64.3

66.9 674.3 193.4 205.6 211.8 1167.1 71.5 148.5

14.5 141.7 351.7 153.5 67.5 176.6 59.8 44.2

0.999997 0.999835 0.999343 0.999701 0.999902 0.999571 0.999947 0.99996

11.1 95.0 276.1 −111.06 −46.5 −144.1 34.0 −28.0

25.3 256.8 578.0 292.8 126.3 301.5 111.4 85.9

a

Error metrics calculated for each set of predictions include root mean square deviation (RMSD), average signed error (), square of the Pearson correlation coefficient (R2), and RMSD of the 20% of complexes with the worst agreement with the reference.

extremely sensitive to the intrinsic atomic radii values,107,125,142 because these determine the position of the dielectric boundary between the low dielectric interior and the high dielectric solvent. Both the GB and the PB models take the intrinsic atomic radii as input−many different radii sets exist.143−147 Thus, one has to carefully consider the choice of the input radii set for any accuracy analysis of these models, including the one presented here. As we shall see below, the appropriate choice is not obvious and depends on the kind of comparison one wants to make and the specifics of the models used. In principle, one could consider using different radii sets for different models; such strategy would make sense if the accuracy metric were an experimental observable such as the binding free energy. In this case, one could choose a set of parameters for each model that gives best agreement with experiment and then compare the outcomes to determine the most accurate model, GB or PB, with respect to the chosen metric. However, for reasons outlined above, our goal here is more limited and focused: to evalute only the electrostatic component of solvation energy relative to the PB reference. We argue that in this case using, for all of the GB models, one and the same radii set employed to obtain the PB reference energies makes sense on the grounds of physical consistency. Indeed, all of the GB models reduce to the Born formula for a single spherical ion, which also corresponds to the exact solution of the Poisson equation within the twodielectric formalism shared by all of the models considered here. Thus, it makes sense that all the models, GB or PB, are required to give identical results for a single isolated spherical ion in water, which implies the same radii set. Here we propose to use the well established Bondi radii set: it is physically wellgrounded and has only 5 main atom types, all of which bodes well for low probability of an unwanted bias toward a specific training set and type of reference used in radii set optimizations. Note that more complex radii sets optimized against hydration free energies of small molecules or amino acids are more likely to be poorly transferable to computation of protein−ligand binding energies.148 Unfortunately, even the

seemingly straightforward choice of Bondi radii needs to be carefully justified for the diverse set of practical GB models tested here. This is because most of these GB models contain multiple fitting parameters used to optimize the agreement with the respective accuracy reference, which most often includes the numerical PB solvation energies. These GB models are intended to be used with specific atomic radii set employed during the parameter optimization, e.g. modifications of Bondi set mbondi, mbondi2, and mbondi3 for GB-HCT, GB-OBC, and GB-neck2, respectively (AMBER), and CHARMM radii for many CHARMM GB models. Among all of the GB models tested here, GBNSR6 is the only model with no recommended radii set, because the model is essentially agnostic to the input intrinsic radii in the same sense as the PB reference is. For all the other GB models, changing the input atomic radii away from the recommended set without an extensive reparametrization effort−which is out of scope here−may result in deterioration of the accuracy of the model relative to the PB reference values based on the same radii set; a possible appearance of accuracy improvement for some structural classes would likely be due to a fortuitous cancellation of error for this class of structures only. To address this subtle issue, we have confirmed that the computed ΔΔGel values are not sensitive to the use of Bondi set instead of the recommended set for the GB models where Bondi set was not listed as one of the recommended options; see the ”Results” section for details. Perfect Effective Born Radii. To compute the perfect122 effective Born radii Ri, a Poisson problem is set up and solved separately for each atom using the numerical PB solver PEP,30 which is especially well suited for the task. In the process, the dielectric boundary of the full molecule is present, but only the charge of that particular atom is nonzero (all other charges are set to zero). The process yields the self-energy ΔGiiel for each atom i, from which the perfect effective Born radius Ri is obtained by inverting the Born equation ΔG elii = 1 − 2 (1/ϵ in − 1/ϵout)R i−1. We set solute dielectric ϵin = 1 and use ϵout = 1000 for the solvent to mimic the conductor limit as 1660


Article


Figure 2. a) The root-mean-square deviations (RMSD) between the electrostatic solvation free energies (ΔGel) predicted by the different generalized Born (GB) models and those predicted by the reference PB for different protein−ligand complex sets. b) The RMSD between the electrostatic binding free energies (ΔΔGel) predicted by the different GB models and those predicted by the reference PB for different complex sets.

discussed previously.126 Here, the perfect radii are used only to further explore the accuracy of GBNSR6; for consistency with the GBNSR6 results without the perfect radii, the corresponding ΔGel (perfect radii) values in Table 6 are computed using eq 4. The corresponding PB reference values of ΔΔGel reported in Table 6 are computed using APBS (rather than PEP) for consistency with the rest of the paper; we verified that these values change insignificantly if computed with PEP instead.

given by the different GB models from those given by the PB equation. The root-mean-square deviation (RMSD), average signed error (), square of the Pearson correlation coefficient (R2), and RMSD of the 20% of complexes with the worst agreement with the reference ΔGel were calculated for each of the four data sets. Overall, the ΔGel values from all the GB models highly correlate with the values from the PB equation. The lowest correlation is for GBMV1 for data set 3 (R2 = 0.9957), but even this correlation is good. Although the R2 between ΔGel given by the PB equation and that given by the different GB models was good for each GB model, the RMSD and average deviations are significant, and significantly different, between different GB models (Table 3 and Figure 2(a)). For instance, while the RMSD of GBNSR6 from the PB is smaller than 39.5 kcal/mol over all four sets, the RMSD of GBMV1 from the PB can be as large as 895.7 kcal/ mol (data set 2). That GBMV1 is less accurate than GBMV2 agrees with similar GB calculations done elsewhere.73 The RMSD of GB-OBC varies significantly over the four sets, from 24.7 kcal/mol for data set 1 to 351.7 kcal/mol for data set 4. Among the AMBER GB models GBNSR6 has the smallest RMSD overall (all four sets), and among the CHARMM GB models GBMV2 had the smallest overall RMSD from the PB; the largest RMSDs in ΔGel for these models are 39.0 kcal/mol (data set 3) and 59.8 kcal/mol (data set 4), respectively. However, for one structure set, namely set 1 (small neutral complexes), it is GB-neck2 that shows the best agreement with the PB. The nonzero average deviations (Table 3) point to the existence of noticeable systematic shifts between the estimates of ΔGel obtained with the PB equation and those obtained with the GB models, which can be due to a number of factors, including the following: (i) Differences in the dielectric boundary definitions used by each model. The PB model typically employs a grid-resolution dependent numerical approximation to the molecular surface, while the analytical GB models mimic this surface to varying degrees of accuracy. The exception is GBNSR6−the version used here employs essentially the same molecular surface, but still subtle details, including surface discretization, differ between it and what is used by the PB solver used here. (ii) The GB Green function is only an approximation to the correct Green function of the Poisson problem for a given dielectric boundary. The errors introduced by the approximate Green function can be tested via

■

RESULTS Poisson−Boltzmann Convergence. To determine whether our PB results were converged with respect to grid spacing, we computed 20 estimates each of ΔGRel, ΔGLel, and ΔGRL el at grid spacings of 0.5 Å with the centers of the solution grid shifted by a random fraction of a grid spacing in each direction and averaged these results. We then compared these numbers to similar numbers obtained at a grid spacing of 0.3 Å. By averaging over these solution grids we attempted to eliminate the error due to grid placement. The root-meansquare deviation (RMSD) between the APBS results with grid spacings of 0.3 Å and those with grid spacings of 0.5 Å are 12.6, 7.9, 2.7, and 2.8 kcal/mol for the RNA-peptide, protein− protein, protein-drug, and small complexes data sets, respectively. These RMSDs between the APBS calculations with grid spacings of 0.3 Å and those with grid spacings of 0.5 Å may at first appear to be large compared to those between the GB models and the APBS calculations at 0.3 Å, but the energies with grid spacings of 0.3 Å are always systematically smaller than those with a grid spacing of 0.5 Å because of the bias present in such calculations.105,107 The correlations between the two sets of APBS calculations are therefore much higher than would be expected from the RMSD alone. The square (R2) of the Pearson’s correlation coefficients between ΔGel given by a grid spacing of 0.3 Å and those given by a grid spacing of 0.5 Å was greater than 0.999 for all four data sets. For ΔΔGel R2 was greater than 0.992 for all four data sets. For any of the GB models, weaker correlations can be attributed to errors in the GB model. Electrostatic Solvation Free Energies (ΔGel). In this section, we briefly review the accuracies of the different GB models in predicting ΔGel for the four sets of complexes (including the complex and receptor and ligand components). Table 3 and Figure 2(a) summarize the deviations of the ΔGel 1661


Article


Figure 3. continued

1662


Article


Figure 3. Correlation between electrostatic binding free energies (ΔΔGel) computed by different GB models relative to the PB. The diagonal lines are x = y.

components, large systematic shifts in ΔGel can cancel out in the binding free energy calculations. In general, we find that GB models that show lower RMSD with respect to the PB ΔGel

perfect (PB-based) effective Born radii, as discussed below. One should notice that since binding free energies depend on the “difference” in ΔGel values of the receptor and ligand 1663


Article

Journal of Chemical Theory and Computation Table 4. Deviation of GB ΔΔGel Values (kcal/mol) from the PB Reference for the Four Data Setsa data set 1 (small complexes) RMSD GBNSR6 GB-HCT GB-OBC GB-neck2 GBSW GBMV1 GBMV2 GBMV3

GBNSR6 GB-HCT GB-OBC GB-neck2 GBSW GBMV1 GBMV2 GBMV3 a

3.9 10.7 7.3 6.2 4.0 9.6 5.0 3.0

R2

data set 2 (protein-drug)

RMS of worst 20%

0.979509 −3.6 0.876606 −10.2 0.946172 −6.7 0.910121 −5.5 0.924397 3.0 0.463363 6.9 0.935667 −4.2 0.913558 0.2 data set 3 (protein−protein)

RMSD

6.71 15.4 10.6 10.6 6.61 15.6 7.7 6.0

4.2 27.7 23.0 13.5 12.8 9.5 7.0 24.4

R2 0.949986 0.764202 0.377215 0.776368 0.774716 0.822145 0.88332 0.638866 data set 4

1.2 −26.0 −19.6 −10.8 −9.53 −5.81 1.8 −20.0 (RNA-peptide)

RMS of worst 20% 5.9 40.3 37.0 24.3 22.2 16.5 11.4 39.2

RMSD

R2

RMS of worst 20%

RMSD

R2

RMS of worst 20%

13.6 50.8 26.2 24.0 19.0 49.1 25.9 16.8

0.985424 0.921176 0.927651 0.952807 0.948775 0.94559 0.99063 0.981712

−9.2 −41.8 −11.9 −11.1 4.3 −44.2 −24.6 12.4

24.7 87.1 46.8 43.9 32.2 71.4 36.4 30.0

7.7 100.5 57.3 65.7 69.6 116.9 80.9 48.2

0.998602 0.91455 0.950862 0.95664 0.892278 0.920653 0.921114 0.924991

4.7 −89.2 −45.0 −56.5 −43.8 −107.0 −64.5 −10.5

12.4 148.1 87.8 99.1 95.4 152.6 123.3 86.3

The error metrics are described in Table 3.

values are also in better agreement with the PB ΔΔGel values, as discussed below. Electrostatic Binding Free Energies (ΔΔGel). Figure 2(b) and Table 5 (bottom) summarize the deviation of different GB models from the PB reference in estimating ΔΔGel, while the breakdown by individual structure test sets is given in Figure 3 and Table 4. We calculate RMSD, average signed error , R2, and the RMSD of the 20% of complexes with the largest deviations between ΔΔGel given by the PB equation and those given by the GB model. Overall, a wide range of accuracies is seen across the GB models tested, as judged by the different accuracy metrics. Some models (GBNSR6 and GBMV3) have very little overall systematic deviation from the PB, Table 5 (bottom), while others deviate substantially, e.g. the older GB-HCT and GBMV1 models. By far the lowest RMS error and the highest R2 correlation are seen for the GBNSR6 model. With regards to specific structures, the GB models also demonstrate a wide range of accuracies. For data set 1 (small, net neutral complexes), all of the GB models show good correlation with the reference PB values, except for GBMV1 for which R2 is 0.46. The low correlation of GBMV1 with the PB ΔΔGel values in set 1 is consistent with the large RMSD and average error of the ΔGel values (RMSD = 135.4 kcal/mol, Table 3). The average error in ΔΔGel for all of the AMBER GB models (GBNSR6, GB-HCT, GB-OBC, and GB-neck2) is negative (also see Figure 3), suggesting that the deviation of AMBER GB models from the used PB reference contains a systematic component for this test set. For data set 1, the best correlations with the PB reference are seen for the GBNSR6 (R2 = 0.98) and GBMV3 (R2 = 0.91) predictions. Overall, the lowest R2 correlation with the PB ΔΔGel across different GB models is seen for data set 2 (protein-drug). Perhaps not surprisingly, GB-OBC, which is among the oldest commonly used general-purpose GB models, yields a very low correlation with the PB (R2 = 0.38). The correlation (R2) is only 0.64 for GBMV3, which is surprising for a model shown to be superior to other CHARMM GB models.73,140 Only

GBNSR6 maintains a high correlation with the PB for this data set, (R2 ≈ 0.95) (Table 4 and Figure 3), suggesting that challenging structures are less forgiving to the approximations made by multiparameter GB models; among the tested models, GBNSR6 makes the fewest approximations relative to the PB. On average, the best agreement with the PB across all of the GB models is achieved for data set 3 (protein−protein). The lowest correlation with the PB is again for the oldest models, GB-HCT and GB-OBC, with (R2 ≈ 0.92), while R2 is as good as ≈0.99 for GBMV2 and GBNSR6. GBNSR6 also yields the smallest RMSD (13.6 kcal/mol). The good agreement with the PB for this set is not unexpected given that for most of the GB models protein structures were always a large part of the training set against which parameters of the models were fit. Data set 4 (RNA-peptide) appears to be the most challenging for all the GB models, with the exception of GBNSR6, with respect to RMSD from the PB (Table 4 and Figure 2). The average R2 is the second worst among all of the data sets, after data set 2. Again, GBNSR6 is an exception for this set: it shows a relatively small deviation from the PB (RMSD = 7.7 kcal/mol, R2 = 0.999). Overall, the large average deviation (and RMSD) from the PB reference in set 4 and the lowest deviation for set 1 correlate well with the average net charge of these two data sets, which carry the largest and smallest average net charges, respectively (Table 2). Our explanation for the correlation is as follows: errors in the GB Green function are scaled by the partial atomic charges, eq 2. Note that all structures in set 4 contain fragments of nucleic acids (RNA), which carry high partial charges on the phosphate groups, absent from all other sets. The correlation with the PB is worst for set 2 (not set 4, with the highest average charge), which suggests a more nuanced story. As discussed above, we conjecture that the GB Green function performs the worst for set 2, leading to the poorest correlation with the PB, but the resulting error values (RMSD and average) are not as high because the structures are not as highly charged (in the limit of zero partial charges, even very 1664


Article


Figure 4. Correlation between the accuracy of the GB models in predicting electrostatic binding affinities (ΔΔGel) and their accuracy in predicting electrostatic solvation free energy components (ΔGel) over all of the four complex sets: a) RMSD in ΔΔGel versus RMSD in ΔGel, b) R2 correlations in ΔΔGel versus R2 correlations in ΔGel, and c) R2 correlations in ΔΔGel versus RMSD in ΔGel.

Table 5. Overall Deviation of the GB ΔGel and ΔΔGel Values (kcal/mol) from Those Computed with the Numerical PB Equation for 60 Complexesa

low correlation with the PB can still result in small average and RMSD errors). Correlation between Accuracy in ΔΔGel and ΔGel Calculations. In general, we find a moderate correlation between the accuracies of ΔGel and those of ΔΔGel, Figure 4. Namely, when comparing a group of GB models that agrees most with the PB to the group that agrees the least, there is a clear correlation between the accuracy of ΔGel and ΔΔGel between these broad groups (but not within each group), whether the accuracy is measured by RMSD or R2 (Figure 4 (a) and (b)). Thus, a very poor agreement of ΔGel with the PB suggests that a good agreement of ΔΔGel is unlikely. Even though the range of R2 in GB ΔGel is very narrow (all R2 values are fairly close to 1 for all the GB models), Figure 4 (b), there is still a correlation between R2 of ΔGel and R2 of ΔΔGel. This correlation suggests that, when interpreting R2 for ΔGel, decimal numbers up to the fourth or even fifth place can be important for distinguishing between the accuracies of different GB models in terms of their ability to predict ΔΔGel. In addition, Figure 4 (c) shows that the accuracy of ΔGel and ΔΔGel measured by RMSD and R2, respectively, also correlate to some extent (Figure 4 (c)). Given that RMSD in ΔGel from the PB can significantly differ across GB models (Table 3), the data shown in Figure 4 (c) suggests that, for GB models on the “opposite ends” of the accuracy spectrum, RMSD can be used instead of R2 to distinguish between the models. One cannot reliably predict the accuracy in ΔΔGel based on that of ΔGel for models where the difference is not very large. For example, while the RMSD in ΔGel for GBMV2 (37.0 kcal/ mol) is smaller than that of GBMV3 (63.4 kcal/mol); the RMSD in ΔΔGel for GBMV2 is larger (40.4 kcal/mol) than that of GBMV3 (26.8 kcal/mol) (Table 5). Another example is the GB-OBC model which is the third least accurate model in ΔGel values (RMSD = 189.8 kcal/mol), yet it is the third most accurate model in ΔΔGel values (RMSD = 32.3 kcal/mol, R2 = 0.956) (Table 5). Source of the Remaining Error in GBNSR6. The error of the GB relative to the exact solution of the corresponding Poisson problem can be decomposed into two independent components: the error in approximating the effect of the same (in our case, molecular surface) dielectric boundary and the error due to the approximate GB Green function.125 The first error can be eliminated by the use of the perfect, that is PB-based effective Born radii−while computationally ineffective, these effective radii ensure that the self-terms in eq 2 and the

ΔGel GBNSR6 GB-HCT GB-OBC GB-neck2 GBSW GBMV1 GBMV2 GBMV3

GBNSR6 GB-HCT GB-OBC GB-neck2 GBSW GBMV1 GBMV2 GBMV3 a

2

RMSD

R

27.2 240.5 189.8 92.1 97.9 577.0 37.0 63.4

0.9999 0.9945 0.9974 0.9994 0.9995 0.9862 0.9999 0.9996

RMSD

R2

RMS of worst 20%

8.7 56.2 32.3 34.0 34.6 61.0 40.4 26.8

0.994889 0.89914 0.955511 0.949725 0.91599 0.804674 0.913058 0.953638

−2.4 −40.2 −19.4 −19.4 −9.4 −35.9 −22.1 −2.8

17.3 113.2 66.0 70.8 73.6 191.2 69.7 54.4

RMS of worst 20%

19.5 −114.6 86.6 −42.4 −61.4 −402.1 15.0 32.3 ΔΔGel

50.2 471.8 395.9 198.1 194.2 1123.3 75.2 123.4

The error metrics are described in Table 3.

corresponding PB values are identical, so that the remaining error is only due to the cross terms, i.e. the GB Green function itself. As the errors of GBNSR6 on our test sets is the lowest, we chose this model to explore the origin of the remaining error. Due to computational constraints, we focused on set 1− the set of smallest structures among the four sets used here. We calculated perfect effective Born radii63,122 for the two complexes from set 1 that showed the largest deviation of GBNSR6 ΔΔGel values from the PB reference. The results, Figure 5, show that, for the two complexes tested here, most of the effective Born radii from GBNSR6 correlate well with the perfect (PB) radii, with the discrepancies seen mostly for “buried” atoms. We attribute most of the discrepancy to subtle differences between the dielectric boundary employed by this version of GBNSR6 (the dielectric boundary is the molecular surface calculated by MSMS149) and the approximation of the boundary employed by typical PB solvers; see ref 136 for further details. One may excpect that reducing/eliminating the 1665


Article


Figure 5. Comparing the inverse effective Born radii calculated using GBNSR6 and PB (perfect radii) for two complexes from data set 1: a) 1bkf and b) 1fkl. The diagonal lines are x = y.

with the reference, so that if a “native” radii set were to lead to an improvement, we would most likely see it on this set. We used the same logic to also evaluate the sensitivity of GB-neck2 to input radii. Here we used the other “worst performing” data set−set 4. We chose to test GB-neck2 on set 4 because GBneck2 uses a different set of parameters for nucleic acids,150 and our data set 4 is the only one that contains nucleic acid fragments. For the test we employed mbondi3 radii (the suggested default for GB-neck2) instead of Bondi radii. As before, we used the matching set of atomic radii for the reference PB calculations. Our results show that the RMSD deviations and R2 correlations change only slightly when using these different atomic radii other than Bondi (Table 7). Thus, all of our conclusions based on the Bondi radii set remain valid.

difference between the dielectric boundaries used by the GB and PB should further reduce the discrepancy between the respective ΔΔGel values. This appears to be the case: an algorithm that produces a molecular surface very similar to that used by the PB solvers was recently implemented136 for GBNSR6 in AMBER, which reduced the ΔΔGel RMSD error with respect to PB by more than a factor of 2, from 3.9 to 1.43 kcal/mol for data set 1 for which the evaluation was performed.136 One might conjecture that direct use of the perfect effective Born in the GB equation would bring the GB accuracy into an even closer accord with the PB but that may not be the case, at least based on the two structures we tested here. Namely, when the perfect radii were used instead of the R6 radii in the GB eq 2, we found no reduction of the GB versus PB error in ΔΔGel, Table 6. The fact that for GBNSR6 the agreement with the PB

Table 7. Sensitivity of Some GB Models to the Choice of Atomic Radii Seta

Table 6. PB ΔΔGel Values (kcal/mol) Compared to GBNSR6 (Bondi Radii) and the GB Using Perfect Effective Born Radiia PDB ID

PB

GBNSR6 (Bondi radii)

GB (perfect radii)

1BKF 1FKL

−0.76 12.28

−9.19 6.23

−10.72 3.98

data set 3 (protein−protein) ΔΔGel (RMSD) Bondi radii GBSW GBMV1 GBMV2 GBMV3

a

When the latter are used in the GB formalism, the remaining error relative to the PB comes from the GB Green function itself.

19.0 49.1 25.9 16.8

ΔΔGel (R2)

CHARMM radii

Bondi radii

29.6 0.949 43.7 0.946 18.4 0.991 20.5 0.982 data set 4 (RNA-peptide) ΔΔGel (RMSD)

does not necessarily improve with the use of the perfect Born radii suggests that for this GB model a significant component of the remaining error in ΔΔGel likely comes from the approximate GB Green function itself and not from the way the effective Born radii are approximated109,133 within this model. It remains to be seen whether more sophisticated GB Green functions123 can further reduce the GB versus PB gap in ligand binding calculations. Insensitivity to the Input Radii Set. We have explicitly confirmed that our key conclusions are not affected by the use of Bondi intrinsic atomic radii instead of recommended defaults for the models where this choice may matter. Here we tested CHARMM GB models (i.e., GBSW, GBMV1, GBMV2, GBMV3) on our data set 3, using CHARMM atomic radii instead of Bondi radii set. We chose this set of structures for the sensitivity test because this set shows the poorest R2 correlation

GB-neck2

CHARMM radii 0.972 0.937 0.971 0.960

ΔΔGel (R2)

Bondi radii

mbondi3 radii

Bondi radii

mbondi3 radii

65.7

64.7

0.957

0.953

Deviations of ΔΔGel values (kcal/mol) relative to the PB values using the same set of radii for both GB and PB. Included here are root mean square deviation (RMSD) and square of the Pearson correlation coefficient (R2).

a

■

CONCLUSIONS High computational efficiency of the methods used in calculation of binding free energies is important, since these calculations involve large structural ensembles. The GB models offer a good compromise: they offer efficient methods to account for solvent without explicitly modeling individual water molecules. Of course, the accuracy of the models is absolutely 1666


Article


free energies generally perform better in binding affinity calculations. (2) The GB models that lead to poor correlations with PB in electrostatic solvation free energies are least accurate in terms of electrostatic binding free energies. (3) Nevertheless, the same accuracy rank order of GB models in the electrostatic solvation free energies does not necessarily hold in the corresponding electrostatic binding affinities predictions;56,63 especially when the differences in accuracies are relatively small. We stress that a close agreement of a GB model with the numerical PB, while being an important metric of the GB accuracy,124 does not, by itself, guarantee an equally good performance of the model compared to experiment. There are several reasons for the lack of a simple direct correlation here. First, the PB itself is a serious approximation to reality−a mean field, linear response treatment of atomistic electrostatic interactions in the presence of discrete water molecules and mobile ions replaced by a structureless continuum within the PB. Equally importantly, the PB accounts for the electrostatics only, while the total binding free energy includes other important contributions, most notably nonpolar and entropic, the latter being particularly difficult to estimate computationally−these contributions have not been considered here. At the same time, a significant disagreement of a GB model with the PB on a certain class of structures is a serious indication of potential problems in practical calculations. Indeed, many publications based on MM(GB)SA and MM(PB)SA have shown a large variability in their predictive power when compared against experimental binding data.156−160 Based on the very different accuracies of the GB models tested here and elsewhere (e.g., refs 151 and 161−163), we think this variability may stem largely from the very different electrostatic binding free energies predicted by the different GB models. The results presented here can therefore be used as a guide for the MM(GB)SA/MM(PB)SA community in deciding which GB models should be avoided, and which are more likely than others to lead to more accurate results. The much improved level of accuracy in more recent GB methods provides an encouraging basis for more reliable GB simulations and opens a whole new range of possibilities for reaching larger system sizes in ligand binding calculations. Future studies will need to evaluate the most promising GB models against appropriate experimental data. In this sense, “electrostatic” mutants in which a charged amino acid is mutated into a neutral one may be particularly useful, as the electrostatic effects should dominate the corresponding free energy changes to the binding free energy. Finally, the detailed description of parameters of many practical GB flavors provided here may serve as a useful reference guide for those wishing to utilize these models in ligand binding computations.

critical. This study is an attempt to systematically investigate the accuracy of several commonly used GB models implemented in major molecular simulation packages in estimating the electrostatic component of the binding free energy, ΔΔGel. The accuracy was evaluated in terms of agreement with the numerical PB reference across a diverse set of 60 biomolecular complexes. Our findings demonstrate that different GB models lead to very different binding affinity predictions, with R2 correlations varying largely among different complex sets, as well as among different GB models. An assessment of the literature reveals a similar behavior of R2 correlations for other biomolecular complexes with varying energy range and spacing (e.g., refs 151−155). Future studies comparing the GB-predicted solvation free energies directly to experiment will need to evaluate to what extent performance of the promising GB models holds in realistic problems where binding free energies differ by only a few kcal/mol. Among all different complex sets tested in this study, the best agreement with the PB was achieved for the protein−protein set, most likely because protein structures have been conventionally used as a large part of the training sets used in the parametrization of most GB models. The correlation with PB is also good for the small complexes, which are relatively small and not highly charged. However, the RNA-peptide and protein-drug complex sets appeared to be most challenging. The RNA-peptide complexes are highly charged, which can directly scale up the errors in GB Green function. The force field parameters for “nonstandard” drug components seem to be least compatible with the GB parameters, leading to a relatively low correlation with PB for the protein-drug test. In general, we found that newer GB models are in better agreement with PB. Among the CHARMM GB models and considering all biomolecular complexes together, GBMV3 shows the highest correlation with the PB results for ΔΔGel, whereas GBMV2 is the most accurate for ΔGel. Moreover, GBMV1, which is an older analytical model, is less accurate than the other CHARMM GB models when all biomolecular complexes are combined into one data set. Among the GB models available in AMBER, GB-neck2, which has the largest number of fitting parameters, shows a good performance (second to GBNSR6) across different test cases. In contrast to the other GB models tested here, GBNSR6 uses only a single additional (to the PB) fitting parameter, yet it shows the best agreement with reference PB among all of the models tested. The closer agreement of GBNSR6 with PB consistently holds over all four complex sets. Moreover, the accuracy of GBNSR6 ΔΔGel estimates does not change significantly if the perfect (PB-based) Born radii are used instead of the R6 radii calculated within GBNSR6. Therefore, we conjecture that the majority of the remaining error in GBNSR6 estimates of ΔΔGel are due to the approximate GB Green function itself, rather than the way the effective Born radii is calculated. Future work will assess if this superior accuracy of the GBNSR6 model in Amber translates into better binding affinity ranking when performing MM(GB)SA studies of biomolecular systems for which solid experimental data is available. We have also explored how the accuracy of the GB models in predicting ΔGel correlates with their ability to predict the binding affinity ΔΔGel. A summary of our finding with respect to the PB accuracy metric is as follows: (1) GB models that tend to be considerably more accurate in electrostatic solvation

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00886. A README file, one folder for each of the four data sets, and *.dat files describing the names of the files in the folder; Amber topology (*.top) and coordinate (*.crd) and *.pqr files are included in each of the four folders (ZIP) 1667


Article


■

(22) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665− 9678. (23) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. J. Phys. Chem. Lett. 2014, 5, 3863−3871. (24) Söderhjelm, P.; Tribello, G. A.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 5170−5175. (25) Proctor, E. A.; Yin, S.; Tropsha, A.; Dokholyan, N. V. Biophys. J. 2012, 102, 144−151. (26) Gao, K.; Yin, J.; Henriksen, N. M.; Fenley, A. T.; Gilson, M. K. J. Chem. Theory Comput. 2015, 11, 4555−4564. (27) Izadi, S.; Aguilar, B.; Onufriev, A. V. J. Chem. Theory Comput. 2015, 11, 4450−4459. (28) Cramer, C. J.; Truhlar, D. G. Chem. Rev. 1999, 99, 2161−2200. (29) Honig, B.; Nicholls, A. Science 1995, 268, 1144−1149. (30) Beroza, P.; Case, D. A. Methods Enzymol. 1998, 295, 170−189. (31) Madura, J. D.; Davis, M. E.; Gilson, M. K.; Wade, R. C.; Luty, B. A.; McCammon, J. A. Biological Applications of Electrostatic Calculations and Brownian Dynamics Simulations. In Reviews in Computational Chemistry; 1994; Vol. 5, pp 229−267, DOI: 10.1002/ 9780470125823.ch4. (32) Gilson, M. K. Curr. Opin. Struct. Biol. 1995, 5, 216−223. (33) Scarsi, M.; Apostolakis, J.; Caflisch, A. J. Phys. Chem. A 1997, 101, 8098−8106. (34) Luo, R.; David, L.; Gilson, M. K. J. Comput. Chem. 2002, 23, 1244−1253. (35) Simonson, T. Rep. Prog. Phys. 2003, 66, 737−787. (36) Onufriev, A. Continuum Electrostatics Solvent Modeling with the Generalized Born Model. In Modeling Solvent Environments, 1st ed.; Feig, M., Ed.; Wiley: USA, 2010; pp 127−165, DOI: 10.1002/ 9783527629251.ch6. (37) Labute, P. J. Comput. Chem. 2008, 29, 1693−1698. (38) Nicholls, A.; Honig, B. J. Comput. Chem. 1991, 12, 435−445. (39) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (40) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219−10225. (41) Im, W.; Beglov, D.; Roux, B. Comput. Phys. Commun. 1998, 111, 59−75. (42) Koehl, P. Curr. Opin. Struct. Biol. 2006, 16, 142−151. (43) Altman, M. D.; Bardhan, J. P.; White, J. K.; Tidor, B. J. Comput. Chem. 2009, 30, 132−153. (44) Li, B.; Cheng, X.; Zhang, Z. SIAM J. Appl. Math. 2011, 71, 2093−2111. (45) Simonov, N. A.; Mascagni, M.; Fenley, M. O. J. Chem. Phys. 2007, 127, 185105. (46) Boschitsch, A. H.; Fenley, M. O.; Zhou, H.-X. J. Phys. Chem. B 2002, 106, 2741−2754. (47) Boschitsch, A. H.; Fenley, M. O. J. Comput. Chem. 2004, 25, 935−955. (48) Mackoy, T.; Harris, R. C.; Johnson, J.; Mascagni, M.; Fenley, M. O. Commun. Comput. Phys. 2013, 13, 195−206. (49) Boschitsch, A. H.; Fenley, M. O. In Computational Electrostatics for Biological Applications: Geometric and Numerical Approaches to the Description of Electrostatic Interaction Between Macromolecules; Rocchia, W., Spagnuolo, M., Eds.; Springer International Publishing: Switzerland, 2015; pp 73−110, DOI: 10.1007/978-3-319-12211-3_4. (50) Fenley, M. O.; Mascagni, M.; McClain, J.; Silalahi, A. R. J.; Simonov, N. A. J. Chem. Theory Comput. 2010, 6, 300−314. (51) Baker, N. A. Curr. Opin. Struct. Biol. 2005, 15, 137−143. (52) Boschitsch, A. H.; Fenley, M. O. J. Chem. Theory Comput. 2011, 7, 1524−1540. (53) Chen, D.; Chen, Z.; Chen, C.; Geng, W.; Wei, G.-W. J. Comput. Chem. 2011, 32, 756−770. (54) Wang, J.; Tan, C.; Chanco, E.; Luo, R. Phys. Chem. Chem. Phys. 2010, 12, 1194−1202. (55) Sarkar, S.; Witham, S.; Zhang, J.; Zhenirovskyy, M.; Rocchia, W.; Alexov, E. Commun. Comput. Phys. 2013, 13, 269−284. (56) Harris, R. C.; Mackoy, T.; Fenley, M. O. J. Chem. Theory Comput. 2015, 11, 705−712.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (M.O.F.). *E-mail: [email protected] (A.V.O.). ORCID

Marcia O. Fenley: 0000-0003-2325-3150 Alexey V. Onufriev: 0000-0002-4930-6612 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.V.O. thanks the NIH (GM076121) for partial support. All authors thank Michael Feig for many helpful discussions of the GBMV family of models and the associated best practices; Carlos Simmerling for suggestions regarding best radii set for GB-neck2; Wonpil Im for a comment on the GBSW model. M.O.F. thanks Dr. Ilisa Lee for doing a literature search on MMPB(GB)SA studies and Dr. Michael E. Zawrotny for maintaining and helping with computer resources at the Institute of Molecular Biophysics at Florida State University.

■

REFERENCES

(1) Chodera, J. D.; Mobley, D. L.; Shirts, M. R.; Dixon, R. W.; Branson, K.; Pande, V. S. Curr. Opin. Struct. Biol. 2011, 21, 150−160. (2) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; et al. J. Am. Chem. Soc. 2015, 137, 2695−2703. (3) Jorgensen, W. L. Science 2004, 303, 1813−1818. (4) Mobley, D. L.; Dill, K. A. Structure 2009, 17, 489−498. (5) Shirts, M. R.; Mobley, D. L.; Brown, S. P. In Structure Based Drug Design, 1st ed.; Merz, K. M., Ringe, D., Reynolds, C. H., Eds.; Lecture Notes in Computer Science; Cambridge University Press: Cambridge, New York USA, 2010; pp 61−85. (6) Kastritis, P. L.; Bonvin, A. M. J. J. J. R. Soc., Interface 2013, 10, 20120835. (7) Lau, A. Y.; Roux, B. A. Nat. Struct. Mol. Biol. 2011, 18, 283−287. (8) Lin, Y.-L.; Meng, Y.; Jiang, W.; Roux, B. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 1664−1669. (9) Gilson, M. K.; Zhou, H.-X. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21−42. (10) Pollard, T. D. Mol. Biol. Cell 2010, 21, 4061−4067. (11) Kollman, P. Chem. Rev. 1993, 93, 2395−2417. (12) Beveridge, D. L.; DiCapua, F. M. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431−492. (13) Åqvist, J.; Medina, C.; Samuelsson, J.-E. Protein Eng., Des. Sel. 1994, 7, 385−391. (14) Srinivasan, J.; Cheatham, T. E., III.; Cieplak, P.; Kollman, P. A.; Case, D. A. J. Am. Chem. Soc. 1998, 120, 9401−9409. (15) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; Cheatham, T. E., III. Acc. Chem. Res. 2000, 33, 889−897. (16) Gohlke, H.; Kiel, C.; Case, D. A. J. Mol. Biol. 2003, 330, 891− 913. (17) Hu, Y.; Sherborne, B.; Lee, T.-S.; Case, D. A.; York, D. M.; Guo, Z. J. Comput.-Aided Mol. Des. 2016, 30, 533−539. (18) Mobley, D. L.; Dill, K. A.; Chodera, J. D. J. Phys. Chem. B 2008, 112, 938−946. (19) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Shirts, M. R.; Dill, K. A. J. Chem. Theory Comput. 2009, 5, 350−358. (20) Shivakumar, D.; Deng, Y.; Roux, B. J. Chem. Theory Comput. 2009, 5, 919−930. (21) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. 1668


Article

Journal of Chemical Theory and Computation (57) Fenley, M. O.; Harris, R. C.; Mackoy, T.; Boschitsch, A. H. J. Comput. Chem. 2015, 36, 235−243. (58) Hoijtink, G.; De Boer, E.; Van der Meij, P.; Weijland, W. Recl. Trav. Chim. Pays-Bas 1956, 75, 487−503. (59) Tucker, S. C.; Truhlar, D. G. Chem. Phys. Lett. 1989, 157, 164− 170. (60) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127−6129. (61) Bashford, D.; Case, D. A. Annu. Rev. Phys. Chem. 2000, 51, 129− 152. (62) Chen, J.; Brooks, C. L., III; Khandogin, J. Curr. Opin. Struct. Biol. 2008, 18, 140−148. (63) Harris, R. C.; Mackoy, T.; Fenley, M. O. Mol. Based Math. Biol. 2013, 1, 63−74. (64) Feig, M.; Brooks, C. L., III Curr. Opin. Struct. Biol. 2004, 14, 217−224. (65) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Chem. Phys. Lett. 1995, 246, 122−129. (66) Haberthür, U.; Caflisch, A. J. Comput. Chem. 2008, 29, 701−715. (67) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 19824−19839. (68) Schaefer, M.; Karplus, M. J. Phys. Chem. 1996, 100, 1578−1599. (69) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. J. Phys. Chem. A 1997, 101, 3005−3014. (70) Edinger, S. R.; Cortis, C.; Shenkin, P. S.; Friesner, R. A. J. Phys. Chem. B 1997, 101, 1190−1197. (71) Jayaram, B.; Liu, Y.; Beveridge, D. L. J. Chem. Phys. 1998, 109, 1465−1471. (72) Ghosh, A.; Rapp, C. S.; Friesner, R. A. J. Phys. Chem. B 1998, 102, 10983−10990. (73) Lee, M. S.; Salsbury, F. R.; Brooks, C. L., III J. Chem. Phys. 2002, 116, 10606−10614. (74) Felts, A. K.; Harano, Y.; Gallicchio, E.; Levy, R. M. Proteins: Struct., Funct., Genet. 2004, 56, 310−321. (75) Romanov, A. N.; Jabin, S. N.; Martynov, Y. B.; Sulimov, A. V.; Grigoriev, F. V.; Sulimov, V. B. J. Phys. Chem. A 2004, 108, 9323− 9327. (76) Dominy, B. N.; Brooks, C. L., III. J. Phys. Chem. B 1999, 103, 3765−3773. (77) David, L.; Luo, R.; Gilson, M. K. J. Comput. Chem. 2000, 21, 295−309. (78) Tsui, V.; Case, D. A. J. Am. Chem. Soc. 2000, 122, 2489−2498. (79) Calimet, N.; Schaefer, M.; Simonson, T. Proteins: Struct., Funct., Genet. 2001, 45, 144−158. (80) Spassov, V. Z.; Yan, L.; Szalma, S. J. Phys. Chem. B 2002, 106, 8726−8738. (81) Simmerling, C.; Strockbine, B.; Roitberg, A. E. J. Am. Chem. Soc. 2002, 124, 11258−11259. (82) Wang, T.; Wade, R. C. Proteins: Struct., Funct., Genet. 2003, 50, 158−169. (83) Nymeyer, H.; García, A. E. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 13934−13949. (84) Onufriev, A.; Bashford, D.; Case, D. A. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (85) Gallicchio, E.; Levy, R. M. J. Comput. Chem. 2004, 25, 479−499. (86) Lee, M. C.; Duan, Y. Proteins: Struct., Funct., Genet. 2004, 55, 620−634. (87) Im, W.; Lee, M. S.; Brooks, C. L., III J. Comput. Chem. 2003, 24, 1691−1702. (88) Nguyen, H.; Perez, A.; Bermeo, S.; Simmerling, C. J. Chem. Theory Comput. 2015, 11, 3714−3728. (89) Lange, A. W.; Herbert, J. M. J. Chem. Theory Comput. 2012, 8, 1999−2011. (90) Tjong, H.; Zhou, H.-X. J. Phys. Chem. B 2007, 111, 3055−3061. (91) Xu, Z.; Cheng, X.; Yang, H. J. Chem. Phys. 2011, 134, 064107. (92) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. J. Chem. Theory Comput. 2016, 12, 5946−5959. PMID: 27748599.

(93) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26, 1668−1688. (94) Brooks, B. R.; Brooks, C. L., III; MacKerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; et al. J. Comput. Chem. 2009, 30, 1545−1614. (95) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. Bioinformatics 2013, 29, 845−854. (96) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781−1802. (97) Eastman, P.; Pande, V. Comput. Sci. Eng. 2010, 12, 34−39. (98) Qin, S.; Zhou, H.-X. Biopolymers 2007, 86, 112−118. (99) Dong, F.; Zhou, H.-X. Proteins: Struct., Funct., Genet. 2006, 65, 87−102. (100) Knight, J. L.; Brooks, C. L., III. J. Comput. Chem. 2011, 32, 2909−2923. (101) Bartlett, R. J.; Musiał, M. Rev. Mod. Phys. 2007, 79, 291−352. (102) Helgaker, T.; Klopper, W.; Tew, D. P. Mol. Phys. 2008, 106, 2107−2143. (103) Fennell, C. J.; Kehoe, C. W.; Dill, K. A. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 3234−3239. (104) Merz, K. M. J. Chem. Theory Comput. 2010, 6, 1769−1776. (105) Harris, R. C.; Boschitsch, A. H.; Fenley, M. O. J. Chem. Theory Comput. 2013, 9, 3677−3685. (106) Harris, R. C.; Boschitsch, A. H.; Fenley, M. O. J. Chem. Phys. 2014, 140, 075102. (107) Sørensen, J.; Fenley, M. O.; Amaro, R. E. In Computational Electrostatics for Biological Applications: Geometric and Numerical Approaches to the Description of Electrostatic Interaction Between Macromolecules; Rocchia, W., Spagnuolo, M., Eds.; Springer International Publishing: Switzerland, 2015; pp 39−71, DOI: 10.1007/978-3319-12211-3_3. (108) Cramer, C. J.; Truhlar, D. G. Acc. Chem. Res. 2008, 41, 760− 768. (109) Aguilar, B.; Onufriev, A. V. J. Chem. Theory Comput. 2012, 8, 2404−2411. (110) Nguyen, H.; Maier, J.; Huang, H.; Perrone, V.; Simmerling, C. J. Am. Chem. Soc. 2014, 136, 13959−13962. (111) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235−242. (112) Katkova, E. V.; Onufriev, A. V.; Aguilar, B.; Sulimov, V. B. J. Mol. Graphics Modell. 2017, 72, 70−80. (113) Dolinsky, T. J.; Nielsen, J. E.; McCammon, J. A.; Baker, N. A. Nucleic Acids Res. 2004, 32, W665−W667. (114) Dolinsky, T. J.; Czodrowski, P.; Li, H.; Nielsen, J. E.; Jensen, J. H.; Klebe, G.; Baker, N. A. Nucleic Acids Res. 2007, 35, W522−W525. (115) Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A. et al. AMBER 2015; University of California, San Francisco, 2015. (116) Hu, C. Y.; Kokubo, H.; Lynch, G. C.; Bolen, D. W.; Pettitt, B. M. Protein Sci. 2010, 19, 1011−1022. (117) Kokubo, H.; Hu, C. Y.; Pettitt, B. M. J. Am. Chem. Soc. 2011, 133, 1849−1858. (118) Kokubo, H.; Harris, R. C.; Asthagiri, D.; Pettitt, B. M. J. Phys. Chem. B 2013, 117, 16428−16435. (119) Harris, R. C.; Pettitt, B. M. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 14681−14686. (120) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (121) Srinivasan, J.; Trevathan, M. W.; Beroza, P.; Case, D. A. Theor. Chem. Acc. 1999, 101, 426−434. (122) Onufriev, A.; Case, D. A.; Bashford, D. J. Comput. Chem. 2002, 23, 1297−1304. 1669


Article

Journal of Chemical Theory and Computation (123) Onufriev, A. V.; Sigalov, G. J. Chem. Phys. 2011, 134, 164104. (124) Onufriev, A. V.; Izadi, S. Wiley Interdisciplinary Reviews: Computational Molecular Science; e1347-n/a, e1347. (125) Onufriev, A. V.; Aguilar, B. J. Theor. Comput. Chem. 2014, 13, 1440006. (126) Sigalov, G.; Scheffel, P.; Onufriev, A. J. Chem. Phys. 2005, 122, 094511. (127) Hawkins, G.; Cramer, C.; Truhlar, D. J. Phys. Chem. 1996, 100, 19824−19839. (128) Schaefer, M.; Van Vlijmen, H.; Karplus, M. Adv. Protein Chem. 1998, 51, 1−57. (129) Tsui, V.; Case, D. Biopolymers 2000, 56, 275−291. (130) Onufriev, A.; Bashford, D.; Case, D. J. Phys. Chem. B 2000, 104, 3712−3720. (131) Onufriev, A.; Bashford, D.; Case, D. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (132) Nguyen, H.; Roe, D. R.; Simmerling, C. J. Chem. Theory Comput. 2013, 9, 2020−2034. (133) Aguilar, B.; Shadrach, R.; Onufriev, A. V. J. Chem. Theory Comput. 2010, 6, 3613−3630. (134) Sigalov, G.; Fenley, A.; Onufriev, A. J. Chem. Phys. 2006, 124, 124902. (135) Mongan, J.; Simmerling, C.; McCammon, J. A.; Case, D. A.; Onufriev, A. J. Chem. Theory Comput. 2007, 3, 156−169. (136) Forouzesh, N.; Izadi, S.; Onufriev, A. V. J. Chem. Inf. Model. 2017, 57, 2505−2513. PMID: 28786669. (137) Homeyer, N.; Gohlke, H. Mol. Inf. 2012, 31, 114−122. (138) Miller, B. R., III; McGee, T. D., Jr; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. J. Chem. Theory Comput. 2012, 8, 3314− 3321. (139) Homeyer, N.; Gohlke, H. J. Comput. Chem. 2013, 34, 965− 973. (140) Lee, M. S.; Feig, M.; Salsbury, F. R., Jr.; Brooks, C. L., III J. Comput. Chem. 2003, 24, 1348−1356. (141) Im, W.; Feig, M.; Brooks, C. L., III Biophys. J. 2003, 85, 2900− 2918. (142) Tjong, H.; Zhou, H.-X. J. Chem. Theory Comput. 2008, 4, 507− 514. (143) Nina, M.; Beglov, D.; Roux, B. J. Phys. Chem. B 1997, 101, 5239−5248. (144) Nina, M.; Im, W.; Roux, B. Biophys. Chem. 1999, 78, 89−96. (145) Swanson, J. M. J.; Adcock, S. A.; McCammon, J. A. J. Chem. Theory Comput. 2005, 1, 484−493. (146) Nicholls, A.; Mobley, D. L.; Guthrie, J. P.; Chodera, J. D.; Bayly, C. I.; Cooper, M. D.; Pande, V. S. J. Med. Chem. 2008, 51, 769− 779. (147) Chen, J. J. Chem. Theory Comput. 2010, 6, 2790−2803. (148) Rankin, K. N.; Sulea, T.; Purisima, E. O. J. Comput. Chem. 2003, 24, 954−962. (149) Sanner, M. F.; Olson, A. J.; Spehner, J. C. Biopolymers 1996, 38, 305−320. (150) Nguyen, H.; Pérez, A.; Bermeo, S.; Simmerling, C. J. Chem. Theory Comput. 2015, 11, 3714−3728. (151) Vorontsov, I. I.; Miyashita, O. J. Comput. Chem. 2011, 32, 1043−1053. (152) Carra, C.; Saha, J.; Cucinotta, F. A. J. Mol. Model. 2012, 18, 3035−3049. (153) Genheden, S.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Ryde, U. J. Phys. Chem. B 2010, 114, 8505−8516. PMID: 20524650. (154) Delgado-Soler, L.; Pinto, M.; Tanaka-Gil, K.; Rubio-Martinez, J. J. Chem. Inf. Model. 2012, 52, 2107−2118. PMID: 22794663. (155) Rastelli, G.; Rio, A. D.; Degliesposti, G.; Sgobba, M. J. Comput. Chem. 2010, 31, 797−810. (156) Wright, D. W.; Hall, B. A.; Kenway, O. A.; Jha, S.; Coveney, P. V. J. Chem. Theory Comput. 2014, 10, 1228−1241. PMID: 24683369. (157) Mikulskis, P.; Genheden, S.; Rydberg, P.; Sandberg, L.; Olsen, L.; Ryde, U. J. Comput.-Aided Mol. Des. 2012, 26, 527−541. (158) Suenaga, A.; Okimoto, N.; Hirano, Y.; Fukui, K. PLoS One 2012, 7, e42846.

(159) Panel, N.; Sun, Y. J.; Fuentes, E. J.; Simonson, T. Frontiers in Molecular Biosciences 2017, 4, 65. (160) Genheden, S. J. Comput.-Aided Mol. Des. 2011, 25, 1085−1093. (161) Ren, X.; Tortorella, M.; Wang, J.; Wang, C. Sci. Rep. 2017, 7, 14934. (162) Ylilauri, M.; Pentikäinen, O. T. J. Chem. Inf. Model. 2013, 53, 2626−2633. PMID: 23988151. (163) Wittayanarakul, K.; Hannongbua, S.; Feig, M. J. Comput. Chem. 2008, 29, 673−685.

1670


Accuracy Comparison of Generalized Born Models ... - ACS Publications

Recommend Documents