New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Dryinginduced hydrophobic polymer collapse

Contributed by David Chandler
Abstract
We have used computer simulation to study the collapse of a hydrophobic chain in water. We find that the mechanism of collapse is much like that of a firstorder phase transition. The evaporation of water in the vicinity of the polymer provides the driving force for collapse, and the rate limiting step is the nucleation of a sufficiently large vapor bubble. The study is made possible through the application of transition path sampling and a coarsegrained treatment of liquid water. Relevance of our findings to understanding the folding and assembly of proteins is discussed.
For nearly a halfcentury, hydrophobic interactions have been considered the primary cause for self assembly in soft matter, and a major source of stability in biophysical assembly (1, 2). Studying these interactions in perhaps their most basic form, we use computer simulations to demonstrate the mechanism for the collapse of a hydrophobic polymer in water. We show that solvent fluctuations induce the transition from the extended coil to the collapsed globule state, where a vapor bubble of sufficient size is formed to envelop and thereby stabilize a critical nucleus of hydrophobic units. This mechanism is different from that usually considered, where coil to globule transitions are attributed to effective interactions between pairs of chain segments and a change in sign of second virial coefficient (3). Rather, the mechanism we find is evocative of the ncluster model, where hydrophobic collapse is produced by solventinduced interactions between a relatively large cluster of segments (4).
As expected from earlier work on the equilibrium theory of hydrophobicity (5), we find that the solvent length scales pertinent to hydrophobic collapse extend over nanometers. We also find that pertinent time scales extend beyond nanoseconds. Given these molecularly large lengths and times, it is understandable that no work before this has provided statistically meaningful computer simulations of the process. Our use of a statistical field model of water allows us to simulate solvent dynamics over large length and time scales that would be impractical to study with purely atomistic simulation. Spatially complex small lengthscale fluctuations are analytically integrated out, thus removing the most computationally costly features from our simulation. Their integration can be performed at the outset because their relaxation is relatively fast (6) and their statistics is essentially Gaussian (7). Only the polymer degrees of freedom and a coarsegrained density field remain. The equilibrium theory for this approach has been detailed in an earlier paper (8). Here, we use a version of the model that is suitably generalized for dynamical applications.
By “small length” we refer to distances smaller than l ≈ 0.3 nm. In the absence of any strong perturbation, such as those that can occur close to a solute, these small lengthscale fluctuations are the only fluctuations of significance. Larger lengthscale fluctuations are generally insignificant in water at ambient conditions. The liquid is relatively cold and incompressible, and spatial correlations extend over only one or two water molecules. But the presence of a hydrophobic polymer can change this situation, making large lengthscale fluctuations important. The liquid lies close to macroscopic vapor–liquid equilibrium, and vaporlike behavior is stabilized in the vicinity of a sufficiently extended hydrophobic surface (9)—a surface formed, for example, by the clustering of hydrophobic groups in a polymer chain. This effect can be captured by the behavior of a coarsegrained density field (5, 10), and such a field is conveniently simulated with a field of binary numbers (8), as we discuss now.
Equilibrium of the Solvent Model
Dividing space into a cubic grid of cells, each with side length l, this binary field is a dynamical variable given by ρ_{l}n_{i}, where ρ_{l} is a constant given by the bulk liquid density, and n_{i} is dynamical given by n_{i} = 1 if cell i contains “liquid,” and n_{i} = 0 if it contains “vapor.” This coarsegrained density is thus the “Ising” field that describes phase coexistence and gas–liquid interfaces (11). The remaining small lengthscale field, δρ(r), is the difference between ρ_{l}n_{i} and the actual density at a point r in cell i. It takes on a continuum of values, obeys Gaussian statistics, and captures the granularity of the solvent (12).
Confining our attention to an ideal hydrophobic solute, namely a polymer consisting of hard spheres, the integration over δρ(r) yields the following Hamiltonian for the coarsegrained field in the grandcanonical ensemble (8): 1 The quantity μ is the solvent chemical potential, ɛ is the nearestneighbor coupling energy parameter of the solvent, the sum over 〈i,j〉 includes all nearest neighbor cells, v_{i} denotes the volume excluded by the solute in cell i, and 2 with c = 2.67 × 10^{8} J/m^{3} (about 60 k_{B}T/nm^{3} at room temperature), is the reversible work to accommodate that excluded volume in the liquid. Eqs. 1 and 2 are approximations to the full results derived in ref. 8. Specifically, in comparison with Equations 19 or 23 of ref. 8, we see that Eq. 1 neglects relatively small intercell interactions. Further, Eq. 2 is a reasonable numerical approximation to the general result for Δμ_{ex}(v_{i}) that follows from Gaussian statistics for δρ(r) (7). The specific value of c is chosen so that Δμ_{ex}(4πR^{3}/3) yields the excess chemical potential for a hard sphere or bubble that excludes solvent from a spherical volume of radius R = 0.5 nm (13).
The predictions of the model for excess chemical potential or solvation free energy of a hard sphere over the range of all possible radii is shown in Fig. 1. It is seen that the agreement between the results of the model and those of an atomistic simulation is reasonable, because the model captures the appropriate freeenergy scaling at small and large lengths and gives a reasonable prediction of the crossover length between the two regimes. For closer agreement between the statistical field model and the atomistic model in the crossover regime, terms neglected in Eq. 1 should be included, as demonstrated in ref. 8.
The crossover occurring around R ∼ 1 nm is a collective effect. It is induced by the entropic cost of constraining small lengthscale fluctuations near a hydrophobic solute. For a given cell i, the cost is Δμ_{ex}(v_{i}). When the v_{i}'s form a connected volume of large enough extent, the net entropic cost is larger than the energetic cost for forming an interface that surrounds the volume. At this point, and for larger volumes, it becomes energetically preferable for n_{i} = 0 for all cells i occupied by or adjacent to the large solute, and free energetics is then dominated by the energy of the interface thus formed. Such interfaces are stable only in the neighborhood of vapor–liquid phase equilibrium. Therefore, crossover behavior can occur only when μ − μ_{coex}≪ k_{B}T, where k_{B}T is Boltzmann's constant times temperature (the thermal energy) and μ_{coex} is the chemical potential at phase coexistence between liquid and vapor. Further, because a pressure–volume term generally contributes to the free energy of a large void, crossover behavior occurs only when pressure is low enough to make this term similarly small. We will see below that the crossover in the scaling behavior of the free energy is closely related to the nucleation of oil–water demixing and the mechanism for the collapse of a hydrophobic chain.
Dynamics of the Solvent and Chain
The excluded volumes v_{i} are dynamical, changing as the chain configuration changes with time. The chain need not move slowly on the time scale for which disturbances in the solvent may relax. As such, reversible work surfaces, like the solvation free energy graphed in Fig. 1, are not entirely pertinent, and the dynamics of the liquid and chain should be studied together. For this purpose, we have constructed a stochastic dynamics that satisfies the basic requirements of time reversibility and the preservation of an equilibrium distribution. With this dynamics, solute particles move through space with continuous variation of coordinates while the coarse grained density evolves simultaneously through its discrete configurations. Because δρ(r) has been integrated out (its effects appear only implicitly), the trajectories can be true to nature only for times greater than those required to relax δρ(r). This time is on the order of 1 ps (6). On this time scale, the appropriate stochastic dynamics is noninertial.
Trajectories of our model are advanced as follows: The polymer, composed of N_{s} = 12 spheres is moved for M_{s} time steps from a configuration {r_{α}} to a new configuration {r} according to propagation rules of Langevin dynamics, in the field of constant coarse grained density variables, {n_{i}}. At the end of those steps, the polymer is held fixed and the coarse grained density is moved to a new configuration {n} by carrying out a full sweep through all cells, applying Glauber Monte Carlo dynamics (17) to each cell with the Hamiltonian H[{n_{i}},{v_{i}}]. The net procedure taking the system from {r_{α}, n_{i}} to {r, n} is then repeated over and over again. The trajectories thus formed are reversible, preserve the norm of the configurational distribution function, and obey detailed balance.
For one of the M_{s} time steps between solvent sweeps, carrying the monomers from their configuration at time t to that at time t + δt_{s}, the Langevin propagation corresponds to 3 where F_{α}({r_{ξ}(t)}, {n_{i}}) is the force acting on monomer α, and γ is the friction constant. The force contains a random part, δF, that is independent of configurational variables. It is the dynamical remnant of the small lengthscale field, δρ(r). We take it to be isotropic, to have a correlation time of zero, and to obey Gaussian statistics with zero mean and variance 〈δF^{2}〉 = 3k_{B}Tγ. We estimate the value of γ by associating the diffusion constant of a single bead moving in water, D_{s} = k_{B}T/γ, to its Stokes estimate which gives γ = 8.39 × 10^{−12} kg/s. The remaining contribution to the force is not random, and is given by −∇_{α}[V({r_{ξ}}) + c∑_{i}v_{i}]. The time step is δt_{s} = 1.4 × 10^{−14} s, which is small enough that Eq. 3 preserves the canonical distribution.
A physically meaningful value of M_{s} coincides with M_{s} ≈ δt_{l}/δt_{s}, where δt_{l} is the physical time associated with a sweep through the coarsegrained density of the solvent. This solvent time should coincide with the correlation time for a density fluctuation of lengthscale l—i.e., δt_{l} = 1/[D(2π/l)^{2})], where D is the selfdiffusion constant of liquid water. This relationship gives δt_{l} = 5.0 × 10^{−13} s, and thus M_{s} = 36.
With these dynamical rules, we show below that the halflife of the extended chain we consider is approximately 10^{−5} s at room temperature, whereas the transient time to move off a threshold between extended coil and compact globule states is of the order of 10^{−10} s. Thus, collapse transitions in this system are rare events (16). Harvesting a representative ensemble of these events without prior knowledge of transition states or mechanism is accomplished with transition path sampling (17). To apply this technique, we use the number of vapor cells, Σ_{i}(1 − n_{i}), as the order parameter to distinguish the solvated coil and globule states.
Results
Fig. 2 illustrates a 1.5ns trajectory in which a chain of 12 hydrophobic units collapses from an extended coil to a compact globule. The volume occupied by each separate unit is typical of that for amino acids in solution. In vacuum, this particular chain will always remain in a coil state because the global minimum in its intrachain potential energy function corresponds to the fully extended conformation. Nevertheless, the collapsed chain is very much the stable state of the solvated chain, as is evident from the freeenergy surface shown in Fig. 3. Although rare, the trajectory illustrated in Fig. 2 is a representative example of the subensemble of all 1.5ns trajectories that do exhibit the collapse transition at temperature T ≈ 300 K. The polymer is initially in the extended coil state, after which a spontaneous fluctuation in chain conformation collapses a section of the polymer. The collapsed section forms a sufficiently large hydrophobic cluster that the formation of a vapor bubble is induced. Eventually, this vapor bubble grows and drives all of the units of the polymer together.
Visual inspections of this and similar dynamical trajectories suggest that the mechanism of the collapse transition arises from an interplay between the size of the polymer and the formation of a vapor bubble. For this reason, we have mapped out the freeenergy landscape in terms of variables manifesting the polymer size and the bubble size. The first of these variables, R, is the squared radius of gyration of the polymer. The other, U, is the volume of the largest vapor bubble, in units of l^{3}. To identify the bubbles in the system, we use the criterion that any two vapor cells (i.e., cells for which n_{i} = 0) belong to the same bubble whenever they are nearest neighbor cells.
The contour plot for the room temperature freeenergy landscape (Fig. 3) shows that the path of lowest free energy from the coil to the collapsed globule is one where, initially, the radius of gyration decreases, while the size of the largest bubble is still essentially zero. In this regime, the solvent still wets the polymer (i.e., n_{i} is mostly 1 for cells occupied by the solute), and the free energy hardly changes. When the radius of gyration becomes small enough, however, a bubble starts to grow. Here, the free energy increases sharply by about 9 k_{B}T where it reaches a saddle point at (U, R) = (98, 23.5 l^{2}). Beyond that saddle point, the bubble grows spontaneously to a size that eventually envelopes the fully collapsed globule. The height of the barrier only weakly depends on chain stiffness because it is due primarily to solvent reorganization.
A close examination of the freeenergy landscape reveals a small barrier of height 2 k_{B}T at (U, R) = (6.5, 70 l^{2}). This feature, not visible on the scale of the graph plotted in Fig. 3, separates the coil from a more compact metastable intermediate. The presence of these two states arises from a competition between the entropy of the chain, which favors the coil state, and weak depletion forces, which favor a more compact state. Weak depletion forces are caused by the reduction in the volume from which the solutes exclude the solvent, when solutes come together while still wet. The attraction between two small hydrophobic objects predominantly arises from this effect (7, 18, 19), although its full description (with the characteristic oscillations in potentials of mean force) requires a more accurate evaluation of Δμ_{ex} (v_{i}) than Eq. 2. Whether evaluated to high accuracy or approximately, this driving force is relatively small, and it is not responsible for spontaneous assembly of hydrophobic units (20). The freeenergy difference is 30 k_{B}T between the coil and the fully collapsed globule, whereas it is only a few k_{B}T between the coil and the intermediate state. The strong driving force for the collapse of the polymer comes from the emptying of cells (i.e., solvent cavitation or drying) and thus the demixing of the hydrophobic units and water. The resulting stabilization follows from the fact that the solvent is close to liquid–vapor equilibrium so that the large lengthscale interfacial free energy of the cavity is far lower than the small lengthscale entropic cost of maintaining a wet state, Σ_{i} Δμ_{ex} (v_{i}). Thus, nucleating the collapsed hydrophobic chain involves the same physical effect that is responsible for the crossover phenomenon discussed above in reference to Fig. 1.
To judge the extent to which R and U describe the pertinent dynamics of the collapse, and are thus good reaction coordinates, we have performed extensive transition path sampling of an ensemble of trajectories, each of length 150 ps. Points taken from these transition trajectories are shown in juxtaposition with the freeenergy surface in Fig. 3. Although ten times shorter than the one illustrated in Fig. 2, the 150ps trajectories are sufficiently long to capture the mechanism of the collapse because 150 ps is significantly longer than the time for a trajectory starting at the dynamical bottleneck to commit to a basin of attraction. Fig. 3 shows that the flow through the bottleneck in this system is primarily due to the solvent, and the solvent in this transition state regime moves on the time scale of picoseconds. To identify the bottleneck (i.e., transitionstate surface), we have located the configurations on each trajectory from which newly initiated trajectories have equal probability of ending in the coil or globule states. These points in configuration space are members of the transition state ensemble (17, 21, 22). We project them onto Fig. 3, where it is seen that points in the transition state ensemble are indeed close to the saddle point in the freeenergy surface. Thus, U and R are the predominant reaction coordinates for this system. The transitionstate ensemble is slightly tilted in the (U, R)plane, showing that the larger the polymer, as measured by its radius of gyration, the larger the size of the critical vapor bubble. However, the scatter of the transitionstate ensemble from a line in that plane is notable, indicating that at least one other variable in addition to U and R plays a quantitative role in the reaction coordinate.
The transition paths shown in Fig. 3 differ from the lowest freeenergy path, and the differences are generally larger than k_{B}T. This behavior shows that the polymer and the solvent move on significantly different time scales when passing through and moving beyond the dynamical bottleneck. After a vapor bubble is nucleated, the polymer does not respond on the time scale at which the bubble grows. In fact, by taking an artificially large value for M_{s} (3,600) to accelerate the dynamics of the chain relative to that of the solvent by a factor of 100, the projected dynamical paths closely follow the lowest freeenergy path.
Discussion
Trajectories of the hydrophobic collapse are generally parallel to U as they pass over the transition state. This finding demonstrates unambiguously that the dynamics of collapse of a hydrophobic polymer in water is dominated by the dynamics of water, specifically the collective emptying of regions of space near a nucleating cluster of hydrophobic species. Further, both the solvent and chain remain out of equilibrium throughout the collapse transition. Real hydrophobic molecules have some affinity to water, and vice versa, because of ubiquitous van der Waals interactions. The ideal hydrophobic chain considered here, however, has no such property because it is composed simply of hard spheres. Nevertheless, we expect our conclusions to remain valid for more realistic models of hydrophobic species because the underlying mechanism of nucleation is largely unaffected by van der Waals interactions (23). Our findings would seem also pertinent to the mechanism of biological assembly, such as protein folding, but to demonstrate so with simulation will require an analogous simulation study of a proteinlike chain. In the past, a statistically meaningful study of such folding with explicit account of solvent dynamics seemed impractical. The current analysis, however, suggests that such studies will become feasible by using a statistical field model like that used herein.
Appropriate models for hydrophobic collapse must account for both the physics of small lengthscale fluctuations and the physics of phase equilibria and interfaces. This accounting is difficult to accomplish without an explicit treatment of the solvent, as we do here. Often, solvation is treated with implicit solvent models in which all solvent degrees of freedom are averaged out in some approximate fashion, assuming that all of these degrees of freedom relax quickly. Some of these models mimic the effects of the solvent through estimates of the pair potential of mean force between solute units. This approach can be correct for describing small lengthscale effects, and therefore can be correct in situations where the solute remains wet. But for collections of hydrophobic solutes large enough to induce drying or cavitation, an implicit model would need to describe correlations of essentially all orders, and not simply pair correlations. For this reason, a recent study of hydrophobic effects in protein folding requires further analysis (24). Other implicit solvent models are based on estimates of exposed surface area (25, 26). This approach can be valid for large enough solutes, provided the time scales of interest are longer than nucleation times (10^{−5} s in the system studied here). But applying it to small unassembled solutes is incorrect because the free energy in this regime scales with excluded volume and not surface area.
The connection we stress between hydrophobic collapse and phase equilibria provides a simple explanation for both cold denaturation and pressure denaturation of proteins. Namely, the tendency for drying and its concomitant stabilization of hydrophobic clustering decreases with decreasing proximity of liquid–vapor equilibrium. Thus, the lowering temperature and the increasing of pressure destabilize hydrophobic collapse because both these actions move ambient water away from its phase equilibrium with vapor. The proximity of phase equilibria is also of relevance to freeenergy barriers for hydrophobic collapse, like that illustrated in Fig. 3. These barriers will increase, for example, with increasing external pressure. Thus, the most direct measure of the importance of hydrophobic collapse in protein folding may therefore come from studying the kinetic effects of changing pressure.
Acknowledgments
We are grateful to Aaron Dinner, David Huang, and Shura Grosberg for helpful comments on an earlier draft of this paper. This work was supported in its initial stages by National Science Foundation Grants 9508336 and 0078458, and in its final stages by the Director (Office of Science, Office of Basic Energy Sciences) of the U.S. Department of Energy (Grant DEAC0376SF00098.
 Accepted March 14, 2002.
 Copyright © 2002, The National Academy of Sciences
References
 ↵
 ↵
 Tanford C
 ↵
 Yu Grosberg A,
 Khokhlov A R
 ↵
 DeGennes PG
 ↵
 ↵
 ↵
 Hummer G,
 Garde S,
 García A E,
 Pohorille A,
 Pratt L R
 ↵
 ↵
 Stillinger F H
 ↵
 ↵
 Chandler D
 ↵
 ↵
 Huang D M,
 Geissler P L,
 Chandler D
 ↵
 Frenkel D,
 Smit B

 Landau D P,
 Binder K
 ↵
 ↵
Bolhuis, P., Chandler, D., Dellago, C. & Geissler, P. L. (2001) Annu. Rev. Phys. Chem., in press.
 ↵
 ↵
 ↵
 Watanabe K,
 Andersen H C
 ↵
 ↵
 Northrup S,
 Pear M R,
 Lee CY,
 McCammon A,
 Karplus M
 ↵
 Huang D M,
 Chandler D
 ↵
 Cheung M S,
 Garcia A E,
 Onuchic J N
 ↵
 ↵
 Vallone B,
 Miele A E,
 Vecchini P,
 Chiancone E,
 Brunori M
 ↵