Efficient Method To Characterize the Context-Dependent

Article pubs.acs.org/JPCB

Efficient Method To Characterize the Context-Dependent Hydrophobicity of Proteins Amish J. Patel*,† and Shekhar Garde*,‡ †

Department of Chemical and Biomolecular Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104, United States Howard P. Isermann Department of Chemical & Biological Engineering, and Center for Biotechnology & Interdisciplinary Studies, Rensselaer Polytechnic Institute, Troy, New York 12180, United States

‡

ABSTRACT: Characterizing the hydrophobicity of a protein surface is relevant to understanding and quantifying its interactions with ligands, other proteins, and extended interfaces. However, the hydrophobicity of a complex, heterogeneous protein surface depends not only on the chemistry of the underlying amino acids but also on the precise chemical pattern and topographical context presented by the surface. Characterization of such context-dependent hydrophobicity at nanoscale resolution is a nontrivial task. The free energy, μex v , of forming a cavity near a surface has been shown to be a robust measure of context-dependent hydrophobicity, with more favorable μex v values suggesting hydrophobic regions. However, estimating μex v for cavities significantly larger than the size of a methane molecule in a spatially resolved manner near proteins is a computationally daunting task. Here, we present a new method for estimating μex v that is 2 orders of magnitude more efficient than conventional techniques. Our method envisions cavity creation as the emptying of a volume of interest, v, by applying an external potential that is proportional to the number of water molecules, Nv, in v. We show that the “force” to be integrated to obtain μex v is simply the average of Nv in the presence of the potential, and can be sampled accurately using short simulations (50−100 ps), making our method very efficient. By leveraging the efficiency of the method to calculate μex v at various locations in the hydration shell of the protein, hydrophobin II, we are able to construct a hydrophobicity map of the protein that accounts for topographical and chemical context. Interestingly, we find that the map is also dependent on the shape and size of v, suggesting an “observer context” in mapping the hydrophobicity of protein surfaces.

■

INTRODUCTION Water plays a crucial role in mediating biomolecular interactions, in particular through nonspecific hydrophobic effects. However, characterizing protein hydrophobicity (and consequently, interactions) is a challenging task, with numerous implications, including on the search for therapeutic remedies1−3 and the development of cellular protein interaction networks.4 While the hydrophobicity of a macroscopic surface is typically characterized by the water droplet contact angle,5 such characterization is not feasible for the rugged, nanoscopic surfaces of proteins and poses a significant challenge.1 Further, hydrophobicity depends not only on the chemistry of the underlying surface but also on the particular context presented by the protein: hydrophobicity is influenced by surface topography6−8 and chemical patterning9−13 and can depend nontrivially on the specific combination of the two.14−18 As a result, approaches that do not incorporate information about how water at an interface responds to the local context, for example, some of the popular hydropathy scales,19,20 are not expected to yield accurate information about specific protein− ligand interactions and may have limited utility for such purposes.12,21,22 © 2014 American Chemical Society

Recent theoretical and simulation studies have shown that water density fluctuations, and fluctuations that result in the formation of a cavity in particular, are enhanced near hydrophobic surfaces, making them a suitable metric for quantifying the context-dependent hydrophobicity of complex surfaces.23−27 In fact, for flat self-assembled monolayer surfaces with a range of headgroup chemistries, Patel et al.28 demonstrated a quantitative connection between the cavity formation free energy near a particular surface and the water droplet contact angle on that surface.29 Such a connection suggests that cavity hydration free energies not only agree with our macroscopic notions of hydrophobicity but also have the added advantage that, unlike contact angles, they can be measured in the vicinity of rugged surfaces and, thereby, be used to characterize protein hydrophobicity. Central to such context-dependent characterization is the ability to compute the hydration free energies, μex v , of cavities of all shapes and sizes, v. For a small (e.g., methane-sized) v, Received: August 15, 2013 Revised: December 8, 2013 Published: January 3, 2014 1564

dx.doi.org/10.1021/jp4081977 | J. Phys. Chem. B 2014, 118, 1564−1573

The Journal of Physical Chemistry B

Article

thermal fluctuations result in the cavity forming spontaneously, and μvex can be estimated easily from typical molecular simulations by performing test particle insertions.30−32 Indeed, Acharya et al. have used the hydration free energies of small methane-sized cavities near a protein to map its hydrophobicity.12 However, because the likelihood of spontaneous cavity formation decreases roughly exponentially with its size, estimating μvex for larger cavities requires computationally intensive biased sampling techniques, such as umbrella sampling.33 Further, the calculation of μex v at a large number of locations on the protein surface, and possibly for several sizes/shapes of v, is required to map its hydrophobicity, making the overall problem computationally intractable. Here we present a method for estimating μex v for cavities of arbitrary shapes and sizes that is about 2 orders of magnitude more efficient than conventional techniques (e.g., umbrella sampling or free energy perturbation). We exploit the fact that creating a cavity of a given shape and size is equivalent to removing all the waters from the corresponding probe volume, v. By applying an external biasing potential, ϕNv, that couples linearly to the number of water molecules, Nv, in v, we are able to systematically empty v and create a cavity (by increasing ϕ). By choosing the unfavorable potential to be linear in Nv, we are able to derive a simple, yet powerful expression that connects μex v to the dependence of the average number of waters, ⟨Nv⟩ϕ, on ϕ. We show that ⟨Nv⟩ϕ decreases monotonically with ϕ and can even be linear in ϕ. As a result, excellent estimates of μex v can be obtained from the knowledge of ⟨Nv⟩ϕ at only a few values of ϕ. In addition, since typical correlation times in liquid water are small (∼1 ps), Nv converges to its equilibrium value quickly (∼50 ps), allowing estimation of μex v using only a few, short simulations. By leveraging the efficiency of the method, we are able to characterize the hydrophobicity of a protein−water interface and find that large parts of the protein surface are superhydrophilic (more hydrophilic than a completely wetting surface). Further, by estimating μex v for different probe volumes, we find that the protein hydrophobicity map can depend on the size and the shape of the cavity. Thus, the cavity hydration free energy at a given location not only captures the local context but also highlights the presence of an “observer” context in the characterization of protein hydrophobicity. In the following section, we first derive the central equation underlying our method. We then present details pertaining to the systems studied and the simulations employed. In the subsequent section, we first use a toy-model system to demonstrate how the method works and to illustrate the reasons underlying its efficiency. We then apply the method to characterize the hydrophobicity of self-assembled monolayer surfaces and to create hydrophobicity maps for the protein hydrophobin II. We conclude with a discussion of opportunities where characterizing the context-dependent hydrophobicity of protein surfaces may provide valuable insights. In the Appendix, we discuss practical issues related to the implementation of the method and ways of minimizing errors.

is the generalized Hamiltonian in the absence of ϕ. Let Qϕ be the partition function associated with Hϕ. Then, the nth derivative of ln Qϕ with respect to −βϕ is the nth cumulant of the water number distribution. Thus, ∂ ln Q ϕ

= ⟨Nv⟩ϕ

∂( −βϕ) ∂ 2 ln Q ϕ 2

∂( −βϕ)

=

(1)

∂⟨Nv⟩ϕ ∂( −βϕ)

= ⟨δNv2⟩ϕ

(2)

and so on. Integrating eq 1, we get Qϕ Q0

⎡ = exp⎢ −β ⎣

