Explicit and implicit modeling of nanobubbles in hydrophobic confinement

Water at normal conditions is a fluid thermodynamically close to the liquid-vapor phase coexistence and features a large surface tension. This combination can lead to interesting capillary phenomena on microscopic scales. Explicitwater molecular dynamics (MD) computer simulations of hydrophobic solutes, for instance, give evidence of capillary evaporation on nanometer scales, i.e., the formation of nanometer-sized vapor bubbles (nanobubbles) between confining hydrophobic surfaces. This phenomenon has been exemplified for solutes with varying complexity, e.g., paraffin plates, coarse-grained homopolymers, biological and solid-state channels, and atomistically resolved proteins. It has been argued that nanobubbles strongly impact interactions in nanofluidic devices, translocation processes, and even in protein stability, function, and folding. As large-scale MD simulations are computationally expensive, the efficient multiscale modeling of nanobubbles and the prediction of their stability poses a formidable task to the ‘nanophysical’ community. Recently, we have presented a conceptually novel and versatile implicit solvent model, namely, the variational implicit solvent model (VISM), which is based on a geometric energy functional. As reviewed here, first solvation studies of simple hydrophobic solutes using VISM coupled with the numerical level-set scheme show promising results, and, in particular, capture nanobubble formation and its subtle competition to local energetic potentials in hydrophobic confinement.


