Explicit Solvent Hydration Benchmark for Proteins with Application to

Article pubs.acs.org/JCTC

Explicit Solvent Hydration Benchmark for Proteins with Application to the PBSA Method Piotr Setny* and Anita Dudek Centre of New Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland S Supporting Information *

ABSTRACT: Explicit and implicit solvent models have a proven record of delivering hydration free energies of small, druglike solutes in reasonable agreement with experiment. Hydration of macromolecules, such as proteins, is to a large extent uncharted territory, with few results shedding light on quantitative consistency between different solvent models, let alone their ability to reproduce real water. In this work, based on extensive explicit solvent simulations employing TIP3P and SPC/E water models we analyze hydration free energy changes between fixed conformations of 5 diverse proteins, including large multidomain structures. For the two solvent models we find better agreement in electrostatic rather than nonpolar contributions (RMSE of 2.3 and 2.7 kcal/mol, respectively), even though absolute values of the latter are typically an order of magnitude smaller. We also highlight the importance of finite size corrections to relative protein hydration free energies, which turn out to be rather large, on the order of several kcal/mol, and are necessary for proper interpretation of results obtained under periodic boundary conditions. We further compare gathered data with predictions of the implicit solvent approach based on the Poisson equation and the surface or volume based nonpolar term. We find definitely lesser consistency than between the two explicit models (RMSE between implicit and TIP3 results of 11.3 and 8.4 kcal/mol for electrostatic and nonpolar contributions, respectively). In the process we determine the value of the protein dielectric constant and the geometric model for the dielectric boundary that provide for the best agreement. Finally, we evaluate the usefulness of surface and volume based models of nonpolar contributions to hydration free energy of large biomolecules.

1. INTRODUCTION

well as provide an environment for successful, ab inito protein folding simulations.15 Nevertheless, in many applications explicit solvent simulations are too inconvenient. They require a complex system setup, the use of periodic boundary conditions, and additional simulation time for equilibration. They are also costly, since a better part of computer power is spent just on handling solvent related degrees of freedom. Accordingly, simplified, implicit solvent methods16−20 are often considered as an alternative for tasks such as binding free energy estimations,21,22 protein structure prediction,23−25 pKa calculation,26 or mutational protein analysis.27 They usually draw from macroscopic physical concepts and divide solvent-related effects into nonpolar and electrostatic contributions. The first term accounts for the creation of cavity within the solvent necessary to accommodate a solute, together with solute−solvent dispersion interaction. The second term describes free energy of solvent polarization due to solute partial charges. In the most typical setup used in biomolecular modeling, the nonpolar term is assumed to be proportional to the area of solvent accessible surface formed around the solute. More complex models28−30 involve volume-based contributions, a separate term for dispersion interactions, and curvature

Physics based modeling of biomolecular systems inevitably requires accounting for the effects of aqueous environment.1,2 Quantified by hydration free energies, water contributions to processes such as protein folding, ligand binding, or transitions between different functional states are often of comparable magnitude to direct (intra)solute interaction energies.3−5 Biological phenomena are typically driven by overall free energy changes in the range of several kcal/mol.4,6 In order to draw meaningful conclusions from computational studies, one has to reach at least a similar level of accuracy in predicting the underlying thermodynamic contributions. This requires comparable precision not only in the description of biomolecules themselves but also their interactions with the solvent. Explicit solvent models based on rigid, fixed charge representation of water molecules7,8 seem to offer an acceptable compromise between the accuracy of energy representation and the extent of obtainable sampling. Certainly, they do not provide an ideal description of an aqueous environment.7,8 Typical problems include the reproduction of temperature dependent properties,9,10 balance between water−water and water-solute dispersion interactions,11 or inconsistent results concerning protein dynamics.12 On the other hand, they are capable of delivering satisfactory predictions of hydration free energies for small molecules13 and neutral amino acids,9,14 as © XXXX American Chemical Society

Received: March 9, 2017 Published: May 12, 2017 A

DOI: 10.1021/acs.jctc.7b00247 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

implicit models on parametrization, in particular the definition of the solute−solvent boundary.66 In the current study we aim at obtaining a possibly reliable set of relative hydration free energies for multiple protein conformations. In addition, we separately consider electrostatic and nonpolar components. We use thermodynamic integration (TI)71 based on explicit solvent molecular dynamics (MD) simulations in which a protein is gradually morphed between two, fixed conformations. The obtained hydration free energy differences are suitable for a straightforward comparison with implicit solvent predictions. We consider a set of five diverse proteins ranging from small, peptide-like chingolin to large, multidomain structures. In addition we analyze hydration of neutral amino acids. Taken together, this gives us the possibility to highlight the challenges which implicit models need to face when transitioning from the description of small, generally water exposed solutes to large macromolecules with solventinaccessible cores. A separate investigation of nonpolar and electrostatic contributions offers a unique opportunity to directly observe what is the proportion between the two effects in shaping the protein hydration free energy landscape and how well each of them is reproduced by an implicit solvent model. Our results are based on two identical, yet independent sets of MD simulations, each employing a different water model: TIP3P72 and SPC/E.73 The assessment of their consistency serves to further validate the numerical results for hydration free energy predictions. In addition, it gives a reference level for the interpretation of discrepancies observed for implicit solvent description. In this respect, we compare explicit solvent results with a Poisson-based model for electrostatic contributions and surface area and volume based estimates of nonpolar effects. We also evaluate how their performance depends on the particular setup of the solute−solvent boundary and the protein dielectric constant. We admit that a thorough evaluation of numerous versions, parametrizations, and implementations of different implicit solvent models is far beyond our study. We believe, however, that the obtained set of explicit solvent free energy estimates, along with provided structures, may serve to do so, thus contributing to further methodological developments.

corrections. The electrostatic term is based on the assumption that the solute constitutes a low dielectric volume which contains fixed, discrete atomic charges and is embedded in high dielectric solvent continuum. The energy of solute charges in the solvent reaction field can be then found by solving the Poisson equation. Augmented with the assumption of electric field dependent ionic density, such an approach gives rise to the popular Poisson−Boltzmann model.31 In turn, its simplification, based on phenomenological interpolation between exact solutions for isolated atomic ions proposed by Born,32 is the foundation of the generalized Born model.33 Implicit models allow fast evaluation of hydration free energies. The price for simplification is a number of important limitations, such as invariance to the sign of solute charge,34 insensitivity of the solute−solvent boundary to physicochemical solute properties and its topography,35,36 or problems with the correct representation of nonpolar contributions.37 Nonetheless, constant efforts are being made to optimize the existing models38−41 and develop them further to address the above deficiencies.19,30,42−51 Their predictions of hydration free energies for small, druglike molecules have been compared with numerous explicit solvent results as well as with experimental data,52−57 not only exposing some obvious problems but also giving a realistic overview of what can be expected. In contrast, the performance of implicit solvent models in describing hydration of large biomolecules such as proteins, that is the kind of solutes they are actually designed for, is comparably less well assessed in a quantitative manner. The obvious reason for this fact is the lack of experimental data concerning macromolecular hydration that would be suitable for direct comparison. Little such data can be also found within computational studies, since obtaining total hydration free energies of large proteins based on explicit solvent simulations is technically difficult. Notable attempts to validate implicit solvent models in the context of protein hydration are based on conformational analysis. They either assess the ability of the polypeptide chain to fold into its native conformation in the implicit solvent environment23,25,58,59 or compare important regions of the free energy landscape derived from explicit and implicit solvent simulations.60−62 Such studies lead to somewhat mixed conclusions. On the one hand, a number of small proteins were shown to be able to fold into correct geometries25 or remain stable in native conformations.23,58 On the other hand, the results concerning apparently less stable structures, for which the actual conformation may be more dependent on the interplay between internal and solvation energies, are less consistent.59 Indeed, more direct evaluation of solvation free energies seems to confirm the general lack of quantitative agreement between explicit and implicit solvent models.60−62 Unfortunately, the above approaches do not allow dissecting the reference, explicit solvent estimates into electrostatic and nonpolar contributions. In part, such insights are instead obtained in studies focusing on only nonpolar37,63 or electrostatic effects. The latter can be evaluated by direct calculations of charging free energies for polypeptides and small proteins in fixed conformations,39,49,64−66 electrostatic contributions to ligand binding,67−69 or salt bridge formation.70 The results generally confirm only qualitative rather than quantitative agreement between explicit and implicit solvent description. They also indicate considerable dependence of

2. METHODS 2.1. Protein Structures. As a source of benchmark structures we selected 5 proteins representing diverse folds, with polypeptide chain lengths ranging from 10 to 238 amino acids (Table 1). Four of them served to generate 4 sets, each Table 1. List of Protein Structures Considered in This Study protein

PDB id

n at. (aa)

sec. struct

n conf

chignolin (CH) protein G (PG) λ-repressor (LR) adenylate kinase (AK) lao protein (LP)

1uao 1qkz 1lmb 4ake/1ake 2lao/1lst

138(10) 247(16) 1258(85) 3341(214) 3649(238)

beta beta/helix α-helix mixed mixed

4 4 10 4 4

containing 4 distinct conformations, and one, λ-repressor, a single set with 10 conformations. In each case, the protein structure was extracted from the respective file from the Protein Databank (PDB)74 and parametrized with the Amber-ff-99 force field,75 with default protonation states assigned by the pdb2gmx GROMACS76 tool. Particular conformations were generated based on original, energy minimized PDB structures, B


Article