∫0

ϕ

⎤ ⟨Nv⟩ϕ′ dϕ′⎥ ⎦

(3)

It can easily be shown that Qϕ/Q0 is also related to the probability, Pv(N), of observing N water molecules in v in the absence of ϕ (i.e., Pv(N) ≡ ⟨δN,Nv⟩0) as Qϕ Q0

=

∑ Pv(N ) e−βϕN

(4)

N

which has the form of a generating function. By defining z ≡ e−βϕ, we get Pv(0) =

Pv(N ) =

Qϕ Q0

and (5)

z → 0, ϕ →∞

1 ∂N ⎛ Q ϕ ⎞ ⎜⎜ ⎟⎟ N! ∂z N ⎝ Q 0 ⎠

(6)

z→0

Combining eqs 3 and 5 with central result:

βμex v

= −ln Pv(0)

30

yields our

Equation 7 is in the spirit of thermodynamic integration, where ⟨Nv⟩ϕ is the force (see eq 1), and the integration is performed over a range of ϕ values from zero (i.e., the unbiased system) to infinity. Increasing ϕ penalizes the presence of water molecules in v, leading to the creation of a cavity in the limit of ϕ → ∞. We note that although we have chosen to focus on evaluating the free energy, μvex, required to empty an observation volume by displacing all water molecules in it, the method presented here is fairly general and can easily be adapted to compute the free energy of driving any order parameter to zero. Simulation Details. We illustrate the use of eq 7 to estimate μex v in interfacial environments using three classes of systems, which in order of increasing complexity are (1) toymodel surfaces that are characterized by the statistics of water number distributions, Pv(N), in their vicinity; (2) selfassembled monolayers (SAMs) of alkyl chains with headgroup chemistries that govern the hydrophobicity of the SAM-water interface;25,34 and (3) the hydration shell of the protein, hydrophobin II.35 To study classes 2 and 3 above, we use molecular dynamics (MD) simulations using a suitably modified version of the GROMACS package,36 which allows us to apply an external potential (see the discussion in the Appendix on “Reweighting” for the exact form of the applied potential).26,37 In all

■

THEORY AND METHODS Theoretical Derivation of the Principal Result. Consider an observation volume, v, described by its shape, size, and location; let Nv be the number of water molecules in v. To systematically empty v and create a cavity, we apply an unfavorable potential, ϕ, that couples linearly to Nv so that the Hamiltonian of the system becomes Hϕ = H0 + ϕNv, where H0 1565



Article

force field Lennard-Jones parameters.46 Locations of the benzene-shaped observation volumes in the protein hydration shell were arrived upon as follows: We first evaluate the instantaneous interface corresponding to the protein heavy (non-hydrogen) atoms following ref 51. We then determine the surface normal at every point on the instantaneous interface. Finally, corresponding to a particular point on the instantaneous interface, the benzene-shaped v is placed with its center on the surface normal, in a plane perpendicular to the normal, and as close to the interface as possible without any of the protein heavy atom centers residing in v. Before performing any simulations with a biasing potential, we first equilibrated an unbiased protein water system for 5 ns. Biased simulations with the application of a linear potential were then run for 1.2 ns, and data were collected at an interval of every 0.1 ps. In most cases, the number of water molecules in the observation volumes reached a steady value within the first 50 ps. We discarded the data from the first 200 ps and used those from the subsequent 1 ns to calculate μex v . The free energy of hydrating the benzene-shaped cavity in bulk water was also computed in a similar manner. To learn whether the hydrophobicity of the underlying protein surface depends on the size and shape of the probe volume used to interrogate it, we also independently computed μex v for small spherical volumes (radius = 2.5 Å) placed at the centers of the benzene-shaped volumes. For these small volumes, μex v could be evaluated simply by running a single 5 ns, unbiased protein water simulation (data collected every 1 ps) and using the test particle insertion method.30

simulations, periodic boundary conditions were employed, the leapfrog algorithm38 was used to integrate the equations of motion with a 2 fs time-step, and the particle mesh Ewald algorithm39 was used to treat long-range electrostatic interactions. Water was represented using the SPC/E model,40 and bonds involving hydrogen atoms were constrained using the SHAKE41 algorithm. In each system, we select an observation volume of a given shape and size near the surface of interest, perform biased simulations at different ϕ values, and estimate μex v using the method derived in the previous section. The cavity hydration free energy, thus obtained, yields insightful information about the hydrophobicity of the underlying surface. System specific details, including those pertaining to the placement of observation volumes at interfaces, are given below. SAM Surfaces. This system presents a simple but realistic class of surfaces that are uniform and flat and for whom contact angles can be readily measured. Our simulation setup is the same as that used in ref 28 and has been described in detail previously.25,34 All calculations were performed in the NVT ensemble with a buffering vapor−liquid interface26,37,42 and at a temperature of 300 K, maintained using the Berendsen thermostat.43 μex v values have been measured previously for cuboid-shaped volumes near SAM surfaces28 using the INDUS method.37 Here, we quantify μex v near these surfaces using the new method to (i) validate the method by resolving differences in the hydrophobicities of these realistic surfaces, and (ii) demonstrate the computational efficiency of our approach. A 2 × 2 × 0.3 nm3 cuboid observation volume is placed either in bulk water or immediately adjacent to the SAM surface such that the smallest dimension (0.3 nm) is perpendicular to the SAM surface (see Figure 3a), and the center of the volume in this dimension is in the plane corresponding to the maximal water density. Hydrophobin II. Hydrophobin II (crystal structure, PDB ID: 2B97),35 is a fungal protein that accumulates as films at interfaces (e.g., the air−water interface) and is known to reverse the wetting properties of the interface. It does so by having two faces, one hydrophobic and one hydrophilic, akin to a Janus particle.44,45 Whereas most folded, globular proteins are expected to be predominantly hydrophilic, hydrophobin II serves as an ideal model protein for our studies because its hydration shell is expected to display a wide range of hydrophobicities. We simulate a single chain of the hydrophobin II dimer, represented using the AMBER 94 force field46 and solvated by roughly 11 000 SPC/E water molecules, in a 7 × 7 × 7 nm3 periodic box. All calculations were performed in the NPT ensemble (T = 300 K, P = 1 bar) using the canonical velocity rescaling thermostat47 and the Parrinello−Rahman barostat.48 To prevent diffusion of the protein while still allowing some thermal motion of the side chains, all the Cα atoms in the protein were position-restrained harmonically with a relatively soft spring constant of 1000 kJ/nm2 in each dimension. To characterize the hydrophobicity of the protein−water interface, we used benzene-shaped volumes, v, in the hydration shell of hydrophobin II that complement the shape of the protein surface (see Figure 4). The C−C and C−H cavity distances were taken to be 1.4 and 1.08 Å, respectively, corresponding to the respective benzene bond lengths in the AMBER force field.46 The effective cavity radii for the C and H atoms were taken to be 3.11 and 2.6 Å, respectively, obtained by following the WCA prescription49,50 and using the AMBER

■

RESULTS AND DISCUSSION Demonstrating the Method Using Toy Model Surfaces. To illustrate how ⟨Nv⟩ϕ responds to ϕ in various interfacial environments, we first consider a family of toy-model surfaces with increasing hydrophobicity. These surfaces are defined by Pv(N) distributions in their vicinity and are parametrized by λ, Pv(N ; λ) = k e−0.5n