INTRODUCTION
The modeling and description of aqueous properties such as water structure, dynamics, and eventually thermodynamics are obviously of fundamental interest as water is the most abundant fluid on our planet, and governs biological evolution and geomechanical and atmospheric phenomena (Ball 1999).Particularly, on microscopic scales, i.e., on length scales on the order of the size of a water molecule (∼ 3Å) to the size of a hydrogen-bonded network ( 1-100 nm), the structural properties of water in bulk and confinement are crucial for the understanding of micro-to macroscale hierarchical processes in our environment.These small length scales, Selected paper presented at the IUTAM Symposium on Swelling and Shrinking of Porous Materials: From Colloid Science to Poromechanics -August 06-10 2007, LNCC/MCT.E-mail: jdzubiel@ph.tum.dehowever, are still difficult to access directly by experiments, and their exploration by theoretical and computational means has become an important and necessary branch of theoretical physical chemistry and biology in the last few decades.
The theoretical modeling of water can be performed by explicitly resolving its atomic and molecular degrees of freedom by quantum-mechanical (QM) methods (Jensen 2006) or classical molecular dynamics (MD) computer simulations (Allen andTildesley 1987, Frenkel andSmit 1996).In the following, we refer to this as an explicit modeling of water.Due to the high computational cost, QM and MD methods are restricted to small systems ranging from ∼ 100 water molecules in QM methods to ∼ 10 5 molecules in MD simulations.Although the latter number seems to be relatively large, the accessible time scales for these system sizes are typically too small to guarantee ergodicity and good sampling of the statistical ensemble.A typical procedure to reduce computational costs is to integrate out the solvent degrees of freedom by introducing effective interactions ("coarse-graining"), see for instance information theory (Hummer et al. 1996), LCW theory (Lum et al. 1999), or the string method (Miller III et al. 2007).
An extremely coarse but efficient approach is to model the water in a continuum manner, i.e., by describing water properties by macroscopic constants only, such as the surface tension and dielectric constant.Those quantities are then assumed to be locally defined in space depending on the particular microscopic aqueous environment.In theoretical biochemistry, this is a common approach to obtain quick estimates of solvation or binding free energies of proteins.(Note that a protein, or in general a solute, can be considered an external, confining potential.)The approach is mainly based on the so-called solvent-accessible surface (SAS) and Poisson-Boltzmann (PB) electrostatics (Roux and Simonson 1999).The SAS, which has to be defined before evaluating any energies, usually serves also as a dielectric boundary.In the following we refer to the continuum modeling of water as an implicit modeling.
Water at normal conditions (i.e., pressure P = 1 bar and temperature T = 300 K) is a fluid thermodynamically close to the liquid-vapor phase transition and features a relatively high surface tension.As a consequence, in strong hydrophobic confinement, capillary evaporation can be triggered -a well known phenomenon in the physics of phase transitions (Kralchevsky and Nagayama 2001).Spelling it out, a stable vapor bubble can form between the confining surfaces (Chandler 2005, Rasaiah et al. 2008).The physical reason behind that phenomenon is that the water can minimize its free energy by evaporating and reducing unfavorable liquid-solid interface area in the hydrophobic environment.As the surface tension of water is large, the thermodynamic volume (V ) work P V for evaporating plays only a minor role on microscopic (∼ nm) scales.Sometimes in the literature the vapor bubble on these scales is called a 'nanobubble'.The nanobubbles in this work however, must not be confused with aqueous surface bub-bles (Attard 2000, Parker et al. 1994) which are due to fluctuations and impurities close to the solid surfaces.These fluctuations, however, may be the trigger to capillary evaporation which usually comes with a nucleation barrier (Huang et al. 2003).
Figure 1 illustrates a simple plate-like confinement with two hydrophobic square plates (length L) in a distance D on a nanometer scale.The free energy difference between the filled state (liquid between the two plates) and the vapor state (nanobubble) can be estimated by simple macroscopic arguments and is where P is the liquid bulk pressure and the vapor pressure is assumed to be zero.γ is the surface tension that is assumed to be the same for all interfaces and curvature effects are neglected.The first term is the thermodynamic work to create a vapor hole against the liquid pressure, while the second term is the interfacial work gained by removing the two plate-liquid interfaces, and the third term the interfacial work paid by forming the four liquid-vapor interfaces.On these scales, the pressure term can be neglected as the bubble volume is small and the interfacial terms dominate.We obtain G 2γ L(2D − L), showing that, for distances D L/2, G < 0, the nanobubble is thermodynamically stable and should have long life times.Given the large surface tension of water, we find for the nanometer plates a large mutual hydrophobic attraction of G γ nm 2 20k B T in agreement with atomistic computer simulations (Koishi et al. 2004).Also see Figure 4 later in this work.
The phenomenon of nanobubble formation in hydrophobic confinement has been confirmed in the last decade by a large number of explicit computer simulations of, e.g., plate-like solutes (Chandler 2005), homopolymers (ten Wolde and Chandler 2002), or channels and pores (Beckstein et al. 2001, Dzubiella and Hansen 2004b, 2005, Rasaiah et al. 2008, Vaitheesvaran et al. 2004).It has been suggested that nanobubbles can play a crucial role in protein stability, folding, and function.Once a nanobubble is formed, it leads to strong and long-ranged hydrophobic forces among the "dry" regions as the system further tries to minimize unfavorable interface area (see for instance the simple plate model above) and can induce quick protein folding and collapse (ten Wolde and Chandler 2002) or stable protein assemblies (Liu et al. 2005).Importantly, the evaporation of liquid in the narrow hydrophobic pore region of ion channel proteins (Anishkin and Sukharev 2004, Beckstein et al. 2001, Sotomayor et al. 2007) seems to be important for the control of physiological ion conductance.The presence of a vapor bubble blocks ion permeation as the ion would have to desolvate to travel through the vapor.It has been recognized that the bubble stability is sensitive to local geometry and electrostatics, implying that it can play a key role in the voltage or mechanical gating of the channel protein, i.e., the electrostatic or geometrical control of particle (ion) permeability (Beckstein et al. 2001, Dzubiella and Hansen 2004b, 2005, Sotomayor et al. 2007).Related to this, nanobubble formation has been made responsible experimentally to block (bio)polymer translocation through narrow hydrophobic solid-state nanopores (Smeets et al. 2006).Direct experimental evidence of nanobubbles is still elusive as they are small, soft, and transient.For this reason, the explicit and implicit theoretical modeling of nanobubbles is crucial for understanding their impact on molecular processes in aqueous solution.It has been systematically demonstrated that the shape and stability of the vapor bubble sensitively competes with local solute geometry, dispersion, or electrostatic potentials (Dzubiella and Hansen 2003, 2004a, b, 2005, Liu et al. 2005, Vaitheesvaran et al. 2004, Zhou et al. 2004).Figure 1 illustrates an example where the charging of hydrophobic plates can destroy the nanobubble by dragging the polar solvent into the region between the two plates.Clearly, the stability depends on the amount of charge and geometry of the solute that determine the local electrostatic field distribution.Modeling and prediction of the nanobubble transition, shape, and stability with coarse-grained models thus poses a formidable challenge to the 'nanophysical' and mathematical communities.Established implicit solvent models, such as those based on SAS definitions (Roux and Simonson 1999) are not suitable, as here the solutesolvent interface is input, i.e., is it predefined.The stability of a nanobubble, however, is a priori not known, and a model is required that predicts nanobubble stability, or, in other words, the solute-solvent interface location.
In this paper we review selected examples in the explicit and implicit modeling of nanobubbles in hydrophobic confinement and present a method -the variational implicit solvent model (VISM) -which is potentially capable of predicting nanobubbles in arbitrary confinement without any a priori assumption of the interface location.It is the first implicit water model that can predict nanobubble formation as it is based on a geometric minimization procedure and does not need the interface location in the confining system as input, as used by conventional methods.In the next section, the mechanisms of and competition between hydrophobic (interfacial) effects and electrostatics are readily exemplified by explicit-water MD simulations of two generic model systems.After that, the VISM is introduced, along with the level-set method, helpful for the numerical evaluation of the VISM equation.Some simple application are discussed thereafter in the last but one section.Parts of this paper have been published elsewhere (Cheng et al. 2007, Dzubiella and Hansen 2003, 2004a, b, 2005, Dzubiella et al. 2006a, b).