Journal of Chemical Theory and Computation and their diverse conformers were extracted from explicit solvent simulations (for a detailed description of the procedure and conformations see the Supporting Information). The set included the following: chignolin (CH): based on a solution NMR structure (PDB id: 1uao), with folded (f), misfolded (m), prefolded (p), and unfolded (u), states; protein G (PG): extracted from a PDB structure (PDB id: 1qkz, chain A, residues 45−61) with beta (b), misfolded (m), unfolded (u), and an artificially generated canonical helix (h) state; λ-repressor (LR): based on a crystal structure (PDB id. 1lmb), with 10 conformations, divided into 3 sets that shared one common member corresponding to the native (n) structure. The sets were divided into the following: a small rms set (LRS) with 3 additional structures (s1−s3), an intercluster set (LRI), with 3 additional structures (i1−i3), and an unfolding set (LRU), with 3 additional structures (u1−u3); adenylate kinase (AK): with open (o) and closed (c) conformations (PDB id: 4ake, 1ake, respectively) and semi (s) and semi2 (s2) states; lao binding protein (LP): with open (o) and closed (c) conformations (PDB id: 2lao, 1lst, respectively) and semi (s) and semi2 (s2) states. We emphasize that particular details of functional states of proteins are of lesser importance here − we focus just on hydration free energies for fixed, protein-like solutes with no aim at characterizing the actual role of the aqueous environment on folding or conformational changes. MD Simulations. MD simulations were carried with the GROMACS 5.0.6 package.76 MD runs serving to generate benchmark structures were conducted in the explicit solvent with the TIP3P water model, in NpT conditions, with temperature 300 or 400 K, and a reference pressure of 1 bar. Periodic boundary conditions were used with the dodecahedral simulation box, providing at least a 12 Å space between protein and the box sides at the starting configuration. All bonds were constrained using the LINCS algorithm, and a time step of 2 fs was used. A 10 Å cutoff was used for shifted van der Waals interactions with the standard Gromacs long-range correction, and the particle mesh Ewald summation method was used for electrostatic interactions, with a 10 Å cutoff in real space and 1.2 Å grid spacing for Fourier space. Following 1000 steps of steepest descent minimization, 100 ps of the NVT run was carried out with harmonic restraints on heavy atoms (with 2.4 kcal/mol/Å2 force constant), and with NpT, an unconstrained production run was conducted with length depending on the structure of interest. 2.2. Free Energy Calculations. Thermodynamic Integration. All free energy differences were estimated between two fixed protein conformations in the explicit solvent environment with a constant number of particles. For a system consisting of a protein fixed in configuration X and water with degrees of freedom Y, free energy reads F(X) = UP(X) − kBT ln

∫ dYe−[UPW(X, Y) + UW(Y)]/ kBT ∫ dYe−UW(Y)/ kBT

ΔF(XA → XB) = UP(XB) − UP(XA) +

∫0

1

dλ

∂UPW(λ) ∂λ

= ΔUP + ΔFH λ

(2)

with the latter term corresponding to the desired change in hydration free energy, ΔFH, and the λ parameter controlling the transition from XA (λ = 0) to XB (λ = 1), and UPW(λ) ≡ UPW(X(λ), Y). Note that as we consider only hydration contributions to free energy, for brevity we put ΔFH ≡ ΔF in the following. A transition path chosen for TI corresponded to gradual, linear transformation (morphing) of the initial structure into the final one. Typically, it was divided into 64 or 128 windows, depending on protein size and RMSD between the two end points. Protein conformation in each such window remained fixed, and only water degrees of freedom were propagated during MD. As only protein−water interaction energies were of interest, nonbonded interactions between protein atoms were switched off in order to avoid spurious energy readings due to protein atom overlaps encountered for nonphysical, intermediate configurations. Solvent Treatment. In order to avoid trapping water molecules in narrow cavities between constrained protein atoms in the course of simulations, an additional step was executed during the preparation stage for each window. Based on initial solvent configuration, either obtained from the GROMACS solvation procedure (the first window) or taken from the last frame of the preceding window, all water molecules closer to any protein atom than 3.5 Å were shifted into a random location within the bulk solvent region (no closer than 4 Å from protein atoms). Possible steric clashes between water molecules were then resolved using 500 steps of steepest descent minimization with soft-core potential turned on for water molecules (as implemented in the standard free energy perturbation GROMACS protocol), followed by up to 5 × 105 minimization steps with standard nonbonded potentials. Such a procedure assured that no contributions from internal hydrated protein cavities would be included in the calculated hydration free energy differences. In our view, handling the energetic effects of such buried solvent requires specialized methods both on the side of explicit solvent simulations77,78 and simplified solvent models.79,80 Their inclusion here would further complicate the analysis of bulk protein hydration. Volume Adjustment. An ensemble of choice to perform hydration free energy calculations corresponds to NpT conditions. Unfortunately, the use of constraints, necessary to fix protein geometries in subsequent windows, precludes correct virial evaluation, making constant pressure simulations unfeasible. In turn, constant volume simulations, involving conformational transitions of large proteins, lead to significant changes in solvent accessible volume. They result in a nonnegligible and hard to estimate pressure−volume term (note that solvent pressure against which protein volume expands/ contracts in a small, fixed simulation box may easily change by thousands of bars due to low water compressibility). To circumvent this issue, we simulated each TI window in a volume adjusted to provide constant water density, of value corresponding to ρ0, the density of the appropriate solvent model (see below) at ambient conditions (T = 300 K, p = 1 bar). Although under such conditions pressure−volume work is not included in the TI formula (eq 2), its contribution remains very small due to the fact that volume changes occur at low,

(1)

where UP and UW denote protein and water potential energies, respectively, UPW denotes the protein water interaction, and kBT is the Boltzmann constant times temperature. A standard derivation of thermodynamic integration formula71 leads to the free energy difference between structures A and B C


Article

Journal of Chemical Theory and Computation ambient pressure. Based on the observed volume changes we estimate it to be at most on the order of 0.1 kcal/mol (the LP case with the largest and smallest observed box size of 89.5 and 89.2 Å, respectively). In order to maintain constant solvent density, prior to starting the first window of TI simulations we performed 1 ns of the NpT run for the initial protein structure with harmonic positional restraints of 2.4 kcal/mol/Å2 on Cα atoms. Its last 500 ps served to evaluate an average system volume, which was then used for the initial equilibration run of 500 ps at NVT conditions, with the now constrained protein. The apparent bulk solvent density, ρ, was then evaluated in a spherical shell of 5 Å thickness and outer diameter corresponding to the simulation box size. Next, the volume of the box was adjusted by the amount δV = N(1/ρ0 − 1/ρ), where N was the number of water molecules in the system. Finally, a new 500 ps of NVT equilibration was performed, and the procedure was iterated until ρ was within 0.1% from the target density, ρ0 (typically, a single iteration is enough). For each subsequent TI window such an equilibration scheme was repeated with an initial volume inherited from the preceding window instead of from the NpT run. Model specific water densities at ambient conditions were evaluated using 1 ns NpT runs, preceded by 100 ps NVT thermalization, with starting configurations corresponding to the (60 Å)3 cube of pure solvent, using the same MD settings as for TI calculations. We evaluated average densities based on the last 500 ps, obtaining ρ0,TIP3 = 0.0329 Å−3 and ρ0,SPCE = 0.0334 Å−3 for TIP3P and SPC/E water models, respectively, in agreement with values reported in the literature.8 Simulation Protocol. For each protein we performed 4 sets of TI simulations: with and without partial charges on protein atoms, each independently for TIP3P and SPC/E water models. Together they served to evaluate total changes in hydration free energy and their nonpolar and electrostatic contributions, as well as to assess the agreement between the two solvent models. All TI simulations were conducted in the NVT setup (see above), using periodic boundary conditions, with cubic simulation boxes. For each protein system the box size provided at least a 12 Å margin of solvent. Water molecules were added using the gmx solvate tool, and no counterions were included. Protein geometries were constrained during MD runs. Water molecules were kept rigid using the SETTLE algorithm. A time step of 2 fs was used. The solvent temperature was kept at 300 K using a Nosé−Hoover thermostat. Nonbonded interactions were treated as described above for free MD simulations. Each TI window was preceded by sets of 500 ps equilibration runs (see above), and production runs were conducted for 2 ns (charged proteins) or 1 ns (proteins with no charges). Simulation frames were saved every 1 ps for analysis. 2.3. Analysis Methods. Thermodynamic Integration. Changes in hydration free energy were evaluated according to the TI formula (eq 2). Transition (linear morphing) from the initial to the final protein conformation was divided into equidistant windows, enumerated with the λ parameter that changed from 0 to 1. Trajectories for each λ value were obtained as specified above. Derivatives of UPW with respect to λ were then evaluated numerically in a postprocessing step, during which the trajectories were reprocessed two times with protein conformation X(λ) changed by δλ = ± 0.002

∂UPW(λ) ∂λ

= λ

1 nλ

nλ

∑ i=1

UPW(λ + δλ)i − UPW(λ − δλ)i 2δλ (3)

where nλ was the number of frames recorded for a given λ value, and UPW was evaluated as a difference between the potential energy of the entire system and protein only, under identical system (periodic) geometry and nonbonded interactions treatment. We checked that the results were numerically stable for δλ ranging from 0.005 to 0.0005. The final integration over λ was performed using the trapezoidal rule. In order to assess random errors of the estimated hydration free energy changes, we used a bootstrap procedure. ΔF values were calculated 105 times based on resampled simulation frames, and the standard deviation, σ, of the resulting distribution was taken as uncertainty of the result. Thermodynamic Cycles. Free energy changes between a set of conformations can be summed up into a closed cycle. Due to the fact that free energy is a function of state, such a sum should be equal to 0. One can benefit from such an observation in three ways. First, having evaluated free energy changes from independent simulations, one can validate the results by checking if their sum is indeed 0 within the limit of estimated uncertainties. For instance, for 4 conformations ABCD, one can consider 4 unique combinations of three-segment cycles (Figure 1) with sums

Figure 1. Equivalence between thermodynamic cycles (up) for three (left) and four (right) structures and systems of points connected by springs and allowed to relax in one dimension (down).

SXYZ = ∑iΔFi, where i ∈ {XY, YZ, ZX}. Their respective uncertainties are σXYZ = (∑iσ2i )1/2, where σi is the individual segment error obtained with the bootstrap procedure described above. Second, in order to make the results suitable for comparison with implicit solvent calculations (which provide absolute hydration free energies and, hence, free energy differences that trivially sum up to 0 for closed cycles), one can adjust them by seeking a solution that provides for the correctly closed cycle subject to minimal perturbation of the original values. The problem is analogous to finding potential energy minimum of a set of points (protein conformations) connected by springs with reference lengths {f 0i } (original free energy differences from simulations, ΔFi) and force constants {ki} (inverse variances of original free energy results, σ−2 i ), allowed to relax in one dimension. For sets of three and four such points (Figure 1), the problem is described by the following equations, in which energy of the system is equivalent to Mean D


Article


solvent model for nonperiodic (ΔFN,i) and periodic (ΔFP,i) systems, of the same geometries and box sizes as in MD simulations:

Square Error (MSE) between the original and regularized free energies MSEABC =

MSEABCD =

3

1 2