2

− λn3[1 − /(n)]

n ≡ (N − n0)/σ

(8)

where / (n) is the Heaviside step function and k is the normalization constant. λ = 0 corresponds to a Gaussian distribution, with mean, n0, and standard deviation, σ. The lifting of the low-N tail with increasing λ for the family of Pv(N; λ) curves shown in Figure 1a mimics the behavior of Pv(N) measured in cuboidal volumes near increasingly hydrophobic surfaces.26,52 In the presence of a linear biasing potential, ϕNv, the statistics of Nv change from Pv(N; λ) to cPv(N; λ) exp (−βϕN), where c is a normalization constant.52,53 ⟨Nv⟩ϕ is the first moment of the biased distribution function and is shown in Figure 1b as a function of ϕ. Because ϕ penalizes the presence of water in v, ⟨Nv⟩ϕ decreases as ϕ increases. In fact, because ⟨δNv2⟩ϕ is always positive, eq 2 confirms that the slope, ∂⟨Nv⟩ϕ/∂ϕ, is always negative, resulting in a monotonic decrease of ⟨Nv⟩ϕ with ϕ. For λ = 0, Pv(N; λ) is Gaussian, and the corresponding ⟨Nv⟩ϕ decreases linearly, as expected from linear response theory. For the more hydrophobic surfaces (higher λ), ⟨Nv⟩ϕ decreases more rapidly with increasing ϕ, resulting in smaller areas under the curve. Thus, the more hydrophobic the surface, the smaller the value of μex v , which is consistent with the ease of cavity formation near 1566



Article

Figure 2. (a) Time dependence of the number of water molecules, Nv, in a cuboid observation volume, v = 2 × 2 × 0.3 nm3, in bulk water obtained using MD simulations of SPC/E water40 in the presence of biasing potentials with ϕ = 0, 2, and 4 kBT (gray lines). The running averages (colored lines) converge to their equilibrium values, ⟨Nv⟩ϕ (dashed lines), in ≈50 ps. (b) Estimates of μex v for the toy-model surfaces (colors are the same as Figure 1) as a function of the number of points, S, used to perform the integration in eq 7. For each S, the set {ϕi} (i = 1, 2, ..., S) was chosen according to eq 9.

Figure 1. (a) A family of curves that mimic water number distributions near increasingly hydrophobic surfaces, and (b) the corresponding response of ⟨Nv⟩ϕ to ϕ. The curves correspond to eq 8 with n0 = 100 and σ = 4, chosen to be similar to those for a cuboid v = 3 × 3 × 0.33 nm3, and 103λ = 0 (black), 2 (blue), 4 (cyan), 6 (magenta), and 8 (red). The arrows indicate the direction of increasing λ (or hydrophobicity). (c) Comparison of μex v calculated using eq 7 with its expected value, −kBT ln Pv(0). The axes range from 175 to 325 kBT.

hydrophobic surfaces.25,26 Figure 1c confirms that for each λ surface, μex v calculated using eq 7 exactly matches its expected value, −kBT ln Pv(0). Why The Method is Efficient: Inherently Small Errors. In the above toy-model example, the entire Pv(N) distributions were known a priori, which permitted ⟨Nv⟩ϕ to be known accurately at all values of ϕ. This results in excellent estimates of the area under the ⟨Nv⟩ϕ vs ϕ curve. In practice, the underlying Pv(N) distributions are not known, and ⟨Nv⟩ϕ is measured by performing biased simulations of finite length at a finite number of ϕ values. As a result, there are two sources of error in our estimate of μex v obtained using eq 7: (i) the error in our estimate of ⟨Nv⟩ϕ due to a finite simulation time and (ii) the numerical integration error resulting from the fact that estimates of ⟨Nv⟩ϕ are available only at a few values of ϕ. Below, we discuss why both these errors are expected to be small in typical situations. Finite Simulation Errors. Errors arising from the uncertainty in ⟨Nv⟩ϕ are generally small because typical water number distributions are unimodal and peaked around their mean value26,31,52 (see caveat below). The observable of interest, ⟨Nv⟩ϕ, is then determined by the most probable region of the water number distribution, resulting in accurate estimates of ⟨Nv⟩ϕ from short simulations (≈50−100 relaxation times). For example, the time dependences of Nv in a cuboidal v in bulk water, shown in Figure 2a, confirm that the running averages indeed converge rapidly (≈50 ps), being governed by the fast (≈1 ps) diffusion of water molecules across the boundaries of v. In contrast, errors in estimating μex v using umbrella sampling originate from insufficient sampling of the tails of biased order parameter distributions in the overlap region between two neighboring windows. The need for sampling such tails makes the umbrella sampling method expensive. The quantity of interest in our method is ⟨Nv⟩ϕ, which is determined by the most probable region of the biased water number distributions, not by the improbable tails, making our method efficient. Caveat An important exception arises when the water number distribution in the presence of the biasing potential becomes bimodal with a large (compared to kBT) barrier separating the low and high-density basins. Under these circumstances, short simulations will not be able to sample phase space ergodically, resulting in systematic errors in the estimates of ⟨Nv⟩ϕ and, consequently, in the estimate of μex v . Fortunately, the conditions under which a biased water number distribution is bimodal are

expected to arise infrequently in characterizing protein hydrophobicity, and even then, the distribution will be bimodal only over a small range of ϕ values. The Appendix contains a detailed discussion on when the biased water number distribution is likely to be bimodal; what its indicators are; how to estimate the resultant error in μex v ; and how to estimate μex v more accurately, if needed. Numerical Integration Errors. The numerical integration errors in estimating μex v using eq 7 depend on the number of biased simulations that are performed, that is, the number of ϕ values at which ⟨Nv⟩ϕ is known, with the error decreasing as more simulations are performed. When the unbiased water number distribution, Pv(N), is Gaussian (e.g., near a hydrophilic surface), the response, ⟨Nv⟩ϕ, is linear in ϕ (see Figure 1, λ = 0 or Figure 3, −OH-terminated SAM). Thus, estimates of ⟨Nv⟩ϕ at only two ϕ values, or estimates of ⟨Nv⟩ϕ, and its derivative, −β⟨δNv2⟩ϕ, at one ϕ value (as was done in ref 31 at ϕ = 0 using an information-theory-based approach), would be sufficient to estimate μex v . However, for hydrophobic surfaces, the underlying

Figure 3. (a) A snapshot of part of the MD system containing a selfassembled monolayer in water. Alkane tail (cyan); −CH3 head groups (blue and white); water (red and white); and a cuboidal volume, v, of interest (white) at the interface are shown. (b) The response of ⟨Nv⟩ϕ (averaged over 250 ps) to ϕ for v near the −CH3, −OCH3, and −OH ex SAM surfaces. The shaded area equals μex v for the −CH3 SAM. (c) μv obtained by using an increasing number of total simulations, S. For each S, the set {ϕi} (i = 1, 2, ..., S), was chosen according to eq 9, and 50 ps of data (after equilibration) from each simulation was used to estimate ⟨Nv⟩ϕ. 1567



Article