EXPLICIT MODELING OF NANOBUBBLES: SELECTED EXAMPLES
A. TWO SPHERICAL HYDROPHOBIC NANOSOLUTES In Figure 2 we show the (cylindrically averaged) density profiles of water around two hydrophobic nanometer sized solutes from MD simulation at normal condi- tions (Dzubiella andHansen 2003, 2004a).The watersphere interaction is hard-core like and the spheres have a radius of R 1 nm and are kept at a fixed surfaceto-surface distance of 4 Å.Additionally, the spheres are homogeneously and oppositely charged with a) z = 0, b) 2, and c) 5 elementary charges.As can be seen in Figure 2 a), a stable nanobubble is developed between the two neutral spheres due to capillary evaporation.The solute-solvent interface has a catenoidal shape and minimizes the area of the enveloping surface, which is smaller than the total solvent-accessible surface area of the two spheres (8π R 2 ).Upon charging, the spheres become increasingly solvated and the nanobubble becomes unstable and vanishes [Fig. 2 b) and c)].Charg-ing the spheres, i.e., adding attraction between the solute and the polar solvent, competes with the hydrophobic effect and clearly suppresses nanobubble formation due to a favorable solvation.
The nanobubble has a strong impact on the mean force (MF) between the two spheres, as shown in Figure 2 in the right plot.The MF displays a significant attraction for the neutral spheres while the attraction decreases for increasing charge.Without charges, the system is strongly hydrophobic and attempts to further minimize solute-solvent interface area by pushing the spheres together.Charging, on the other hand, attracts solvent to the solute surface and the solutes like to stay in solution.This simple system already exemplifies the subtle competition between hydrophobicity and electrostatics in molecular systems, which might play a role in macromolecular assemblies on nanometer scales.

B. NANOPORES AND ION CHANNELS
Another relevant generic model system for studying hydrophobic effects is that of a cylindrical pore with length L and radius R connecting two reservoirs of water and ions, as shown in Figure 3.In explicit-water MD simulations (Allen et al. 2002, Beckstein et al. 2001, Dzubiella and Hansen 2004b, Rasaiah et al. 2008) it was shown that the pore can exhibit capillary evaporation, i.e., nanobubble formation, with a strong dependence on the geometry (R, L) of the pore.This can be readily understood by the simple macroscopic arguments presented in the introduction.By evaporating from the cylinder interior, the solvent saves unfavorable interactions with the hydrophobic cylinder walls ( G ∝ −R L) while new interfaces have to be created at the cylinder mouths ( G ∝ R 2 ).In the shown example, the pore radius is ∼ 0.5 nm while the length is about 1.5 nm and a stable nanobubble is observed.Consequently, ions in the reservoir do not seem to permeate the channel as they have to pay a large solvation free energy when crossing the bubble.Under the bottom line, the channel, although sterically large enough, is blocked by the nanobubble for particle permeation.
However, it was further demonstrated that an ionic concentration imbalance between the two reservoirs can give rise to a large electric field E across the membrane that drives the polar solvent into the pore.It was observed that the electric field necessary for wetting of the pore had to be larger than a critical field E c .Once the pore was solvated, ions could be easily dragged through it by the help of the driving field, as shown in Figure 3 b).Subsequent permeation events lower the concentration gradient and the electric field E falls below E c , whereupon the pore empties again from water and closes to ion permeation.The nanobubble acts hereby as an 'electrostatic switch' and may be relevant for biological voltage gating or biotechnological applications as previously proposed (Dzubiella andHansen 2004b, 2005).Furthermore, explicit MD simulations of atomistically resolved ion channel proteins revealed that this mechanism may indeed play a key role in the electrostatic control of water and ion permeation in biological membranes (Sotomayor et al. 2007).