∑ ki(fi

1 2

∑ ki(fi

− f i0 )2 + λ 0(f1 − f2 − f3 )

i=1

P,e N → P,i ΔFAN,e − ΔΔFBN → P,i = ΔFAP,e→ B + FSC → B = ΔFA → B + ΔΔFA (6)

(4)

The idea is illustrated by the thermodynamic cycle depicted in Figure 2. The full cycle

6

− f i0 )2 + λ1(f3 − f1 − f2 )

i=1

N,e → i ΔFAN,e − ΔFAN,i + ΔFAP,i + ΔFAP,i → e + ΔFAP,e→ B → B = ΔFA

+ λ 2(f6 − f4 − f1 ) + λ3(f4 − f2 − f5 ) + λ4(f6 − f3 − f5 )

+ ΔFBP,e → i − ΔFBP,i + ΔFBN,i + ΔFBN,i → e

(5)

where λi are Lagrange multipliers. Minimization of MSE can be performed analytically by solving a set of linear equations. It results in a set of equilibrium lengths {f opt i }, which are equivalent to regularized free energy differences that fulfill the requirement of the vanishing sum over any cycle. They also provide for minimal deviation from the original results, weighted by their accuracies. Furthermore, the resultant spring constant between any two points can be interpreted as the inverse variance of the respective, regularized free energy change and, hence, serves as its error estimate (see the Supporting Information). The third benefit of the above approach is that each of the final, regularized free energy differences incorporates information from all simulations and, hence, has decreased uncertainty with respect to the original result. It is worth noting, that the regularization procedure effectively decreases the number of independent free energy differences for n protein conformations from n(n − 1)/2 (ΔF for all possible pairs, evaluated in independent simulations) to n − 1 (ΔF between any conformation of choice and all the others). Accordingly, having four conformations for each of the considered proteins, we performed six independent TI calculations of pairwise free energy differences and applied the above regularization procedure. In the case of the λrepressor, three such sets, sharing one common, native structure, were treated in the same manner. Free Energy Contributions. We separately estimated hydration free energy differences for proteins with and without partial charges. They correspond to the total hydration free energy differences, ΔF, and nonpolar contributions, ΔFnp, respectively. Their difference gives the magnitude of electrostatic contribution, ΔFel = ΔF − ΔFnp, which is directly comparable with the results of PB or GB calculations. Uncertainties of electrostatic contributions obtained with TI simulations were evaluated by error propagation from the total and nonpolar contributions. Finite Size Corrections (FSC). The use of periodic boundary conditions in simulations together with the treatment of Coulomb interactions based on the Ewald summation method leads to finite size effects.81−85 The first-order effect is associated with the change of net solute charge and is not relevant for the considered simulations. Nevertheless, the large difference in dipole and higher order moments of two protein conformations may alone lead to large discrepancies. The desired hydration free energy change from state A to B in the nonperiodic, explicit solvent system, ΔFN,e A→B, can be derived from the TI result obtained under periodic boundary P,e conditions, ΔFA→B , by applying the FSC. Here, we use a numerical scheme for FSC.84 For each of the two end states we evaluate the differences in electrostatic contributions: ΔΔFN→P,i = ΔFP,i − ΔFN,i, obtained with the implicit, Poisson-based

(7)

Figure 2. Thermodynamic cycle used for the definition of finite size corrections to free energies based on explicit solvent simulations under PBC.

can be simplified to eq 6 by the assumption that differences between explicit and implicit solvent descriptions arising due to short ranged, solute−solvent boundary effects cancel out between periodic and nonperiodic treatments,84 resulting in ΔFN,e→i + ΔFP,i→e = 0. Indeed, FSC evaluated for test protein systems proved to be insensitive to the solvent boundary definition used for solving the Poisson equation (molecular vs solvent accessible surface). In order to test the correctness of the above FSC we conducted additional simulations to evaluate a) hydration free energy changes for two pairs of monovalent atomic ions (with the same and opposite charges of 1e magnitude) as a function of their separation, using two simulation boxes of drastically different sizes and b) all ΔF values for the PG case, using smaller and larger than standard simulation boxes. In both cases, the application of FSC proved necessary to reach the agreement of hydration free energies obtained under different PBC. The detailed results are provided in the Supporting Information. The magnitude of FSC is obviously sensitive to the length of the simulation box. Due to different densities of TIP3P and SPC/E bulk water, box sizes in simulation of identical protein structures were slightly different (up to 0.5 Å). The resulting differences in FSC for several test cases were below 0.1 kcal/ mol, and hence, we used box sizes from TIP3P simulations for all corrections. Solvent Distribution. Spatial water density maps (number density of water oxygen atoms) were prepared using the VMD86 volmap tool, with (0.5 Å)3 voxels. They served to analyze the number and location of high density surfacial hydration spots, surface areas, and volumes of solvent excluded E


Article

Journal of Chemical Theory and Computation Table 2. Explicit Solvent Results and Their Uncertainties (Subscript) in kcal/mol A→B CH f, u (T) f, u (S) f, m (T) f, m (S) f, p (T) f, p (S) u, m (T) u, m (S) u, p (T) u, p (S) m, p (T) m, p (S) PG b, h (T) b, h (S) b, m (T) b, m (S) b, u (T) b, u (S) h, m (T) h, m (S) h, u (T) h, u (S) m, u (T) m, u (S) LRU n, u1 (T) n, u1 (S) n, u2 (T) n, u2 (S) n, u3 (T) n, u3 (S) u1, u2 (T) u1, u2 (S) u1, u3 (T) u1, u3 (S) u2, u3 (T) u2, u3 (S) LRI n, i1 (T) n, i1 (S) n, i2 (T) n, i2 (S) n, i3 (T) n, i3 (S) i1, i2 (T) i1, i2 (S) i1, i3 (T)

ΔF̃P,e

FSC

ΔFN,e

ΔFnp

ΔFel

−115.20.6 −116.40.6 −15.40.3 −15.80.3 −40.90.2 −41.00.2 99.90.6 98.60.6 74.20.6 74.60.6 −25.10.2 −25.70.3

−1.6 −1.6 1.1 1.1 0.6 0.6 2.7 2.7 2.2 2.2 −0.4 −0.4

−116.80.4 −117.10.4 −14.50.2 −14.80.2 −40.10.2 −40.50.2 102.30.4 102.40.4 76.70.4 76.60.4 −25.60.2 −25.80.2

0.80.3 0.90.3 −1.10.2 −1.20.2 1.70.1 1.70.1 −2.00.3 −2.10.3 0.90.3 0.80.3 2.80.1 2.90.2

−117.60.5 −118.00.5 −13.30.3 −13.60.3 −41.80.2 −42.30.2 104.30.5 104.40.5 75.80.5 75.80.5 −28.40.2 −28.70.3

−124.11.1 −122.81.2 −124.30.8 −124.60.8 −197.11.2 −197.21.2 −0.41.1 −1.41.1 −72.90.9 −72.81.0 −71.81.0 −73.11.1

−1.4 −1.4 −1.0 −1.0 2.8 2.8 0.4 0.4 4.1 4.1 3.8 3.8

−125.40.7 −124.80.8 −125.50.7 −125.50.8 −194.00.7 −194.20.8 −0.20.7 −0.70.7 −68.60.7 −69.40.8 −68.50.7 −68.70.7

3.30.5 2.90.6 3.50.5 2.80.6 6.60.6 7.10.6 0.30.5 −0.10.5 3.40.5 4.20.5 3.10.5 4.30.5

−128.61.0 −127.71.0 −129.10.8 −128.20.8 −200.60.9 −201.21.0 −0.40.9 −0.60.9 −72.00.9 −73.60.9 −71.60.9 −73.00.9

157.91.4 155.31.4 80.21.0 82.81.0 30.11.5 34.51.5 −76.21.1 −77.81.1 −116.13.3 −115.61.8 −45.91.5 −42.21.5

4.7 4.7 1.3 1.3 −2.3 −2.3 −3.4 −3.4 −7.0 −7.0 −3.6 −3.6

161.00.9 161.20.9 81.20.8 82.00.8 30.51.1 35.31.0 −79.80.9 −79.20.9 −130.51.3 −125.91.1 −50.71.7 −46.71.0

2.80.7 3.10.7 −2.60.6 −3.50.7 17.10.8 20.90.8 −5.40.6 −6.50.6 14.30.8 17.80.9 19.70.8 24.40.8

158.21.2 158.21.2 83.81.0 85.51.2 13.41.3 14.41.3 −74.41.0 −72.71.1 −144.81.5 −143.81.4 −70.41.3 −71.01.3

103.60.8 105.30.8 −83.00.9 −82.10.9 53.20.5 53.70.5 −187.10.9 −182.91.0 −50.10.8

0.3 0.3 −0.5 −0.5 0.1 0.1 −0.8 −0.8 −0.2

104.40.5 104.90.5 −82.90.6 −81.40.6 53.10.4 53.70.4 −186.90.6 −186.40.6 −50.10.6

1.90.4 1.90.4 −1.80.4 −1.70.4 0.50.3 0.50.3 −3.70.4 −3.60.4 −1.40.4

102.20.7 103.00.7 −81.00.7 −79.80.7 52.60.5 53.20.5 −183.20.7 −182.80.8 −49.60.7

a

A→B

ΔF̃P,e

FSC

ΔFN,e

ΔFnp

ΔFel

i1, i3 (S) i2, i3 (T) i2, i3 (S) LRS n, s1 (T) n, s1 (S) n, s2 (T) n, s2 (S) n, s3 (T) n, s3 (S) s1, s2 (T) s1, s2 (S) s1, s3 (T) s1, s3 (S) s2, s3 (T) s2, s3 (S) AK c, s (T) c, s (S) c, o (T) c, o (S) c, s2 (T) c, s2 (S) s, o (T) s, o (S) s, s2 (T) s, s2 (S) o, s2 (T) o, s2 (S) LP c, o (T) c, o (S) c, s (T) c, s (S) c, s2 (T) c, s2 (S) o, s (T) o, s (S) o, s2 (T) o, s2 (S) s, s2 (T) s, s2 (S)

−52.20.8 133.50.9 135.70.8

−0.2 0.7 0.7

−51.20.6 135.90.6 135.20.6

−1.40.4 −2.30.4 −2.20.4

−49.80.7 133.60.7 133.00.7

34.20.7 36.40.7 42.40.7 42.40.7 −64.80.7 −67.10.7 8.10.6 7.30.6 −99.40.4 −100.50.4 −106.20.6 −107.60.6

1.2 1.2 1.4 1.4 0.9 0.9 0.3 0.3 −0.3 −0.3 −0.6 −0.6