statistics, as quantified by Pv(N), are Gaussian only near the mean and have fat, low-N tails.23,26,27,52 As a result, ⟨Nv⟩ϕ varies sensitively with ϕ, and to obtain the best estimate of μex v , it is important to sample ϕ values in the region over which ⟨Nv⟩ϕ varies rapidly. Indeed, there is freedom to choose not only the number of biased simulations to be performed but also the ϕ values at which they are performed. To this end, an adaptive strategy, wherein information from existing biased simulations is used to select the ϕ value for the next simulation, might be useful. For example, if S biased simulations are to be performed to estimate μex v , we propose the following adaptive scheme to select a set of ϕ values ({ϕi}, i = 1, 2, ..., S). ϕ1 = 0, and the estimate of the slope, ∂⟨Nv⟩ϕ/∂ϕ = −β⟨δNv2⟩ϕ at ϕi informs the value of ϕi+1 as follows: βϕi + 1 = βϕi + [⟨Nv⟩0 /S]/⟨δNv 2⟩ϕi

speed-up is expected to be greater near a hydrophilic surface than near a hydrophobic surface for two reasons: (i) The dependence of ⟨Nv⟩ϕ on ϕ is simpler (approximately linear) near a hydrophilic surface compared with that near a hydrophobic surface (sigmoidal), so fewer simulations are ex needed to estimate μex v accurately. (ii) μv is larger near a hydrophilic surface, and Pv(0) = exp(−βμex v ) is therefore smaller, which makes the umbrella sampling method less efficient. In particular, harmonic potentials centered at smaller Nv values are needed to sample the region near N = 0 accurately using umbrella sampling. Thus, more umbrella simulations (i.e., windows) will be required to obtain μex v near a hydrophilic surface. Further, as the size of v increases, so does the average number of waters in v, and a larger number of umbrella simulations are again required. Thus, our method is expected to be increasingly more efficient near more hydrophilic surfaces, as well as for larger v’s near all surfaces. Application: Painting the Protein Hydrophobin II. We use the cavity formation free-energy-based measure to paint hydrophobicity maps of a protein hydrophobin II (PDB ID: 2B97),35 for benzene-shaped and spherical probes (see Figure 4). As described in the Methods section, we first identify an exact location of the observation volume, v, in the protein hydration shell. To calculate μex v at that location, we run multiple simulations with biasing potentials that are linear in Nv and integrate the response, ⟨Nv⟩ϕ, with respect to ϕ using eq 7. We similarly calculate μex v for different locations of v that span the protein hydration shell. The number of probe locations and their concentrations in specific protein hydration shell regions may be selected depending on the problem of interest. In our calculations, we used 154 evenly spaced locations of v. For the selected locations, we use the ratio of the hydration shell, μex v , normalized by its value in bulk water, μex , as our measure of bulk hydrophobicity. For extended flat surfaces, this measure can be directly related to the water droplet contact angle on the surface.28 For example, for a thin cuboid volume with a large ex cross-sectional area near a flat surface, μex v /μbulk can be shown to reduce to 0.5(1 + cos θ) and ranges from 0 (hydrophobic, θ = π) to 1(hydrophilic, θ = 0).29 Such an analogy for protein surfaces is not necessary, but provides a way to think about spatially resolved wettability from a macroscopic frame of reference. ex Figure 4a shows the distribution of μex v /μbulk obtained from data at 154 different locations. The distribution is peaked roughly near 1, suggesting that the protein surface is largely hydrophilic, and the tails of the distribution correspond to hydrophobic (ratio < 1) and increasingly hydrophilic (ratio > 1) regions. Because of an abundance of hydrophilic groups, typical protein surfaces are expected to display a relatively narrow range of hydrophilic/phobic behavior. Thus, we do not paint the protein surface directly using the μex v data normalized ex by μbulk . Instead, to tease out useful information using hydrophobicity maps, we use coloring schemes that assign a number between 0 (red, hydrophobic) and 1 (blue, hydroex philic) corresponding to the local value of μex v /μbulk. The functional forms of the two coloring schemes that we employ, 1-erf and 2-erf, are also shown in Figure 4a. The 1-erf coloring scheme highlights differences across the most probable region ex of the P(μex v /μbulk) distribution, capturing broader changes in hydrophobicity, whereas the 2-erf scheme focuses on the tails of the distribution, highlighting the most hydrophobic and the most hydrophilic regions on the protein surface.

(9)

For linear response, ⟨δNv2⟩ϕi will be constant, and {ϕi} selected using eq 9 will be equidistant, with ⟨Nv⟩ϕS ≈ 0. For the sigmoidal response characteristic of hydrophobic surfaces, the chosen ϕ values will be spaced close together (i.e., small ϕi+1 − ϕi) where the slope is steep (−∂⟨Nv⟩ϕ/∂ϕ = β⟨δNv2⟩ϕ is large) and farther apart where it is shallow. Although better adaptive schemes can be developed, Figure 2b shows that as few as three ϕ values picked according to eq 9 are sufficient to obtain excellent estimates of μex v and differentiate surfaces of different hydrophobicities (i.e., λ values). Thus, the relatively simple, monotonic dependence of ⟨Nv⟩ϕ on ϕ (Figure 1b) means that the integral in eq 7 can be evaluated numerically using estimates of ⟨Nv⟩ϕ from just a few simulations, making our method particularly efficient. Application: Self-Assembled Monolayer Surfaces. Figure 3 illustrates the use of eq 7 to estimate cavity hydration free energies near self-assembled monolayers (SAMs) with three different headgroup chemistries: −CH3, −OCH3, and −OH. The response, ⟨Nv⟩ϕ, to the potential, ϕNv, applied in a 2 × 2 × 0.3 nm3 cuboidal volume located at these interfaces is shown in Figure 3b. Near the hydrophilic −OH surface, the response is approximately linear, suggesting that the statistics of the unbiased water number distribution must be close to Gaussian, consistent with previous studies.26,52 Near the −OCH3 SAM, the response is linear for small ϕ but dips at ϕ ≈ kBT, suggesting that the underlying statistics are Gaussian near the mean but display a fat tail at lower N values.52 Near the −CH3 surface, the response is qualitatively similar to that near the −OCH3 SAM but becomes sigmoidal at a smaller ϕ value. Even though the average Nv values at ϕ = 0, ⟨Nv⟩0, are similar near each surface, the response to ϕ is markedly different near each surface at larger ϕ and depends on surface hydrophobicity: the more hydrophobic the interface, the smaller the value of ⟨Nv⟩ϕ. That is, the unfavorable biasing potential, ϕNv, is more effective in displacing waters from the vicinity of surfaces that are more hydrophobic. As a result, the area under the ⟨Nv⟩ϕ-vsϕ curve is smallest near the hydrophobic −CH3 SAM, and it is easiest to solvate a cavity in its vicinity. Importantly, as shown in Figure 3c, as few as 2−3 short (50 ps) simulations at values of ϕ prescribed by eq 9 are sufficient to obtain excellent estimates of μex v . In contrast, typical umbrella sampling calculations for this v require 15−20 windows with roughly 500 ps simulations for each window to ensure sufficient overlap, implying that our method is about 2 orders of magnitude more efficient. The exact speed-up in calculating μex v will depend on the particular size, shape, and location of v. The 1568



Article