A. VARIATIONAL IMPLICIT SOLVENT MODEL (VISM)
Consider an assembly of arbitrary solutes surrounded by solvent in a region W. V is the region occupied by the solutes (the cavity region empty of solvent) and we identify the solute-solvent interface to be the boundary of region V, and denote it by = ∂V.We further define a volume exclusion function 1 else, characterizing space devoid of solvent.We assume that the position of each solute atom r i and the solute conformation is fixed.Thus, the solutes can be considered as an external potential to the solvent without any degrees of freedom.Coupling the VISM model the solute's degrees of freedom is work in progress (Cheng et al. 2009).In this continuum solvent model, the solvent density distribution is given by a sharp-interface distribution ρ( r ) = ρ 0 v( r ), where ρ 0 is the bulk density of the solvent.
The basic idea of VISM is to express the solvation free energy G as a functional of the cavity geometry and obtain the 'optimal' interface by minimization.It has been proposed that (Dzubiella et al. 2006a, b) where the first term is proportional to the volume of V, is the energy of creating a cavity in the solvent against the difference in bulk pressures P between the liquid and vapor phase, P = P l − P v (typically one can assume P v 0 at normal conditions and for nanoscale volumes).The second term describes the energetic cost due to the solvent rearrangement around the cavity, i.e., near the solute-solvent interface , in terms of a function γ ( r , ) with dimensions of free energy/surface area and depending on the particular topology of the solute-solvent interface.The approximation γ ( r , ) = γ 0 [1 − 2τ H ( r )] based on scaled-particle theory assumptions (Reiss 1965) accounts for the interface curvature dependence of the surface tension (Dzubiella et al. 2006a, b), where γ 0 is the constant solvent liquid-vapor surface tension for a planar interface, τ is a constant first-order curvature correction coefficient, and is the local mean curvature in which κ 1 ( r ) and κ 2 ( r ) are the local principal curvatures.The third term is the total energy of the non-electrostatic, van der Waals type, solute-solvent interaction.The potential is the sum of U i that describes the interaction of the ith solute atom (with N total atoms) centered at r i with the surrounding solvent.Each term U i includes the short-ranged repulsive exclusion and the long-ranged attractive dispersion interaction between each solute atom i at position r i and a solvent molecule at r .Classical solvation studies typically represent the interaction U i as an isotropic Lennard-Jones (LJ) potential, with an energy scale , length scale σ , and center-tocenter distance r .
The last term G ele [v] is the electrostatic energy due to charges possibly carried by solute atoms and the ions in the solvent and is modeled by a functional term on the Poisson-Boltzmann (PB) level (Dzubiella et al. 2006b).We refer to literature (Dzubiella et al. 2006a, b) and a recent work on a more mathematical derivation (Che et al. 2008) for further details and references.
The solvation free energy G is the reversible work to solvate the solute and is given by G = G[v min ] − G 0 , with some arbitrary reference energy G 0 and v min ( r ) obtained from minimization at every point of the boundary , leading to the main result where we introduced the Gaussian curvature K = κ 1 κ 2 .This partial differential equation (PDE) for the optimal exclusion function v min ( r ), or equivalently the optimal solute-solvent interface min , is expressed in terms of pressure, curvatures, short-range repulsion, and dispersion, all of which have dimensions of energy density.It can also be interpreted as a mechanical balance between the forces per surface area generated by each of the particular contributions.Note that the electrostatic contribution is not included in (5); see the literature (Che et al. 2008, Dzubiella et al. 2006b) for details.