35.00.5 37.00.5 43.00.5 44.30.5 −64.80.5 −63.80.8 8.00.4 7.30.4 −99.70.3 −100.80.3 −107.70.4 −108.10.4

−3.30.3 −4.60.4 −0.70.3 −1.30.3 −1.70.3 −2.60.3 2.60.2 3.30.3 1.60.2 2.00.2 −1.00.2 −1.30.3

38.30.5 41.60.6 43.70.5 45.60.6 −63.00.5 −61.20.6 5.40.5 4.10.5 −101.30.4 −102.80.4 −106.70.5 −106.80.5

−224.81.6 −221.41.6 −43.31.8 −40.51.9 −139.92.5 −135.21.0 174.91.3 176.21.0 93.61.0 90.41.0 −86.71.3 −88.81.4

−1.8 −1.8 −1.8 −1.8 0.8 0.8 0.0 0.0 2.6 2.6 2.6 2.6

−226.61.2 −224.20.9 −48.81.2 −46.71.0 −132.01.3 −132.80.8 177.71.0 177.40.8 94.50.8 91.40.8 −83.21.0 −86.10.9

9.10.5 17.10.6 9.70.5 18.00.6 5.70.5 11.60.5 0.60.5 0.90.5 −3.40.5 −5.50.5 −4.00.5 −6.30.5

−235.71.3 −241.31.1 −58.51.3 −64.71.1 −137.71.4 −144.41.0 177.21.1 176.61.0 97.90.9 96.90.9 −79.21.1 −79.81.1

188.51.7 197.41.3 159.11.3 161.50.9 25.31.3 32.10.9 −35.01.4 −34.40.9 −173.72.7 −172.30.9 −129.90.7 −131.10.5

−3.1 −3.1 −4.3 −4.3 −0.8 −0.8 −1.2 −1.2 2.3 2.3 3.5 3.5

188.51.1 196.20.8 152.30.9 158.20.7 25.20.9 30.10.6 −36.21.0 −38.20.7 −163.31.1 −166.30.6 −127.10.7 −128.10.4

−0.90.4 2.80.4 34.20.4 39.50.4 37.30.3 42.80.3 35.00.4 36.70.4 38.20.3 39.90.3 3.10.3 3.30.3

189.31.1 193.60.9 118.11.0 118.70.8 −12.10.9 −12.70.7 −71.21.1 −74.90.8 −201.41.1 −206.30.7 −130.20.7 −131.40.5

ΔF̃ P,e − raw changes in hydration free energy estimated by TI under PBC, FSC − finite size corrections (eq 6), ΔFN,e − regularized (see Methods) change in hydration free energy for nonperiodic, explicit solvent system, ΔFnp, ΔFel − nonpolar and electrostatic contributions, respectively. (T), (S) − results for TIP3P and SPC/E solvent models, respectively.

a

regions for the considered protein structures. The latter were defined as regions, r, with solvent density ρ(r) < 0.5ρ0. Their measurements were carried out using Chimera87 volumetric tools. 2.4. Amino Acids Hydration. Hydration free energies of 16 neutral amino acids, capped with N-terminal acetyl and Cterminal amide groups, were calculated using the standard free energy perturbation protocol available in GROMACS. Energyminimized conformations were positioned in the center of cubic simulations boxes providing at least a 12 Å solvent

margin. Nonpolar contributions to hydration free energies were evaluated during NpT runs at 1 bar pressure and temperature of 300 K, using stepwise annihilation of positionally restrained solutes (harmonic potential with force constant of 2.4 kcal/ mol/Å2 was applied to all heavy atoms) with no partial atomic charges (21 windows, each of 500 ps, preceded by 100 ps equilibration), with the use of soft-core potentials as implemented in GROMACS. Electrostatic contributions were evaluated by the stepwise charging procedure (10 windows of 1 ns each, preceded by 100 ps equilibration) of solutes in F


Article

Journal of Chemical Theory and Computation constrained geometries at NVT conditions with T = 300 K, and volume was taken as an average from 500 ps of NpT simulation of restrained solutes at 1 bar. All settings for MD runs were the same as described above. 2.5. Implicit Solvent Calculations. Solutions of the Poisson equation (note that no counterions were present in the reference, explicit solvent TI calculations) were obtained with the CHARMM88 C39b2 PBEQ module and the APBS program.89 The former was used for the calculations of finite size corrections, owing to its capability to handle periodic boundary conditions. The latter, owing to its computational efficiency, was used for the calculation of implicit solvent hydration free energies with various parameter sets considered for comparison with TI results. The grid used for solving the Poisson equation for finite size corrections covered the size of the respective simulation box with the protein structure placed exactly as in the case of MD simulations. Grid resolution was 0.2 Å. The protein dielectric constant, ϵp, was set to 1. Hydration free energy under periodic and nonperiodic conditions was calculated as a difference between protein charging free energies in vacuum (ϵv = 1) and ”water”, for which the dielectric constant was ϵw = 78.5. As a solute−solvent boundary we considered van der Waals protein surface, as well as the molecular surface obtained with solvent probe of 1 Å radius. We found no difference, however, in the magnitude of finite size corrections obtained with either of the two (the results differed at most by 0.1 kcal/mol). The atomic radii for protein were taken as corresponding to the minimum of the LJ potential from respective AMBER99 force field atom types. We note that no neutralizing background charge density was considered for calculations under periodic conditions as suggested by Rocklin et al.,84 since in our case the correction was evaluated for two conformations with equal net charges, leading to the cancellation of net-charge contributions. Calculations of absolute protein hydration free energies with APBS were performed with a focusing technique. The finest grid, with resolution 0.2 Å, extended 8 Å past solute atoms, and a coarse grid, with the same number of points, extended by a further 10 Å. For such settings we observed the orientationdependent spread of hydration free energy values in the range of up to 1 kcal/mol for the largest structures (AK, LP). The default Debye−Hückel boundary potential with κ = 0 was used. Hydration free energy was calculated as a difference between protein charging in dielectric medium of ϵw = 78.5 and vacuum. For the sake of comparison with TI results, a set of protein dielectric constants ranging from 0.1 (see discussion for rationale behind using such unphysical values) to 3.0 with the step of 0.1 was used. The dielectric boundary was defined as a molecular surface, and a set of solvent probe radii ranging from 0 to 1.8 Å was considered. For protein atomic radii we considered both the values from the AMBER force field (minimum of the atomic LJ potential) referred to as ”Amberbased radii” in the following and from the Bondi set.90 For the estimation of nonpolar contributions, surface areas and volumes of the solvent-excluded region were evaluated with the MSMS program91 based on solvent-accessible and molecular surface definitions using the range of solvent probe radii from 0 to 1.8 Å.

order of 100 kcal/mol (Table 2). Particular values appear not to depend much on protein size nor RMSD between conformational states. While the larger size or the larger magnitude of conformational change tends to lead to more pronounced hydration free energy changes, by no means does this constitute a rule. For instance, the small PG structure exhibits ΔF variations of similar scale to more than 10 times larger LP. Also, hydration free energy differences between structurally similar members of the LR small-rms set are in some cases larger than for definitely more diversified members of the LR unfolding set. Most likely this results from the fact that major contributions to hydration free energy originate from solvent-exposed charged and polar amino acid side chains. Small variations in their instantaneous spatial arrangement are responsible for the large noise in ΔF values, that dominates over systematic contributions due to global conformational changes. TI simulations resulted in random errors (see the Methods section) on the order of 1−2 kcal/mol, which typically led to relative uncertainties not exceeding 0.02. Similar relative errors are achieved in calculations of absolute hydration free energies for small, druglike solutes based on alchemical molecular dynamics simulations.57 Insight into systematic errors, for instance related to solvent interactions with unphysical solute geometries occurring at intermediate TI windows, can be gained by the analysis of thermodynamic cycle closures. Having 6 independent ΔF values for each set of 4 protein conformations one can consider 4 nonredundant, threesegment cycles (see the Methods section). The obtained sums, SXYZ, of ΔF values turned out to be consistent with 0 within three standard deviations, σXYZ, in all cases for uncharged proteins and in most cases for charged proteins (Figure 3). Deviations were observed for TIP3P

Figure 3. Closures of thermodynamic cycles for nonpolar contributions (left) and total hydration free energy changes (right). |∑iΔFi| with the expected value of 0 is plotted as a function of cycle lengths (∑i|ΔFi|). For each protein case, average values over 4 independent cycles are presented. Error bars correspond to an average of respective standard deviations (see Methods). Empty symbols denote cases in which at least one cycle was found with the sum greater than 3 standard deviations. Note the change in scale between plots on the left and the right side.

simulations of charged LRu, AK, and LP sets, where in each case discrepancies were found in 2 out of 4 cycles. Their investigation revealed that such pairs of unclosed cycles always shared a common segment, which was then assumed to suffer from systematic problems. For the sake of the subsequent regularization procedure, uncertainties of such unsettled segments were increased (hence their weights with respect to

3. RESULTS 3.1. Explicit Solvent Simulations. Convergence. Explicit solvent simulations revealed differences between hydration free energies of protein conformations that were typically on the G


Article


to the overall change in hydration free energy, and the subscript indicates the results obtained with a particular water model; electrostatic and nonpolar components are denoted in the text as ΔFel and ΔFnp, respectively). Hydration free energy changes evaluated with the two different solvent models are fairly similar, with root-mean-square error (RMSE) for all 42 collected ΔF values of 2.4 kcal/mol (Figure 4), average