contact angle of zero.29 Interestingly, for large parts of the protein (blue patches in Figures 4b and c), μex v values are larger than μex bulk, implying that these regions of the protein surface can be thought of as being “superhydrophilic”, that is, they are even more hydrophilic than typical fully wetting surfaces. This highlights another advantage of the μex v -based mapping: while the contact-angle-based measure of hydrophilicity is bounded by a contact angle of zero, μex v can differentiate between ex ex ex hydrophilic (μex v ≈ μbulk) and superhydrophilic (μv > μbulk) surfaces. Using the same data, the 2-erf color scheme highlights patches with extreme values of μvex, thereby identifying hydrophobic (and hydrophilic) hot-spots (Figure 4d). Although the 1-erf scheme identified a large hydrophobic face of the protein (Figure 4b), it is only weakly hydrophobic, and the 2-erf color scheme reveals that specific locations on that face are significantly more hydrophobic. Those locations are likely to have a higher propensity to bind hydrophobic ligands and may serve as nucleation sites in protein aggregation. An experimental test of our calculations would involve mutating the hydrophobic hot-spots to hydrophilic or ionic residues, which would significantly hamper the ability of hydrophobin II to bind to hydrophobic interfaces. Even a single mutation in this region may have a dramatic influence on hydrophobin II function. Observer Context: Hydrophobicity Map Depends on the Probe. We also investigate how the size and shape of the cavity probes affect the resulting hydrophobicity maps by constructing maps for benzene-shaped (nonspherical) as well as small spherical (radius = 2.5 Å) probes (see Methods). Figure 4e shows several instances in which the two hydrophobicity color ex maps disagree significantly, indicating that μex v /μbulk (and, therefore hydrophobicity at a particular location on the protein surface) depends on the size and shape of the probe. Although more extensive calculations would have been needed to explore this aspect completely, the present calculations highlight an important aspect of hydrophobicity maps: namely, the “observer context”, that is, the hydrophobicity of a protein surface depends on the size/shape of the cavity used to probe ex it. Although the general trend for μex v /μbulk is similar in both cases, the scatter in Figure 4f means that certain locations on the protein surface that strongly bind benzene will not bind small spherical solutes and vice versa. We note that our cavity probes are hard objects, without dispersive interactions. Realistic solutes, such as benzene or methane, also interact with the protein and water via van der Waals interactions. For specific solutes of interest, the effect of dispersive interactions can be incorporated efficiently using ideas from perturbation theory.33,54 Nevertheless, it is wellappreciated that the cavity formation free energy constitutes the dominant contribution to the hydrophobic solvation free energy. Thus, the mapping of hydrophobin II as well as the presence of observer context is not expected to be significantly affected by the addition of van der Waals interactions.

Figure 4. Mapping the hydrophobicity of the protein hydrophobin II (shown in black) using spatially resolved cavity hydration free energies. ex (a) The distribution of μex v /μbulk values for the benzene-shaped probe volume (gray bars) and a Gaussian fit to the distribution. Basis of the color schemes 1-erf and 2-erf are also shown. On the basis of its μex v / μex bulk value, the coloring scheme assigns to every v on the protein surface a number between 0 and 1 that determines its color. Small ex values of μex v /μbulk (hydrophobic) are assigned numbers closer to 0 and ex are painted red, whereas large values of μex v /μbulk (hydrophilic) are assigned numbers closer to 1 and are painted blue. Intermediate numbers (near 0.5) are painted white (see color bar between panels b and c). (b, c) Hydrophobicity maps of two sides of the protein using a benzene-shaped v (shown by sticks) using the 1-erf scheme. Panel d shows a map using the 2-erf color scheme to highlight locations corresponding to the extreme values of μex v . (e) Comparison of the maps obtained using the benzene-shaped v with that using a smaller spherical v. Both are colored using the 1-erf scheme. (f) Scatter in the ex μex v /μbulk values for the benzene-shaped and the spherical v.

Hydrophobicity maps that emphasize different aspects of protein hydrophobicity are shown in Figure 4(b−e). The hydrophobicity map shown in Figure 4b using the 1-erf scheme shows a relatively large hydrophobic region identified by the large contiguous red patch. The presence of this large patch is consistent with the function of hydrophobin II: it is thought to help the protein bind to vapor−liquid and other hydrophobic surfaces during spore formation.35 Figure 4c shows the 1-erf mapping of this protein in another orientation, highlighting that much of the rest of the protein surface is covered with a ex waterlike hydrophilic surface (white), with μex v ≈ μbulk. Previous 25,28 ex ex studies of μv near flat surfaces have shown that μex v ≈ μbulk corresponds to surfaces that are fully wetting, that is, have a

■

OUTLOOK We have presented an efficient computational method to develop hydrophobicity maps of proteins and other complex heterogeneous nanoscopic surfaces. The method uses contextdependent cavity formation free energies in the vicinity of proteins as a measure of hydrophobicity. It includes the “observer context” because the free energies are dependent on the size, shape, and orientation of the probe. Our method is 1569



Article

Figure 5. Demonstrating the asymptotic behavior of ⟨Nv⟩ϕ as ϕ → ∞ for the toy-model surfaces. (a) For all surfaces, log⟨Nv⟩ϕ + βϕ plateaus at large ϕ, confirming that ⟨Nv⟩ϕ ∼ e−βϕ. (b) The plateau values are indeed given log[Pv(1)/Pv(0)], as shown in eq 11. (c) The negative of the slope of log⟨Nv⟩ϕ + βϕ, given by ⟨δNv2⟩ϕ/⟨Nv⟩ϕ − 1 (see eq 14), must approach 0 as ϕ → ∞. The color scheme is the same as that in Figure 1.

encapsulated in reverse micelles,62,63 the use of molecular probes combined with NMR64 or IR spectroscopy,65 and the multiple solvent crystal structures or solvent mapping method,66 which employs differential X-ray crystallography of neat protein crystals with those of nonpolar ligands bound to proteins.

potentially useful in obtaining fundamental insights and also has applied relevance. The efficiency of our method may allow systematic sampling of the space of parameterssolute curvature, chemical heterogeniety, topography, etc.to facilitate a fundamental understanding of the context dependence of hydrophobicity. Such an understanding may lead to the development of predictive correlations relating properties of a complex surface to its context-dependent hydrophobicity. An understanding of the context dependence of hydration and hydrophobicity is also important in numerous applications. A few mutations at specific locations in the complementarity determining regions of antibodies have been shown to change their propensity to aggregate.55 Locations for such mutations are typically arrived at after a long screening process through sequence space. An understanding of the context dependence of hydration may help focus the locations and residues for such specific mutations. The Protein Data Bank currently contains over 90 000 structures. Those structures, however, have limited information about the hydration of proteins (except for the presence of crystal waters). The efficiency of our method and its amenability to parallelization may also allow the development of a hydration database containing hydrophobicity maps of many proteins or even classes of proteins for a spectrum of probe solutes. Hydrophobicity maps of a class of proteins may reveal aspects that are important to the interactions and function of that class. Further, studying classes of proteins that perform certain desirable functions could reveal the underlying patterns responsible and inform the biomimetic design of functional surfaces, for example, nonfouling surfaces,56 or surfaces that are able to bind and detect ligands/toxins of interest. The most direct application of our method is to estimate binding preferences of a hydrophobic solute of a given shape and size to a protein of interest. Binding of solutes to proteins is important in drug design3 as well as in chromatographic separations through ligand design.57 Although computational methods can provide in-depth information on protein interactions in a highly efficient manner, complementary experimental studies of protein− ligand interactions will be needed to validate the findings. To this end, recent experimental studies that report precise measurements of hydration or binding free energies and entropies using techniques such as isothermal titration calorimetry,22,58 microfluidic binding assays,59 and single molecule AFM pulling60,61 are promising. Other innovative experimental strategies that aim to characterize these interactions indirectly include NMR measurements on proteins