B. THE LEVEL-SET METHOD
The level-set method is a stable numerical technique for the description and evolution of 2D surfaces in 3D space and perfectly suited for solving eq. ( 5).One of the major advantages of the level-set method is its easy handling of topological changes of surfaces during the surface evolution.For instance, the merge or break of bubbles can be captured by level-set calculations.

9
The starting point of the level-set method is to identify the surface in three-dimensional space as the zero level-set of a function φ = φ( r ) [see literature (Cheng et al. 2007) for references] = r : φ( r ) = 0 .In other words, the surface consists exactly of those points r at which the function φ vanishes.The function φ = φ( r ) is called a level-set function of the surface .The unit normal vector n at the interface , the mean curvature H , and the Gaussian curvature K can all be expressed in terms of the level-set function φ: where H e(φ) is the 3 × 3 Hessian matrix of the function φ whose entries are all the second order partial derivatives ∂ 2 i j φ of the level-set function φ, and adj (H e(φ)) is the adjoint matrix of the Hessian H e(φ).The basic idea is now to track the motion of the moving surface (t) by evolving the level-set function φ( r , t) and its zero level-set at each time t according to the so-called level-set equation, where v n = d r (t)/dt • n is the normal velocity.As in common practice, we define the level-set normal velocity to be proportional to the negative of the first variation of the Gibbs free energy v n = −μδG[v]/δv( r ) yielding with ( 5) where μ is the interface mobility.For the equilibrium solution of (5), the particular form of μ is arbitrary.If one is interested in the right dynamics and time scales, however, μ has to be further defined.Although only valid for spherical bubbles, it was shown by a simple analysis of the Rayleigh-Plesset equation that μ = R/4η, where R is the radius of mean curvature of the bubble and η the dynamic solvent viscosity (Dzubiella 2007).
For complex geometries, one could naively assume that a simple generalization μ = (4ηH ) −1 holds, but here a more thorough derivation based on the hydrodynamic Stokes equations is necessary.Interestingly, however, if one assumes dominance of capillary pressure (∝ γ 0 H ) in the driving of the nanobubble interface motion, one immediately sees that the mean curvature cancels out and the absolute interface velocity is given by v n ∼ γ 0 /η, which is of the order of 1Å/ps (100m/s).This is in accord with previous explicit-water computer simulations of nanobubble collapse (Huang et al. 2003, Liu et al. 2005, ten Wolde and Chandler 2002, Zhou et al. 2004).

C. SELECTED VISM EXAMPLES
The following selected examples illustrate that complex cavities with different interface curvature, as typical in molecular structures, can be efficiently tackled by the VISM together with the level-set method.For more quantitative details of the systems including the calculation of solvation free energies, see previous work (Cheng et al. 2007).In all of our examples, we focus on water close to normal conditions (T = 298K and P = 1bar).All solutes are modeled by fixed assemblies of identical and uncharged spheres with hard-core or LJ interactions.

Two spherical hydrophobic nanosolutes
We applied the VISM to study the simple two sphere system discussed in section 2 (Dzubiella et al. 2006a).
The results for the MF are shown in Figure 2 together with the explicit modeling results and show good agreement.In particular, VISM captures the trend that nanobubble formation is increasingly suppressed (Dzubiella et al. 2006a) and the strong attraction in the MF reduces with charging of the solutes.Thus, VISM captures the subtle competition between hydrophobic (interfacial) and electrostatic effects.

Two parallel nanometer-sized paraffin plates
Here we consider the strongly hydrophobic system of two parallel paraffin plates as investigated in explicit water MD simulations (Koishi et al. 2004) As can be seen in Figure 4, we have almost perfect agreement with the MD simulation results.Our variational implicit solvent model captures the nanobubble formation, see the interface snapshots in Figure 4 for d < 15Å, a consequence of the energetic desire of the system to minimize the total interface area.For larger distances the interface breaks, i.e., the bubble ruptures and the equilibrium interface is given by two isolated plates having almost no mutual interaction.