the remaining 5 unchanged segments decreased) such that for 1 cycles containing those segments σXYZ = 3 SXYZ . Interestingly, the analysis of cycle closures depicted in Figure 3 shows that the convergence of nonpolar free energies is somewhat better for the TIP3P water model than for the SPC/ E model, while the opposite is true for simulations involving charged proteins. One plausible explanation for the TIP3P advantage in nonpolar simulations is its markedly faster diffusion rate,7 which may provide for better overall sampling. At the same time, due to its slightly smaller size (by comparison of Lennard-Jones parameters), TIP3P water may be more prone to being trapped by rough surface elements of intermediate, nonphysical rigid structures encountered in our simulations. We note that since our simulation setup involved nonstandard and rigid solute geometries, the above observations concerning the results for charged proteins should not be automatically extended to simulations of flexible, well behaved solutes. Finite Size Corrections. Most current MD simulations are conducted in periodic boundary conditions (PBC), employing lattice sum methods for the evaluation of electrostatic interactions. This naturally involves artifacts related to artificial system periodicity. Their two dominant components arise (1) due to solute interactions with its own copies and neutralizing background charge density and (2) due to undersolvation resulting from reduced solvent polarization.81−84 Both those contributions have opposite signs and cancel out to a large extent when the overall electrostatic energy of the system is considered.82,84 Importantly, however, electrostatic contributions to hydration free energies depend only on the second component, hence do not benefit from such cancellation. It necessitates the use of FSC for the sake of proper comparison of hydration free energy results based on different PBC setups or nonperiodic conditions. Since the leading contribution to undersolvation occurs in processes involving the change in net solute charge, hydration free energies of neutral, small molecules obtained under PBC are not much affected. Indeed, FSC estimated for neutral amino acids considered here are typically on the order of 0.01 kcal/ mol. The situation is different in the case of protein conformational changes. Here, large variations in dipole and higher order moments of charge distribution lead to substantial finite size effects on the order of several kcal/mol, even though the net charge remains intact (Table 2). Typically, the omission of FSC leads to underestimation of relative hydration free energies for states with higher dipole moments in simulations with PBC. Given the variety of protein structures it is not possible to assume that such states generally correspond to less folded, extended conformations, as has been observed in some studies.62,82 For instance, such is indeed the case for LR conformation u3, with minimal content of secondary structure and high dipole moment. Its solvation free energy with respect to the native, compact state is indeed underestimated (not favorable enough) by 7 kcal/mol. The opposite is observed, however, for PG, where hydration of the unfolded state relative to the compact, yet polarized helix state is overestimated by 4 kcal/mol (dipole moments of all considered protein conformations are given in the Supporting Information). TIP3P vs SPC/E. For the sake of all comparisons in the following, the direction of hydration free energy changes (Table 2) was always taken such that ΔFTIP3 ≥ 0 (as above, ΔF refers

Figure 4. Root mean square errors between the estimates of free energies and their contributions based on TIP3P and SPC/E simulations (Table 2). Star symbols indicate cases in which two sets of values are statistically equivalent with the p-value for the χ2 test greater than 0.1. Agreement upon removal of some outlier cases is indicated by brackets − see text for details.

difference (⟨ΔFTIP3 − ΔFSPCE⟩) of −0.5 kcal/mol, and the slope of linear fit to ΔFTIP3 vs ΔFSPCE of 1.00 (see the Supporting Information for detailed plots). Discrepancies for particular cases turn out to be not proportional to the magnitude of ΔF values (correlation between |ΔFTIP3 − ΔFSPCE| and ΔFTIP3, measured by the square of Pearson’s coefficient yields R2 = 0.02) nor to the RMSD between the two conformations (R2 = 0.01) but rather scale with the system size (correlation between RMSE for each protein and the number of its atoms is R2 = 0.78). As evidenced by the problems with the closure of thermodynamic cycles for the two largest proteins (Figure 3), at least in part, this may be attributed to lower accuracy of TI in the case of conformational changes that involve a larger number of atoms. The RMSE of nonpolar contributions (2.7 kcal/mol) is comparable to the one obtained for the total free energy changes, in spite of the fact that the magnitudes of the latter are on average more than 10 times larger than the magnitudes of the former (⟨ΔFnp⟩ = 7.0 kcal/mol, while ⟨ΔF⟩ = 93.3 kcal/ mol). As before, discrepancies between the two solvent models apparently increase with the solute size. Notably, however, in this case satisfactory convergence of TI calculations across closed cycles was achieved for all considered proteins (Figure 3). Together with the average difference between ΔFnp for TIP3P and SPC/E models of −1.4 kcal/mol, and the slope of linear fit equal to 1.12, this may indicate that the SPC/E model provides for systematically larger nonpolar free energy contributions. Indeed, such an observation is further corroborated by independent estimates of its surface tension being ∼1.2 larger than of that of TIP3P water.8 Strikingly, differences in nonpolar contributions to hydration free energies do not seem to propagate to electrostatic components, even though ΔFel are indirectly calculated as ΔFel = ΔF − ΔFnp and hence are expected to accumulate discrepancies from the two independent terms. The effect is particularly evident for larger structures (LR, AK, LP), for which inconsistencies of electrostatic contributions are markedly smaller than of nonpolar contributions. Furthermore, an average difference between TIP3P and SPC/E-derived ΔFel estimates of −0.9 kcal/mol, and the slope of linear fit of 1.01, indicate a better overall agreement than in the case of the H


Article


Figure 5. RMSEs for the Poisson-based electrostatic hydration component with respect to TIP3P ΔFel, for neutral amino acids (aa), individual protein structures, and all protein structures (ALL). Numbers in the upper-right corners denote individually optimal (rs, ϵ) pairs and corresponding RMSEs in kcal/mol. Dashed lines follow minimal RMSEs for the entire rs range for neutral amino acids and serve as a reference for the remaining plots. White dots and neighboring numbers denote the position of global RMSE minimum, (rs, ϵ) = (1.2, 1.1), and respective RMSEs for each protein structure.

agreement was observed only for ΔFel, for two out of three sets (LRU and LRI). In the case of electrostatic contributions for AK and LP proteins, it was possible to identify individual outlier conformations (AK:closed, LP:open), as all free energy changes involving those conformations were consistently different between TIP3P and SPC/E models by the same offset (see Table 2). Their removal and evaluation of χ2 based on ΔFel values for the 3 remaining conformations and two degrees of freedom resulted in p-values of 0.61 and 0.29 for AK and LR, respectively. The above analysis again indicates that the best agreement between TIP3P and SPC/E hydration was achieved for electrostatic contributions. We note, however, that uncertainties for indirectly calculated ΔFel values are a bit larger than for independent ΔF and ΔFnp, and hence the conditions to fulfill the χ2 test are less strict. The division of hydration free energy between electrostatic and nonpolar contributions is a somewhat arbitrary one, since it depends on the existence of an unphysical, uncharged state. Nevertheless, it is a useful concept for understanding and modeling of hydration effects. As indicated by our results, the consistency in reproduction of electrostatic effects between two different explicit solvent models is quite remarkable, given the wide range of observed ΔFel values. Apparently, the requirements to model solvent polarization due to solute charges are less demanding than to properly describe the process of cavity creation and dispersion interactions. While subtle factors such as water structure and hydrogen bonding next to the solute surface are just a background for the strong electrostatic response, they are the essence for much weaker, nonpolar effects. Accordingly, one can speculate that the nonpolar effects are more dependent on the solute size, surface topography, and the solvent model, as evidenced by the worse overall agreement between ΔFnp for TIP3P and SPC/E.

nonpolar component. Discrepancies here also scale with the solute size; however, a very good agreement of even large ΔFel, obtained for smaller proteins (PG, LR cases), suggests that systematic problems in the reproduction of bulk electrostatics are not significant. Instead, we believe that inconsistencies arise due to accumulation of model-dependent differences in the interaction of water particles with local peculiarities of the protein structure. In each of the AK and LP conformations we observed a few isolated spots of extremely high local solvent density (ρ > 6ρ0, Supporting Information), which typically occurred in the vicinity of charged surface groups forming tight pockets, capable of sequestering individual water molecules. The probability of finding such nontrivial hydration sites obviously increases with increasing size and complexity of protein structure, possibly explaining the size dependence of ΔFel differences. We admit, however, that precise investigation of the nature of underlying, model dependent, solvent effects, as well as their relevance for flexible solute structures, is beyond our current study, and hence the conclusion regarding the source of discrepancies in ΔFel is in part speculative. We further assessed the agreement between free energy estimates obtained with the two solvent models by considering the χ2 test. For each set of 4 protein conformations we calculated the statistics 6

χ2 =

∑ i=1

(ΔFTIP3, i − ΔFSPCE, i)2 2 2 σTIP3, i + σSPC/E, i

(8)

based on all 6 available free energy differences and evaluated pvalues assuming 3 degrees of freedom (note that 3 degrees of freedom are lost due to the regularization procedure). The results for the two solvent models were deemed equivalent, if the p-value >0.1. Such equivalence was observed for CH and PG sets for all three free energy components. For the LR, the I


Article

Journal of Chemical Theory and Computation 3.2. Implicit Solvent Results. Electrostatic Contributions. For a possibly direct comparison between explicit and implicit solvent models we conducted MD simulations without counterions. We checked in a set of additional simulations of the PG structure with 3 neutralizing counterions that the similarity of MD results for ΔF with and without counterions was well within the estimated uncertainties (as per Table 2). Most likely salt-dependent effects would become non-negligible for higher ionic concentrations. The use of the pure aqueous solvent, however, allows us to consider the implicit model based just on the Poisson equation, thus avoiding complications in the interpretation of hydration free energy results due to ionic effects. Implicit solvent estimates obtained with different combinations of rs and the protein dielectric constant (ϵ) were compared with ΔFel derived from TIP3P simulations. In Figure 5 we show the results obtained for the molecular surface (MS)based solute−solvent boundary,92 with protein atomic radii taken from the Bondi set (the results for Amber-based radii are qualitatively the same − see the Supporting Information for the respective plot). The lowest joint RMSE of 11.3 kcal/mol for all 42 protein ΔFel values (Figure 5, ALL) is achieved for rs = 1.2 Å, and ϵ = 1.1. For the results based on the atomic radii derived from the Amber force field we get the minimum joint RMSE of 12.5 kcal/mol for rs = 1.0, and ϵ = 1.2. Interestingly, other studies also report the best consistency of results obtained for rs ∼ 1 Å, rather than for the most often used, default value of 1.4 Å.64,84 It is evident that no single (rs, ϵ) pair provides universally optimal results; however, some trends in the behavior of electrostatic solutions can be identified. Generally, smaller structures (amino acids, CH) and cases with no substantial rearrangement of protein core (LRS and LRI) show lesser sensitivity to rs values. An extreme case here is the amino acids set, for which RMSE ∼ 1 kcal/mol was obtained for the entire considered rs range, subject to minor changes in the ϵ value. We believe that the range of usable rs values for large, complex structures is narrowed down by two different mechanisms. In the limit of small rs, the occurrence of false regions of high dielectric within the solvent inaccessible protein core makes for a much overestimated overall electrostatic effect (note in Figure 5 that as rs decreases, the optimal ϵ values increase, so that the electrostatic part of hydration free energy is scaled down).66 Even more importantly, such artificial solvent distribution is not properly modulated by conformational changes, leading to improper scaling of contributions from particular atomic partial charges due to the spurious solvent accessibility pattern (this is also evidenced by low correlation of MS and MD-based solvent surface areas for small rs, see the Supporting Information). Obviously, the importance of the above effects increases with the fraction of truly solventinaccessible atoms in the solute structure and the degree to which solvent accessibility is modified as a result of conformational change. In turn, in the limit of large rs, discrepancies may be related to the inability of MS to adequately sample the details of solute surface. The existence of surface peculiarities, such as deep narrow pockets, is particularly likely between neighboring, yet topologically distant, solvent exposed solute fragments. Though such configurations are not limited to large proteins, the likelihood of their occurrence scales up with structure size and complexity. An isolated example of this type of situation is a likely reason for unproportionally large RMSE observed for PG given its small structure size. Implicit solvent errors arising for rs

> 0.75 Å (Figure 5, PG) are caused by underestimated hydration free energy for the misfolded state: all ΔFel involving this state indicate its undersolvation by ∼16 kcal/mol with respect to explicit solvent results (note that TI3P and SPC/E results are in good agreement for PG). Upon exclusion of this outlier conformation, the RMSE based on the remaining 3 conformations, evaluated at rs = 1.2 Å and ϵ = 1.1, drops from 11.6 to 3.8 kcal/mol. The analysis of explicit solvent distribution around the misfolded structure reveals a highdensity hydration site (ρmax ≥ 6ρ0) located in a pocket, that is tight enough to remain buried under the MS-based dielectric boundary for solvent probes with rs > 0.5 Å. Its ”hydration” occurring for smaller rs results in the misfolded state being no longer an outlier and in better overall agreement of implicit solvent results with ΔFel from MD. Nonpolar Contributions. For the evaluation of nonpolar contributions we considered solvent accessible surface (SAS)93 based on Bondi and Amber atomic radii, with rs of 1.0 and 0.7 Å, respectively, as such combinations provided the most consistent correlations between the areas of constructed surfaces and explicit solvent boundaries (see the Supporting Information). We assumed that the nonpolar term scales linearly with SAS area (SASA) and sought for the proportionality constant, γ, that results in the minimal RMSE with respect to explicit solvent ΔFnp. This lead to γ of 0.014 and 0.016 kcal/mol/Å2 for the Bondi and Amber radii, respectively. Although the particular value of such microscopic surface tension depends on the assumed physical nature of the nonpolar term (cavity contributions only vs cavity plus dispersion), the kind of explicit solvent model, or definition of the solute−solvent boundary, the obtained values are comparable to those considered in other studies under similar assumptions: 0.01729and 0.01562 kcal/mol/Å2. Since we observed no significant differences between the predictions based on Bondi and Amber radii, in the following we focus on the former ones. The table summarizing all obtained values is in the Supporting Information. The agreement of SASA-based predictions with TIP3P estimates of the nonpolar contributions is poor (Figure 6, left).

Figure 6. Comparison of nonpolar terms based on SASA: Fnp ∼ γA (left) and SEV: Fnp ∼ pV (right) using solute surface area (A) and volume (V) derived for Bondi radii, with ΔFnp obtained for the TIP3P water model.

While correlations between the magnitude of surface area changes and explicit solvent ΔFnp are observed for individual proteins, the proportionality constant is apparently not transferable. The universally optimal value mentioned above leads to the overall RMSE of 8.6 kcal/mol and squared Pearson’s correlation coefficient, R2 = 0.40. Since the selected rs provides for good correlation of SASA and explicit solvent boundary areas (⟨R2⟩ = 0.93, see the Supporting Information), such results indicate that the poor performance of the SASAJ


Article


in the case of differences between explicit solvent models, discrepancies between TIP3P and implicit model predictions do not correlate with the magnitudes of individual free energy differences (R2 = 0.01) nor the RMSD of conformational changes (R2 = 0.01). Strong positive correlation (R2 = 0.81) is instead observed between RMSE for each structure and its size (the number of amino acids). An apparent anomaly here is the PG structure which is small but has an exceptionally large RMSE, dominated by contributions from the electrostatic component. Most likely it is a consequence of inadequate representation of the dielectric boundary for the misfolded conformation, as discussed above. Unfortunately, the progression of RMSE as the function of structure size, at least as observed based on our data, is rather fast. While reasonable accuracy is achieved for the smallest CH structure (RMSECH = 2.4 kcal/mol for ⟨|ΔF|⟩CH = 63 kcal/ mol), RMSE in the range of 10 kcal/mol is already reached for the midsize LR structure (Figure 7, B). Such an error is comparable to the energetic effect of strong protein−ligand binding (nanomolar range) or protein folding free energies.4 Discrepancies between implicit and explicit solvent models for the considered data set are on average dominated by electrostatic contributions (Figure 7, B). It should be kept in mind, however, that |ΔFel| are much larger compared to |ΔFnp|. It translates to average relative errors on the order of 0.1 and 1.0 for the electrostatic and nonpolar terms, respectively. The observation of such large relative uncertainties of the SASA based nonpolar term indicates that in many cases it may not contribute useful information. Indeed, if hydration free energy differences are evaluated for small structures (CH) or protein conformations sharing a similar core structure (LRs, LRi), the electrostatic term alone results in similar RMSE (differences up to 0.2 kcal/mol) with respect to the total ΔF, as the sum of two implicit solvent components. For cases with ΔFnp ≳ 5 kcal/mol (they constitute ∼25% of the considered cases, representing large conformational changes in PG or LRu sets, or transitions between conformations of the two largest proteins) the SASAbased nonpolar component indeed shifts the results in a qualitatively correct manner; however, the overall low accuracy precludes quantitative interpretation of the results. As discussed above, the exchange of the surface-based nonpolar term into the volume-based term hugely improves the results for LP, for which the electrostatic component becomes the limiting factor in such a case. Still, however, no systematic improvement is observed for other protein sets, leaving the quality of predictions for total free energy changes at a practically unchanged level.

based nonpolar term is primarily due to the lack of generally preserved proportionality between surface areas and ΔFnp, rather than flaws in SAS construction. Interestingly, seemingly a much better agreement with MD results was obtained for the solvent excluded volume (SEV)based nonpolar term. Here, considering volumes enclosed by SAS generated with rs = 0.7 Å (such a probe radius provided for the best overall correlation of SEV with protein volumes determined from explicit solvent boundary − see the Supporting Information) we got RMSE of 4.0 kcal/mol, and R2 = 0.88 for the proportionality constant p = 0.029 kcal/mol/ Å3 (Figure 6, right). Though the potential usefulness of the volume-based nonpolar term, with an almost identical proportionality constant, has been already noted in the context of small protein structures,63 conclusions regarding its favorable performance should not be too hastily drawn based on our data. Good agreement of SEV-based estimates with MD-based ΔFnp, in particular relatively high R2, is due to four data points representing the results for LP protein, for which ΔFnp > 30 kcal/mol. Their removal from the data set results in a drop of R2 to 0.43 and places the SEV-based model at the same, relatively poor level as the SASA-based model. Total Free Energy Changes. For the comparison with ΔFTIP3 we focused on the combination of Poisson-based electrostatic and SASA-based nonpolar terms − the most popularly used setup. We used the parameters providing for the globally lowest RMSE, as determined above for Bondi atomic radii (Figure 7, A). The results including volume-based

Figure 7. A: Comparison of total changes in hydration free energy (ΔF) for TIP3P simulations and Poisson + SASA implicit solvent. B: RMSE between TIP3P free energies and Poisson + SASA (S) or Poisson + SEV (V) models for individual protein structures and joint RMSE for all 42 conformational changes (ALL). For comparison, shaded areas of the bar plots correspond to the RMSE between TIP3P and SPC/E solvent models.

4. CONCLUSIONS Using explicit solvent simulations we calculated conformationdependent hydration free energy changes and their nonpolar and electrostatic components for a set of 5 diverse proteins. The structures ranged from small, peptide-like chignolin to large, multidomain globular proteins, thus covering the spectrum extending from compounds with generally solvent exposed atoms to macromolecules with large, solvent inaccessible cores. In order to increase the consistency of the results and make them suitable for direct comparison with implicit solvent estimates, we structured the considered conformational transitions into sets of closed cycles. Benefiting from the requirement of zero net free energy change over each cycle, we regularized the original data using our introduced error relaxation procedure. We also addressed the problem of

nonpolar contributions are presented for comparison in Figure 7, B, but they are qualitatively similar. Also, no important differences were found for implicit solvent estimates based on protein atomic radii taken from the Amber force field (see the Supporting Information). As can be seen in Figure 7, A, implicit and explicit solvent results correlate very well (R2 = 0.96); however, due to the wide range of the energy scale spanned by ΔF values, the overall RMSE remains at the level of 12.4 kcal/mol, roughly 5 times larger than that between TIP3P and SPC/E results. Similarly as K


Article

Journal of Chemical Theory and Computation finite size corrections to hydration free energies, required in simulations that use periodic boundary conditions and lattice based electrostatics. We found that they can be on the order of several kcal/mol in the case of macromolecular solutes, even if only conformational change is considered and the net charge remains constant. We note that finite size effects will be much larger in the case of electrostatic contributions to absolute hydration free energies. As seen from the hydration perspective, proteins are distinct from typical small, druglike solutes not just due to their size, but perhaps more importantly, because they are highly flexible, contain multiple solvent exposed charged groups, and form complex solvent accessible surfaces. The most extreme manifestation of flexibility is the transition from unfolded to folded state. Typically it is strongly opposed by hydration free energy (note strong hydration preference to unfolded states for CH and PG, Table 2) as during folding multiple hydrophilic groups are withdrawn from the solvent and form intramolecular interactions. Somewhat less obviously, even small movements of solvent exposed, charged side chains are able to produce rather large variations in hydration free energy (tens of kcal/ mol) between seemingly invariable conformational states (e.g., LRS set). This implies that hydration free energy estimates for macromolecules based on a single conformation may give, to a large extent, spurious results. What makes protein hydration particularly complicated, however, is the complexity of their surface, which forms multiple narrow pockets and crevasses. Such irregularities are capable of sequestering individual water molecules, which gives rise to hydration free energy contributions that are dependent on the molecular nature of the solvent, and become sensitive to subtle details of a particular solvent model. As evidenced by our results, explicit solvent hydration represented by TIP3P and SPC/E water models is consistent when it comes to the reproduction of even large free energy changes for less complex proteins (CH, PG, LR), irrespectively to the magnitude of conformational change. In the case of larger structures (AK, LP), discrepancies start to arise, in spite of the similar range of ΔF. In part they may be attributed to lower accuracy of TI calculations that involve a larger number of transformed atoms. On the other hand, we observed systematic differences in nonpolar contributions, with ΔFnp for the SPC/E model being typically larger than for TIP3P. Given that nonpolar effects scale with solute size, such a systematic effect, not easily recognizable for small molecules typically studied to date, may noticeably influence model-dependent differences in hydration of larger systems. Somewhat surprisingly, indirectly calculated electrostatic contributions turned out to be more consistent than nonpolar terms. They showed very good agreement of even large ΔFel obtained for small to medium proteins (CH, PG, LR) and smaller RMSE than ΔFnp for the two larger structures. We interpret those results as indicative of little systematic difference between the two solvent models in terms of bulk electrostatics. We ascribe the discrepancies observed for larger structures to the accumulation of small differences arising due to water interaction with local peculiarities of protein surface. We note, however, that further studies are required to pinpoint the actual source of model-dependent irregularities in hydration in such cases. The assembled data set gives the unique possibility of direct, quantitative comparison of explicit and implicit solvent predictions of hydration free energy changes for macro-

molecular solutes. In this respect we have found that RMSE between explicit (TIP3P) and implicit (Poisson equation plus surface area based nonpolar term) models is roughly 5 times larger than between both considered explicit models. In most cases it is enough to rank protein conformations according to hydration free energy consistently with explicit solvent reference; however, it does not allow meaningful estimates of measurable solvent contributions. Moreover, it is worth stressing that we considered only relative hydration free energies, and much larger inaccuracies are expected for absolute values, although fortunately those are not needed for most practical purposes. Reasoning based on implicit solvent predictions for large structures is further complicated by large sensitivity of the results to the parameters of the model. Concerning electrostatic contributions, the observation of the optimal internal dielectric constant being very close to 1 might have been expected, given that such value was also used in our explicit solvent simulations. Major difficulty lies in the fact that there is no a priori general rule for defining the solute−solvent boundary. Focusing exclusively on molecular surfaces, we found that there is only a narrow range of solvent probe radii, centered around rs = 1 Å which leads to possibly optimal results. Apparently such size of the solvent probe provides a reasonable balance between the amount of spurious regions of high dielectric inside the solventinaccessible core and good representation of the details of the actually solvated surface. Furthermore, this range is dependent on the set of protein atomic radii and cannot be easily deduced based on results for generally solvent exposed small solutes, such as amino acids, because their hydration bears little sensitivity to solvent probe radius. The problem may be further complicated by the need to implicitly account for the asymmetry in hydration of opposite charges, for instance by small, electric field dependent adjustments of the surface.34 Predictions for the nonpolar term, either based on solvent accessible surface or solvent excluded volume, did not give reliable agreement with explicit solvent results. Although our determined optimal values for surface tension and pressure coefficients were similar to those reported earlier in the literature, the problem lies in the lack of general correlation between the magnitude of nonpolar effects observed in explicit solvent simulations with solute surface area or volume. Unfortunately, unlike small solutes for which nonpolar contributions are insignificant, justifying their crude, practically qualitative modeling, large protein structures apparently induce nonpolar effects reaching tens of kcal/mol. Their proper estimation becomes equally important to modeling of electrostatic effects, if one aims at implicit solvent approaches that would deliver usable predictions. Finally, there remains a question whether the results of explicit solvent simulations can be regarded as a meaningful reference point. Given years spent on optimizing and finetuning their parameters, it seems that currently available rigid, nonpolarizable explicit solvent models, while not perfect, provide for coherent description of biomolecular hydration at the level compatible with contemporary fixed-charge atomistic force fields.94 From a theoretical standpoint, they account for effects that are typically not built into implicit approaches, such as hydrogen bonding, dielectric saturation, electrostriction,36 or charge asymmetry.34 Improving implicit solvent approaches to a similar level of fidelity would be a challenging, yet rewarding task, and hopefully our quantitative insights into macromolecular hydration may prove helpful. It seems, however, L


Article


(14) Shirts, M. R.; Pande, V. S. Solvation free energies of amino acid side chain analogs for common molecular mechanics water models. J. Chem. Phys. 2005, 122, 134508. (15) Piana, S.; Klepeis, J. L.; Shaw, D. E. Assessing the accuracy of physical models used in protein-folding simulations: Quantitative evidence from long molecular dynamics simulations. Curr. Opin. Struct. Biol. 2014, 24, 98−105. (16) Roux, B.; Simonson, T. Implicit solvent models. Biophys. Chem. 1999, 78, 1−20. (17) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161−2200. (18) Feig, M.; Brooks, C. L. Recent advances in the development and application of implicit solvent models in biomolecule simulations. Curr. Opin. Struct. Biol. 2004, 14, 217−224. (19) Baker, N. A. Improving implicit solvent simulations: a Poissoncentric view. Curr. Opin. Struct. Biol. 2005, 15, 137−143. (20) Chen, J.; Brooks, C. L.; Khandogin, J. Recent advances in implicit solvent-based methods for biomolecular simulations. Curr. Opin. Struct. Biol. 2008, 18, 140−148. (21) Honig, B.; Nicholls, A. Classical electrostatics in biology and chemistry. Science 1995, 268, 1144−1149. (22) 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. Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889−897. (23) Onufriev, A. V.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-Scale Conformational Changes with a Modified Generalized Born Model. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (24) Chen, J.; Brooks, C. L., III Implicit modeling of nonpolar solvation for simulating protein folding and conformational transitions. Phys. Chem. Chem. Phys. 2008, 10, 471−481. (25) Nguyen, H.; Maier, J.; Huang, H.; Perrone, V.; Simmerling, C. L. Folding simulations for proteins with diverse topologies are accessible in days with a physics-based force field and implicit solvent. J. Am. Chem. Soc. 2014, 136, 13959−13962. (26) Bashford, D.; Karplus, M. pKa’s of ionizable groups in proteins: atomic detail from a continuum electrostatic model. Biochemistry 1990, 29, 10219−10225. (27) Massova, I.; Kollman, P. A. Computational alanine scanning to probe protein-protein interactions: A novel approach to evaluate binding free energies. J. Am. Chem. Soc. 1999, 121, 8133−8143. (28) Gallicchio, E.; Levy, R. M. AGBNP: an analytic implicit solvent model suitable for molecular dynamics simulations and high-resolution modeling. J. Comput. Chem. 2004, 25, 479−499. (29) Wagoner, J. A.; Baker, N. A. Assessing implicit models for nonpolar mean solvation forces: the importance of dispersion and volume terms. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 8331−8336. (30) Dzubiella, J.; Swanson, J. M. J.; McCammon, J. A. Coupling nonpolar and polar solvation free energies in implicit solvent models. J. Chem. Phys. 2006, 124, 084905. (31) Sharp, K. A.; Honig, B. Calculating total electrostatic energies with the nonlinear Poisson-Boltzmann equation. J. Phys. Chem. 1990, 94, 7684−7692. (32) Born, M. Volume and heat of hydration of ions. Eur. Phys. J. A 1920, 1, 45−48. (33) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129. (34) Mobley, D. L.; Barber, A. E.; Fennell, C. J.; Dill, K. A. Charge asymmetries in hydration of polar solutes. J. Phys. Chem. B 2008, 112, 2405−2414. (35) Wang, L.; Berne, B. J.; Friesner, R. A. Ligand binding to proteinbinding pockets with wet and dry regions. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 1326−1330.

that optimizing implicit models any further, toward even better reproduction of experimental data, would first require improvements in the physical description on the solute side.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00247. Details concerning the regularization procedure, methodology validation for finite size corrections, descriptor values for considered protein conformations, plots of explicit solvent results, detailed numerical values for implicit solvent results, information concerning the available benchmark data set (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Piotr Setny: 0000-0002-0769-0943 Funding

This work was supported by the Foundation for Polish Science grant Homing-Plus 2012-6/13, EMBO IG 3051/2015, and NCN Sonata UMO-2013/11/D/ST4/02821. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Levy, Y.; Onuchic, J. N. Water mediation in protein folding and molecular recognition. Annu. Rev. Biophys. Biomol. Struct. 2006, 35, 389−415. (2) Ball, P. Water as an active constituent in cell biology. Chem. Rev. 2008, 108, 74−108. (3) Lazaridis, T.; Karplus, M. Thermodynamics of protein folding: a microscopic view. Biophys. Chem. 2002, 100, 367−395. (4) Baldwin, R. L. Energetics of Protein Folding. J. Mol. Biol. 2007, 371, 283−301. (5) Baron, R.; Setny, P.; McCammon, J. A. Water in cavity-ligand recognition. J. Am. Chem. Soc. 2010, 132, 12091−12097. (6) Makhatadze, G. I.; Privalov, P. L. Energetics of Protein Structure. Adv. Protein Chem. 1995, 47, 307−425. (7) Guillot, B. A reappraisal of what we have learnt during three decades of computer simulations on water. J. Mol. Liq. 2002, 101, 219−260. (8) Vega, C.; Abascal, J. L. F. Simulating water with rigid nonpolarizable models: a general perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663−19688. (9) Hess, B.; van der Vegt, N. F. A. Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models. J. Phys. Chem. B 2006, 110, 17616−26. (10) Ashbaugh, H. S.; Collett, N. J.; Hatch, H. W.; Staton, J. A. Assessing the thermodynamic signatures of hydrophobic hydration for several common water models. J. Chem. Phys. 2010, 132, 124504. (11) Best, R. B.; Zheng, W.; Mittal, J. Balanced Protein - Water Interactions Improve Properties of Disordered Proteins and NonSpecific Protein Association. J. Chem. Theory Comput. 2014, 10, 5113− 5124. (12) Florova, P.; Sklenovsky, P.; Banas, P.; Otyepka, M. Explicit Water Models Affect the Specific Solvation and Dynamics of Unfolded Peptides While the Conformational Behavior and Flexibility of Folded Peptides Remain Intact. J. Chem. Theory Comput. 2010, 6, 3569−3579. (13) Mobley, D. L.; Guthrie, J. P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. J. Comput.-Aided Mol. Des. 2014, 28, 711−720. M


Article


(57) Mobley, D. L.; Wymer, K. L.; Lim, N. M.; Guthrie, J. P. Blind prediction of solvation free energies from the SAMPL4 challenge. J. Comput.-Aided Mol. Des. 2014, 28, 135−150. (58) Lwin, T. Z.; Zhou, R.; Luo, R. Is Poisson-Boltzmann theory insufficient for protein folding simulations? J. Chem. Phys. 2006, 124, 034902. (59) Maffucci, I.; Contini, A. An Updated Test of AMBER Force Fields and Implicit Solvent Models in Predicting the Secondary Structure of Helical, β-Hairpin, and Intrinsically Disordered Peptides. J. Chem. Theory Comput. 2016, 12, 714−727. (60) Zhou, R.; Berne, B. J. Can a continuum solvent model reproduce the free energy landscape of a beta -hairpin folding in water? Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12777−12782. (61) Roe, D. R.; Okur, A.; Wickstrom, L.; Hornak, V.; Simmerling, C. L. Secondary structure bias in generalized born solvent models: Comparison of conformational ensembles and free energy of solvent polarization from explicit and implicit solvation. J. Phys. Chem. B 2007, 111, 1846−1857. (62) Cumberworth, A.; Bui, J. M.; Gsponer, J. Free energies of solvation in the context of protein folding: Implications for implicit and explicit solvent models. J. Comput. Chem. 2016, 37, 629−640. (63) Lee, M. S.; Olson, M. A. Comparison of volume and surface area nonpolar solvation free energy terms for implicit solvent simulations. J. Chem. Phys. 2013, 139, 044119. (64) Lee, M. S.; Olson, M. A. Evaluation of Poisson solvation models using a hybrid explicit/implicit solvent method. J. Phys. Chem. B 2005, 109, 5223−5236. (65) Aguilar, B. H.; Shadrach, R.; Onufriev, A. V. Reducing the secondary structure bias in the generalized Born model via R6 effective radii. J. Chem. Theory Comput. 2010, 6, 3613−3630. (66) Onufriev, A. V.; Aguilar, B. H. Accuracy of continuum electrostatic calculations based on three common dielectric boundary definitions. J. Theor. Comput. Chem. 2014, 13, 1440006. (67) Zhang, L. Y. U.; Gallicchio, E.; Friesner, R. A.; Levy, R. M. Solvent Models for Protein - Ligand Binding: Comparison of Implicit Solvent Poisson and Surface Generalized Born Models with Explicit Solvent Simulations. J. Comput. Chem. 2001, 22, 591−607. (68) Pearlman, D. A. Evaluating the molecular mechanics poissonboltzmann surface area free energy method using a congeneric series of ligands to p38 MAP kinase. J. Med. Chem. 2005, 48, 7796−7807. (69) Izadi, S.; Aguilar, B. H.; Onufriev, A. V. Protein-Ligand Electrostatic Binding Free Energies from Explicit and Implicit Solvation. J. Chem. Theory Comput. 2015, 11, 4450−4459. (70) Salari, R.; Chong, L. T. Desolvation Costs of Salt Bridges across Protein Binding Interfaces: Similarities and Differences between Implicit and Explicit Solvent Models. J. Phys. Chem. Lett. 2010, 1, 2844−2848. (71) Straatsma, T. P.; McCammon, J. A. Multiconfiguration thermodynamic integration. J. Chem. Phys. 1991, 95, 1175−1188. (72) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (73) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269−6271. (74) Berman, H. M. The Protein Data Bank. Nucleic Acids Res. 2000, 28, 235−242. (75) Wang, J.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 2000, 21, 1049−1074. (76) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1−2, 19−25. (77) Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The statistical-thermodynamic basis for computation of binding affinities: a critical review. Biophys. J. 1997, 72, 1047−1069.

(36) Muddana, H. S.; Sapra, N. V.; Fenley, A. T.; Gilson, M. K. The electrostatic response of water to neutral polar solutes: implications for continuum solvent modeling. J. Chem. Phys. 2013, 138, 224504. (37) Levy, R. M.; Zhang, L. Y.; Gallicchio, E.; Felts, A. K. On the nonpolar hydration free energy of proteins: surface area and continuum solvent models for the solute-solvent interaction energy. J. Am. Chem. Soc. 2003, 125, 9523−9530. (38) Banavali, N. K.; Roux, B. Atomic radii for continuum electrostatics calculations on nucleic acids. J. Phys. Chem. B 2002, 106, 11026−11035. (39) Swanson, J. M. J.; Adcock, S. A.; McCammon, J. A. Optimized Radii for Poisson - Boltzmann Calculations with the AMBER Force Field. J. Chem. Theory Comput. 2005, 1, 484−493. (40) Mongan, J.; Simmerling, C. L.; McCammon, J. A.; Case, D. A.; Onufriev, A. V. Generalized Born model with a simple, robust molecular volume correction. J. Chem. Theory Comput. 2007, 3, 156− 169. (41) Nguyen, H.; Roe, D. R.; Simmerling, C. L. Improved generalized born solvent model parameters for protein simulations. J. Chem. Theory Comput. 2013, 9, 2020−2034. (42) Yu, Z.; Jacobson, M. P.; Josovitz, J.; Rapp, C. S.; Friesner, R. a. First-Shell Solvation of Ion Pairs: Correction of Systematic Errors in Implicit Solvent Models. J. Phys. Chem. B 2004, 108, 6643−6654. (43) Gavryushov, S. Dielectric saturation of the ion hydration shell and interaction between two double helices of DNA in mono- and multivalent electrolyte solutions: Foundations of the epsilon-modified poisson-boltzmann theory. J. Phys. Chem. B 2007, 111, 5264−5276. (44) Abrashkin, A.; Andelman, D.; Orland, H. Dipolar PoissonBoltzmann equation: Ions and dipoles close to charge interfaces. Phys. Rev. Lett. 2007, 99, 077801. (45) Koehl, P.; Orland, H.; Delarue, M. Beyond the PoissonBoltzmann model: Modeling biomolecule-water and water-water interactions. Phys. Rev. Lett. 2009, 102, 087801. (46) Cheng, L.-T.; Wang, Z.; Setny, P.; Dzubiella, J.; Li, B.; McCammon, J. A. Interfaces and hydrophobic interactions in receptorligand systems: A level-set variational implicit solvent approach. J. Chem. Phys. 2009, 131, 144102. (47) Purisima, E. O.; Sulea, T. Restoring charge asymmetry in continuum electrostatics calculations of hydration free energies. J. Phys. Chem. B 2009, 113, 8206−8209. (48) Li, L.; Li, C.; Zhang, Z.; Alexov, E. On the dielectric ”constant” of proteins: Smooth dielectric function for macromolecular modeling and its implementation in DelPhi. J. Chem. Theory Comput. 2013, 9, 2126−2136. (49) Mukhopadhyay, A.; Aguilar, B. H.; Tolokh, I. S.; Onufriev, A. V. Introducing Charge Hydration Asymmetry into the Generalized Born Model. J. Chem. Theory Comput. 2014, 10, 1788−1794. (50) Zhou, S.; Cheng, L. T.; Dzubiella, J.; Li, B.; McCammon, J. A. Variational implicit solvation with Poisson-Boltzmann theory. J. Chem. Theory Comput. 2014, 10, 1454−1467. (51) Zhou, S.; Sun, H.; Cheng, L. T.; Dzubiella, J.; Li, B.; McCammon, J. A. Stochastic level-set variational implicit-solvent approach to solute-solvent interfacial fluctuations. J. Chem. Phys. 2016, 145, 054114. (52) Mobley, D. L.; Dill, K. A.; Chodera, J. D. Treating entropy and conformational changes in implicit solvent simulations of small molecules. J. Phys. Chem. B 2008, 112, 938−946. (53) Guthrie, J. P. A blind challenge for computational solvation free energies: introduction and overview. J. Phys. Chem. B 2009, 113, 4501−4507. (54) Geballe, M. T.; Skillman, a. G.; Nicholls, A.; Guthrie, J. P.; Taylor, P. J. The SAMPL2 blind prediction challenge: introduction and overview. J. Comput.-Aided Mol. Des. 2010, 24, 259−79. (55) Geballe, M. T.; Guthrie, J. P. The SAMPL3 blind prediction challenge: transfer energy overview. J. Comput.-Aided Mol. Des. 2012, 26, 489−496. (56) Knight, J. L.; Brooks, C. L. Surveying implicit solvent models for estimating small molecule absolute hydration free energies. J. Comput. Chem. 2011, 32, 2909−2923. N


Article

Journal of Chemical Theory and Computation (78) Michel, J.; Tirado-Rives, J.; Jorgensen, W. L. Prediction of the water content in protein binding sites. J. Phys. Chem. B 2009, 113, 13337−13346. (79) SZMAP, Version 1.1.1; OpenEye Scientific Software, Inc.: 2013. (80) Setny, P. Prediction of water binding to protein hydration sites with discrete, semi-explicit solvent model. J. Chem. Theory Comput. 2015, 11, 5961−5972. (81) Hummer, G.; Pratt, L. R.; García, A. E. Ion Sizes and Finite-Size Corrections for Ionic-Solvation Free Energies. J. Chem. Phys. 1997, 107, 9275−9277. (82) Hünenberger, P. H.; McCammon, J. A. Effect of artificial periodicity in simulations of biomolecules under Ewald boundary conditions: A continuum electrostatics study. Biophys. Chem. 1999, 78, 69−88. (83) Hünenberger, P. H.; McCammon, J. A. Ewald artifacts in computer simulations of ionic solvation and ion-ion interaction: A continuum electrostatics study. J. Chem. Phys. 1999, 110, 1856−1872. (84) Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects. J. Chem. Phys. 2013, 139, 184103. (85) Lin, Y. L.; Aleksandrov, A.; Simonson, T.; Roux, B. An overview of electrostatic free energy computations for solutions and proteins. J. Chem. Theory Comput. 2014, 10, 2690−2709. (86) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (87) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera − A Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−1612. (88) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; 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.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: The biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545−1614. (89) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (90) Bondi, A. van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441−451. (91) Sanner, M. F.; Olson, A. J.; Spehner, J.-C.; Sanner, M. D. Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 1996, 38, 305−320. (92) Connolly, M. L. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983, 221, 709−713. (93) Lee, B.; Richards, F. M. The interpretation of protein structures: estimation of static accessibility. J. Mol. Biol. 1971, 55, 379−400. (94) Lindorff-Larsen, K.; Maragakis, P.; Piana, S.; Eastwood, M. P.; Dror, R. O.; Shaw, D. E. Systematic validation of protein force fields against experimental data. PLoS One 2012, 7, e32131.

O


Explicit Solvent Hydration Benchmark for Proteins with Application to

Recommend Documents