■

APPENDIX

Reweighting

Nv changes discontinuously as molecules cross the boundaries of v. As a result, employing a Hamiltonian, which is a function of Nv (e.g., Hϕ = H0 + ϕNv) would lead to impulsive forces in MD simulations.26,37 To circumvent this challenge, we perform our simulations in an ensemble defined by H′ϕ = H0 + ϕÑ v, where Ñ v is a coarse-grained variable that is strongly correlated with Nv and varies continuously across v. For the definition of Ñ v that we employ here, see ref 37. Thermal averages in the ensemble defined by Hϕ can then be obtained by reweighting as ⟨X ⟩ϕ =

⟨X exp[−βϕ(Nv − Nṽ )]⟩′ϕ ⟨exp[−βϕ(Nv − Nṽ )]⟩′ϕ

(10)

where ⟨...⟩′ϕ represents an average in the ensemble defined by Hϕ′ . We note that the right hand side of eq 10 involves averages of exponentials. From a practical standpoint, it can require a long time for these averages to converge,67 especially if the exponents are large (≫1). This is typically not the case, because Nv is strongly correlated with Ñ v, and the difference between Ñ v and Nv that appears in the exponents is typically small. In addition, by reducing the coarse-graining length used to define Ñ v,37 this difference can be made even smaller. The exponents in eq 10 also depend on βϕ. Fortunately, for large values of ϕ, ⟨Nv⟩ϕ is already very small and does not contribute appreciably to μex v (see eq 7). Further, as shown below, an understanding of the asymptotic behavior of ⟨Nv⟩ϕ in the limit of large ϕ obviates the need for simulations at large ϕ values. Extrapolating to ϕ → ∞

In practice, when calculating μex v using eq 7, ⟨Nv⟩ϕ is computed at a few values of ϕ, and the integral is estimated using standard numerical integration procedures. However, the integral in eq 7 runs from 0 to ∞, necessitating an extrapolation of ⟨Nv⟩ϕ from the maximum value of ϕ at which a simulation was run, ϕ = ϕmax, to ϕ = ∞. Further, as discussed above, estimates of ⟨Nv⟩ϕ obtained by reweighting become less reliable for large values of ϕ. Here, we resolve both these issues by shedding light on the 1570



Article

Figure 6. (a) The susceptibility, ∂⟨Nv⟩ϕ/∂(−βϕ) = ⟨δNv2⟩ϕ, for the toy-model surfaces. For hydrophobic surfaces, there is a peak in susceptibility, which increases with the hydrophobicity of the surface. (b) The biased distribution function, cPv(N; λ)e−βϕ*N, where ϕ* corresponds to the peak in the susceptibility, ∂⟨Nv⟩ϕ/∂(−βϕ), shown in part a, for the two most hydrophobic surfaces; λ = 6 × 10−3 (magenta) and λ = 8 × 10−3 (red). Although the biased distribution function is unimodal for moderately hydrophobic surfaces, it can become bimodal for highly hydrophobic surfaces, with barriers separating the high- and low-density basins. (c) Even for the most hydrophobic surface (λ = 8 × 10−3), the biased distribution function is bimodal only for values of ϕ near ϕ* = 2.57kBT. For ϕ values that correspond to the susceptibility being half its peak value (βϕ = 2.41, 2.68), the biased water number distribution is no longer bimodal.

asymptotic behavior of ⟨Nv⟩ϕ as ϕ → ∞. Equation 6 for N = 1 is ⎡ ⎛ Q ⎞⎤ ∂ ϕ ⎟⎟⎥ = Pv(0) lim [e βϕ⟨Nv⟩ϕ ] Pv(1) = lim ⎢ ⎜⎜ ϕ →∞ z → 0⎢ z Q ∂ ⎣ ⎝ 0 ⎠⎥⎦

Bimodal distributions with large barriers separating the lowand high-N basins would result in broken ergodicity, leading to erroneous estimates of ⟨Nv⟩ϕ and, consequently, of μex v . For v near hydrophilic surfaces, the unbiased water number distribution, Pv(N), is approximately Gaussian.26,52 As a result, the biased water number distributions are also expected to be nearly Gaussian52 or, at the very least, unimodal for the entire range of ϕ values. However, water near an extended hydrophobic surface sits at the edge of a dewetting transition,23,24,26−28,52,53 and this proximity to the dewetting transition is reflected in a fat tail in the unbiased water number distribution, the sensitive response of water density to unfavorable perturbations,52,52 and a corresponding peak in the susceptibility (for the toy-model surfaces, see Figures 1a, 1b, and 6a, respectively). It is important to note that because of the molecular volumes considered, however, the dewetting transition is not a true “thermodynamic” phase transition, so we expect coexistence between the wet and dewet (high and low N, respectively) basins only when v is sufficiently large and is situated near a very hydrophobic surface. Even then, the biased water number distribution is expected to be bimodal only near ϕ = ϕ*, the ϕ value corresponding to the peak in susceptibility. In Figure 6b, we show the biased water number distribution, at ϕ = ϕ*, for the two most hydrophobic toy-model surfaces. The distribution is bimodal for only the most hydrophobic surface (λ = 8 × 10−3), and the barrier separating the high and low N basins is small. This suggests that significantly larger (with more than 100 waters) volumes or more hydrophobic surfaces (or both) would be required to break ergodicity. In Figure 6c, we show the biased water number distribution, near ϕ = ϕ* for the λ = 8 × 10−3 hydrophobic surface. In addition to ϕ = ϕ*, biased distribution functions are shown for ϕ = ϕ1/2, ϕ values corresponding to the susceptibility being half its peak value. For the λ = 8 × 10−3 hydrophobic surface, the biased distribution functions at ϕ = ϕ1/2 are unimodal, confirming that the distributions are bimodal only for a narrow range of ϕ values near ϕ*. Even though biased water number distributions are expected to be unimodal under most circumstances, we recommend that the susceptibility, ⟨δNv2⟩ϕ, be estimated for each biased simulation. While a nearly constant susceptibility indicates an underlying hydrophilic surface and a unimodal biased distribution, a sharp peak in susceptibility suggests an underlying hydrophobic surface with the possibility of a

(11)

implying that as ϕ → ∞, ⟨Nv⟩ϕ → [Pv(1)/Pv(0)] exp(−βϕ). In Figure 5a, we demonstrate the validity of this asymptotic behavior by plotting the logarithm of eβϕ⟨Nv⟩ϕ as a function of βϕ for the toy-model surfaces. In all cases, the curves plateau at large ϕ. Further, in Figure 5b, we show that the plateau values are consistent with eq 11. Assuming that ⟨Nv⟩ϕ follows its asymptotic form for ϕ > ϕmax, the residual integral from ϕmax to ∞ can easily be estimated as

∫ϕ

∞

⟨Nv⟩ϕ dϕ = kBT ⟨Nv⟩ϕmax

max

(12)

Thus, a more practical version of eq 7 is μvex =

∫0

ϕmax

⟨Nv⟩ϕ dϕ + kBT ⟨Nv⟩ϕmax