Helical alkane chains
In order to check if our level-set method works for more complex shaped molecules, we study model helical alkane chains.We investigate the solvation of two different configurations, one is more loosely packed with 32 atoms, and the other uses more tightly packed 22 atoms with hardly room for water in the helical core.
The results are plotted in Figure 5, which include a comparison of our level-set calculation and a typical SAS type surface constructed by taking the envelope of all the LJ spheres.Though they look quite similar, our level-set result leads to a much smoother solute-solvent interface, and a nanobubble in the tightly packed helix, a result from the minimization of the interface area based on the energy functional.In Figure 5, on the right side, we show the same two helical chains but using a color code that indicates the value of mean curvature at each point of the solute-solvent interface.The curvature ranges between positive and negative values showing also concave parts (blue).The highest positive curvature (deep red) is roughly given by the inverse of the length σ of one LJ sphere.According to the energy functional, the concave parts give large interfacial energy penalties (mainly entropic) but the latter might be compensated by large favorable dispersion attraction (mainly enthalpic), which result in total solvation energies close to small alkanes (Cheng et al. 2007).Note that the concave parts of the interface are on the inner side of the loosely packed helix in contrast to the convex outer parts.This finding suggests a qualitatively different hydration of the helix based on local curvature, i.e., geometry.

CONCLUSIONS
In conclusion, capillary effects of water on nanometer scales (e.g., nanobubble formation) have strong relevance in a broad range of confining systems as found in bio(techno)logical or nanofluidic processes.As demonstrated by explicit studies of these effects in generic model systems, a subtle competition between interfacial and energetic contributions (e.g., electrostatics, dispersion) governs the stability and occurrence of nanobubbles and, consequently, the behavior and function in any possible nano-devices.The implicit modeling of these intricate effects poses a formidable theoretical and numerical challenge for the 'nano'-physical and mathematical community.The presented VISM together with the level-set numerics shows high potential for that and, after further refinement, might be a suitable candidate for the implicit modeling and prediction of nanobubbles in hydrophobic confinement, and for molecular solvation studies in general.Palavras-chave: solvatação, hidrofobicidade, nano-bolhas, modelo de água implícito, simulações de dinâmica molecular.

Fig. 1 -
Fig.1-Sketch of a hydrophobic plate-like confinement on nanometer scales.Left: the plates trigger capillary evaporation and the region between the two plates is devoid of water (nanobubble).Right: Charging of the plates will drag the polar water between the plates and induces filling, i.e., bubble rupture.

Fig. 2 -
Fig. 2 -Left: Cylindrically averaged water density profiles around two hydrophobic (hard-core) solutes from explicit-water MD simulations (Dzubiella and Hansen 2003, 2004a).The solutes are fixed at a surface-to-surface distance of 4Å.Dark regions show low density of water while light regions show high water density.The solutes (spheres) are homogeneously and oppositely charged with z = a) 0, b) 2, and c) 5 elementary charges.Right: Mean force between the two solutes vs. their distance.The hydrophobic range and depth of attraction decreases with an increasing charge.Symbols are MD results (squares: z = 0, diamonds: z = 2, triangles: z = 5), while lines are the according VISM results.
Fig. 3 -MD simulation snapshots of a cylindrical hydrophobic pore (dashed lines) connecting to reservoirs explicitly filled with water and ions (red and green spheres) (Dzubiella and Hansen 2004b, 2005).a) The hydrophobic pore triggers capillary evaporation and is blocked for ion permeation.b) A field E larger than a critical field E c induces water permeation and subsequent tunneling of ions.The nanobubble acts as an 'electrostatic switch'.

Fig. 4 -
Fig. 4 -Comparison of level-set and MD calculations for the twoplate system.The dotted line is the pmf by MD simulations and circles indicate values of pmf by the level-set calculations using τ = 0.9Å.Three configurations of the system are also shown in their respective separations.
Fig. 5 -A comparison of the level-set (left) and SAS (middle) description of helical polymer chains with 32 atoms (above) and 22 atoms (below), respectively.Right: color represents mean curvature.The units of mean curvature are in Å −1 .