(13)

Because estimates of ⟨Nv⟩ϕ get progressively noisier for larger values of ϕ as a result of reweighting errors, a judicious choice of ϕmax is one that is large enough that eβϕ⟨Nv⟩ϕ has reached its plateau, but not so large that reweighting results in large errors in ⟨Nv⟩ϕ. A convenient way to ensure that eβϕ⟨Nv⟩ϕ has reached its plateau entails estimating the slope of log⟨Nv⟩ϕ + βϕ, which approaches 0 as ϕ → ∞. Using eq 2, we can show that the slope is given by 1 ∂⟨Nv⟩ϕ ∂ [log⟨Nv⟩ϕ + βϕ] = +1 ∂(βϕ) ⟨Nv⟩ϕ ∂(βϕ) ⎡ ⟨δN 2⟩ ⎤ v ϕ = −⎢ − 1⎥ ⎢⎣ ⟨Nv⟩ϕ ⎥⎦

(14)

For the toy-model system, Figure 5c confirms that the slope goes to 0 for large ϕ. Thus, ϕmax may be chosen to be the smallest ϕ for which ⟨δNv2⟩ϕ/⟨Nv⟩ϕ − 1 is less than some predetermined small (≪1) cutoff. Loss of Ergodicity and Resultant Errors in ⟨Nv⟩ϕ

Here, we discuss when water number distributions in the presence of a linear biasing potential, ϕNv, are bimodal. 1571



Article

(14) Zhou, R.; Huang, X.; Margulis, C. J.; Berne, B. J. Hydrophobic Collapse in Multidomain Protein Folding. Science 2004, 305, 1605− 1609. (15) Liu, P.; Huang, X.; Zhou, R.; Berne, B. J. Observation of a Dewetting Transition in the Collapse of the Melittin Tetramer. Nature 2005, 437, 159−162. (16) Giovambattista, N.; Debenedetti, P. G.; Rossky, P. J. Enhanced Surface Hydrophobicity by Coupling of Surface Polarity and Topography. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15181−15185. (17) Surblys, D.; Yamaguchi, Y.; Kuroda, K.; Nakajima, T.; Fujimura, H. Analysis on Wetting and Local Dynamic Properties of Single Water Droplet on a Polarized Solid Surface: A Molecular Dynamics Study. J. Chem. Phys. 2011, 135, 014703. (18) Rotenberg, B.; Patel, A. J.; Chandler, D. Molecular Explanation for Why Talc Surfaces Can be Both Hydrophilic and Hydrophobic. J. Am. Chem. Soc. 2011, 133, 20521−20527. (19) Kyte, J.; Doolittle, R. F. A Simple Method for Displaying the Hydropathic Character of a Protein. J. Mol. Biol. 1982, 157, 105−132. (20) Black, S. D.; Mould, D. R. Development of Hydrophobicity Parameters To Analyze Proteins which Bear Post- or Cotranslational Modifications. Anal. Biochem. 1991, 193, 72−82. (21) Jamadagni, S. N.; Godawat, R.; Garde, S. Hydrophobicity of Proteins and Interfaces: Insights from Density Fluctuations. Ann. Rev. Chem. Biomol. Eng. 2011, 2, 147−171. (22) Snyder, P.; Lockett, M.; Moustakas, D.; Whitesides, G. Is It the Shape of the Cavity, or the Shape of the Water in the Cavity? Eur. Phys. J. Special Top. 2013, 1−39. (23) Lum, K.; Chandler, D.; Weeks, J. D. Hydrophobicity at Small and Large Length Scales. J. Phys. Chem. B 1999, 103, 4570−4577. (24) Chandler, D. Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640−647. (25) Godawat, R.; Jamadagni, S. N.; Garde, S. Characterizing Hydrophobicity of Interfaces by Using Cavity Formation, Solute Binding, and Water Correlations. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15119−15124. (26) Patel, A. J.; Varilly, P.; Chandler, D. Fluctuations of Water near Extended Hydrophobic and Hydrophilic Surfaces. J. Phys. Chem. B 2010, 114, 1632−1637. (27) Varilly, P.; Patel, A. J.; Chandler, D. An Improved CoarseGrained Model of Solvation and the Hydrophobic Effect. J. Chem. Phys. 2011, 134, 074109. (28) Patel, A. J.; Varilly, P.; Jamadagni, S. N.; Acharya, H.; Garde, S.; Chandler, D. Extended Surfaces Modulate Hydrophobic Interactions of Neighboring Solutes. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 17678− 17683. (29) The hydration free energy for a thin cuboid-shaped cavity with cross-sectional area Ac near a flat surface with contact angle θ is given ex −1 ≡ kBT is the thermal by βμex v = βμbulk − γAc(1 − cos θ), where β energy (kB is the Boltzmann constant and T is the temperature), the subscript v in μex v represents the properties of the volume (size, shape, and location) that must be emptied to form the cavity, μex bulk is the cavity hydration free energy in bulk water, and γ is the vapor−liquid surface tension. (30) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808−2812. (31) Hummer, G.; Garde, S.; Garcia, A. E.; Pohorille, A.; Pratt, L. R. An Information Theory Model of Hydrophobic Interactions. Proc. Nat. Acad. Sci. 1996, 93, 8951−8955. (32) Siebert, X.; Hummer, G. Hydrophobicity Maps of the N-Peptide Coiled Coil of HIV-1 gp41. Biochemistry 2002, 41, 2956−2961. (33) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (34) Shenogina, N.; Godawat, R.; Keblinski, P.; Garde, S. How Wetting and Adhesion Affect Thermal Conductance of a Range of Hydrophobic to Hydrophilic Aqueous Interfaces. Phys. Rev. Lett. 2009, 102, 156101. (35) Hakanpaa, J.; Linder, M.; Popov, A.; Schmidt, A.; Rouvinen, J. Hydrophobin HFBII in Detail: Ultrahigh-Resolution Structure at 0.75 Å. Acta Crystallogr. Sect. D 2006, 62, 356−367.

bimodal biased distribution. Such a possibility should be explored further by performing biased simulations in the forward and backward directions (or alternatively, by starting from a high-N and a low-N state). Hysteresis in the corresponding response functions, ⟨Nv⟩ϕ vs ϕ, indicates broken ergodicity resulting from large barriers separating the low- and high-N basins. The difference in the estimates of μex v from the forward and backward response curves also provides us with an estimate of the error in μex v . If the error is unacceptably large, it may be more efficient to perform umbrella sampling using parabolic biasing potentials with a sufficiently large spring constant so that the biased system commits to either the liquid or the vapor basin and results in a unimodal biased water number distribution.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected] Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank David Chandler and Patrick Varilly for numerous helpful discussions and acknowledge financial support in part from NSF(CBET-1159990 and DMR-1207411).

■

REFERENCES

(1) Granick, S.; Bae, S. C. A Curious Antipathy for Water. Science 2008, 322, 1477−1478. (2) Mobley, D. L.; Dill, K. A. Binding of Small-Molecule Ligands to Proteins: “What You See” Is Not Always “What You Get”. Structure 2009, 17, 489−498. (3) 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. (4) Ito, T.; Chiba, T.; Ozawa, R.; Yoshida, M.; Hattori, M.; Sakaki, Y. A Comprehensive Two-Hybrid Analysis To Explore the Yeast Protein Interactome. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 4569−4574. (5) de Gennes, P. G. Wetting: Statics and Dynamics. Rev. Mod. Phys. 1985, 57, 827−863. (6) Giovambattista, N.; Lopez, C. F.; Rossky, P. J.; Debenedetti, P. G. Hydrophobicity of Protein Surfaces: Separating Geometry from Chemistry. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 2274−2279. (7) Mittal, J.; Hummer, G. Interfacial Thermodynamics of Confined Water Near Molecularly Rough Surfaces. Faraday Discuss. 2010, 146, 341−352. (8) Daub, C. D.; Wang, J.; Kudesia, S.; Bratko, D.; Luzar, A. The Influence of Molecular-Scale Roughness on the Surface Spreading of an Aqueous Nanodrop. Faraday Discuss. 2010, 146, 67−77. (9) Giovambattista, N.; Debenedetti, P. G.; Rossky, P. J. Hydration Behavior Under Confinement by Nanoscale Surfaces with Patterned Hydrophobicity and Hydrophilicity. J. Phys. Chem. C 2007, 111, 1323−1332. (10) Hua, L.; Zangi, R.; Berne, B. J. Hydrophobic Interactions and Dewetting between Plates with Hydrophobic and Hydrophilic Domains. J. Phys. Chem. C 2009, 113, 5244−5253. (11) Argyris, D.; Cole, D. R.; Striolo, A. Hydration Structure on Crystalline Silica Substrates. Langmuir 2009, 25, 8025−8035. (12) Acharya, H.; Vembanur, S.; Jamadagni, S. N.; Garde, S. Mapping Hydrophobicity at the Nanoscale: Applications to Heterogeneous Surfaces and Proteins. Faraday Discuss. 2010, 146, 353−365. (13) Wang, J.; Bratko, D.; Luzar, A. Probing Surface Tension Additivity on Chemically Heterogeneous Surfaces by a Molecular Approach. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 6374−6379. 1572



Article

(36) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 435−447. (37) Patel, A. J.; Varilly, P.; Chandler, D.; Garde, S. Quantifying Density Fluctuations in Volumes of All Shapes and Sizes Using Indirect Umbrella Sampling. J. Stat. Phys. 2011, 145, 265−275. (38) Frenkel, D.; Smit, B. Understanding Molecular Simulations: From Algorithms to Applications, 2nd ed.; Academic Press: New York, 2002. (39) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (40) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (41) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (42) Miller, T.; Vanden-Eijnden, E.; Chandler, D. Solvent CoarseGraining and the String Method Applied to the Hydrophobic Collapse of a Hydrated Chain. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14559− 14564. (43) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (44) Chen, Q.; Yan, J.; Zhang, J.; Bae, S. C.; Granick, S. Janus and Multiblock Colloidal Particles. Langmuir 2012, 28, 13555−13561. (45) Kumar, A.; Park, B. J.; Tu, F.; Lee, D. Amphiphilic Janus Particles at Fluid Interfaces. Soft Matter 2013, 9, 6604−6617. (46) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (47) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (48) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (49) Huang, D. M.; Chandler, D. The Hydrophobic Effect and the Influence of Solute−Solvent Attractions. J. Phys. Chem. B 2002, 106, 2047−2053. (50) Chandler, D.; Weeks, J. D.; Andersen, H. C. Van der Waals Picture of Liquids, Solids, and Phase Transformations. Science 1983, 220, 787−794. (51) Willard, A. P.; Chandler, D. Instantaneous Liquid Interfaces. J. Phys. Chem. B 2010, 114, 1954−1958. (52) Patel, A. J.; Varilly, P.; Jamadagni, S. N.; Hagan, M. F.; Chandler, D.; Garde, S. Sitting at the Edge: How Biomolecules Use Hydrophobicity to Tune Their Interactions and Function. J. Phys. Chem. B 2012, 116, 2498−2503. (53) Chandler, D.; Varilly, P. Lectures on Molecular- and Nano-Scale Fluctuations in Water; Cornell University Library; DOI: arXiv 1101.2235 [cond-mat.soft], 2011. (54) Weeks, J.; Chandler, D.; Andersen, H. Role of Repulsive Forces in Forming the Equilibrium Structure of Simple Liquids. J. Chem. Phys. 1971, 54, 5237−5247. (55) Perchiacca, J. M.; Tessier, P. M. Engineering AggregationResistant Antibodies. Annu. Rev. Chem. Biomol. Eng. 2012, 3, 263−286. (56) White, A.; Jiang, S. Local and Bulk Hydration of Zwitterionic Glycine and its Analogues through Molecular Simulations. J. Phys. Chem. B 2011, 115, 660−667. (57) Freed, A. S.; Garde, S.; Cramer, S. M. Molecular Simulations of Multimodal Ligand−Protein Binding: Elucidation of Binding Sites and Correlation with Experiments. J. Phys. Chem. B 2011, 115, 13320− 13327. (58) Snyder, P. W.; Mecinović, J.; Moustakas, D. T.; Thomas, S. W.; Harder, M.; Mack, E. T.; Lockett, M. R.; Héroux, A.; Sherman, W.; Whitesides, G. M. Mechanism of the Hydrophobic Effect in the

Biomolecular Recognition of Arylsulfonamides by Carbonic Anhydrase. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 17889−17894. (59) Ma, H.; Horiuchi, K. Y.; Wang, Y.; Kucharewicz, S. A.; Diamond, S. L. Nanoliter Homogenous Ultra-High Throughput Screening Microarray for Lead Discoveries and IC50 Profiling. Assay Drug Dev. Technol. 2005, 3, 177−187. (60) Li, I. T.; Walker, G. C. Signature of Hydrophobic Hydration in a Single Polymer. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 16527−16532. (61) Garde, S.; Patel, A. J. Unraveling the Hydrophobic Effect, One Molecule at a Time. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 16491− 16492. (62) Nucci, N. V.; Pometun, M. S.; Wand, A. J. Site-Resolved Measurement of Water-Protein Interactions by Solution NMR. Nat. Struct. Mol. Biol. 2011, 18, 245−249. (63) Nucci, N. V.; Pometun, M. S.; Wand, A. J. Mapping the Hydration Dynamics of Ubiquitin. J. Am. Chem. Soc. 2011, 133, 12326−12329. (64) Franck, J. M.; Scott, J. A.; Han, S. Nonlinear Scaling of Surface Water Diffusion with Bulk Water Viscosity of Crowded Solutions. J. Am. Chem. Soc. 2013, 135, 4175−4178. (65) King, J. T.; Kubarych, K. J. Site-Specific Coupling of Hydration Water and Protein Flexibility Studied in Solution with Ultrafast 2D-IR Spectroscopy. J. Am. Chem. Soc. 2012, 134, 18705−18712. (66) Mattos, C.; Bellamacina, C. R.; Peisach, E.; Pereira, A.; Vitkup, D.; Petsko, G. A.; Ringe, D. Multiple Solvent Crystal Structures: Probing Binding Sites, Plasticity and Hydration. J. Mol. Biol. 2006, 357, 1471−1482. (67) Pohorille, A.; Jarzynski, C.; Chipot, C. Good Practices in FreeEnergy Calculations. J. Phys. Chem. B 2010, 114, 10235−10253.

1573


Efficient Method To Characterize the Context-Dependent

Recommend Documents