Are hydrodynamic interactions important in the kinetics of hydrophobic collapse?

Are hydrodynamic interactions important in the kinetics of hydrophobic collapse?

Jingyuan Li Department of Chemistry, Columbia University, 3000 Broadway, MC 3103, New York, NY 10027 Chinese Academy of Sciences Key Lab for Biomedical Effects of Nanomaterials and Nanosafety, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China    Joseph A. Morrone Department of Chemistry, Columbia University, 3000 Broadway, MC 3103, New York, NY 10027    B. J. Berne Department of Chemistry, Columbia University, 3000 Broadway, MC 3103, New York, NY 10027

We study the kinetics of assembly of two plates of varying hydrophobicity, including cases where drying occurs and water strongly solvates the plate surfaces. The potential of mean force and molecular-scale hydrodynamics are computed from molecular dynamics simulations in explicit solvent as a function of particle separation. In agreement with our recent work on nanospheres [J. Phys. Chem. B 116, 378 (2012)] regions of high friction are found to be engendered by large and slow solvent fluctuations. These slow fluctuations can be due to either drying or confinement. The mean first passage times for assembly are computed by means of molecular dynamics simulations in explicit solvent and by Brownian dynamics simulations along the reaction path. Brownian dynamics makes use of the potential of mean force and hydrodynamic profile that we determined. Surprisingly, we find reasonable agreement between full scale molecular dynamics and Brownian dynamics, despite the role of slow solvent relaxation in the assembly process. We found that molecular scale hydrodynamic interactions are essential in describing the kinetics of assembly.

I Introduction

Hydrophobic interactions are important in structural biology and nanoscience where they are often the major driving force in self-assembly. Berne et al. (2009); Wallqvist and Berne (1995); Chandler (2005); Lum et al. (1999); Jamadagni et al. (2011); Huang et al. (2003); Willard and Chandler (2008); Patel et al. (2010); Zangi (2011) For example, hydrophobic interactions are acknowledged to play a major role in the formation and folding of proteins. Zhou et al. (2004); Liu et al. (2005) A classic example of hydrophobic assembly is when two plates that bear little attraction to water are separated by less than some critical separation, capillary evaporation occurs in the inter-plate region, thereby driving association.

Although a great deal of work has concentrated on the structural and thermodynamic aspects of hydrophobic assembly, less is known about the kinetic mechanisms involved in the association process. Prior studies have investigated the rate of evaporation of solvent confined between two hydrophobic surfaces, Bolhuis and Chandler (2000); Leung et al. (2003); Huang et al. (2003); Luzar (2004); Sharma and Debenedetti (2012) although the full association pathway remains largely unexplored. One possible route to modeling the kinetics is within the framework of Brownian (Smoluchowski) dynamics, a coarse grained description in which the solvent degrees of freedom are not explicitly treated and where the dynamics of the heavy bodies are instead modeled in a stochastic bath. The effect of solvent is encoded in the potential of mean force (PMF) between bodies and the friction coefficient. If hydrodynamic interactions (HI) are included, then the friction coefficients depend on the positions of the heavy bodies. Within this framework, one can simulate protein diffusion and association using the appropriate equations together with the potential of mean force and the (position dependent) friction tensor as input. Ermak and McCammon (1978) For example, such an approach has been utilized to study crowded cellular environments, where the diffusion of proteins is thought to slow down considerably due in part to hydrodynamic interactions. Ando and Skolnick (2010)

Most typically, hydrodynamic interactions that are utilized in Brownian dynamics simulation are computed within a continuum description of the solvent. Brady (1988); Brady et al. (1988) This description, however breaks down for inter-particle length scales of 1-2 nm. Thomas and McGaughey (2009); Bocquet and Charlaix (2010) On such length scales, fluctuations which accompany solvent layering or capillary drying could have an enormous effect on hydrodynamic interactions. In recent work, we have shown how the nature of solvent confined between nanoscopic spheres determines the role of HI in assembly. Morrone et al. (2012) There is a growing interest in the behavior of water molecules in confined geometries. Giovambattista et al. (2007); Rasaiah et al. (2008); Baron et al. (2010) The diffusion constant and the hydrogen bond lifetime of water molecules in confined environments are distinct from those in the bulk. Li et al. (2005) Furthermore, similar signatures have been recently observed in simulations of crowded protein solutions. Harada et al. (2012) We believe such effects can potentially play an important role in the transport and association of nanoscopic bodies.

In this paper, we extend our previous study of nanospheres to the study of hydrodynamic interactions and association kinetics of plates of varying hydrophobicity. We calculate the friction on two parallel plates as a function of their separation from molecular dynamic simulations with explicit water using a technique previously employed for spherical solutes. Bocquet et al. (1997); Morrone et al. (2012) We also calculate the potential of mean force between the two plates as a function of their separation. Both the molecular-scale effects of hydrodynamic interactions and the free energy profile are found to be highly correlated with the behavior of water molecules in the inter-plate region, in agreement with our prior work. Morrone et al. (2012)

We presently consider three types of graphene-like plates: plates with a “fully” attractive solute-solvent interaction potential, plates with a “reduced” attractive solute-solvent potential, and plates with a purely repulsive solute-solvent potential. For the plates with either a reduced attractive potential or a purely repulsive potential, capillary evaporation occurs when the inter-plate separation is smaller than a critical value. In contrast, capillary evaporation does not occur for plates with a fully attractive solute-solvent potential. As the fully attractive plates approach each other, we observe solvent layering with water molecules being eventually expelled from the inter-plate region due to steric repulsion. Both dewetting in the former cases and solvent layering in the latter case greatly affect the solvent fluctuations in the inter-plate region and give rise to molecular scale effects in the hydrodynamic interactions. We compute the spatial dependence of the friction tensor and investigate the relation between the static and dynamic fluctuations of water density in the inter-plate region and this property.

Given the potential of mean force and position dependent friction coefficient we compute the rate of diffusion controlled association by means of Brownian dynamics (BD) simulations. The predictions of BD are tested against the results garnered from a set of molecular dynamics (MD) simulations of assembly in explicit solvent. This serves as a test of the validity of the Markovian approximation inherent in Brownian dynamics simulation, and the ability of coarse-grained stochastic dynamics to adequately capture the kinetic effects associated with the bath. We find that there is reasonable agreement between the two simulations when molecular scale hydrodynamic interactions are included in the BD, but marked disagreement when hydrodynamic interactions are neglected. The observation that BD with hydrodynamic interactions is in agreement with the MD result is somewhat surprising since, as we show, slow solvent fluctuations play an important role in the assembly process, and we would thus expect non-Markovian effects to be of importance. The agreement between BD and MD indicates that the non-Markovian effects may be encoded in the spatial dependence of the friction coefficient in an average way.

Ii System and Method

We modeled the hydrophobic plate as a graphene-like sheet with an area of nm as shown in Fig. 1. The plate size was chosen so as to facilitate both nanometer scale drying transitions and comparison with our prior work on nanospheres. In order to vary the hydrophobicity we studied three types of plates, all with a carbon-carbon bond length of nm, and a Lennard-Jones diameter of nm but with different interaction potential well depths: (a) full Lennard-Jones (LJ) interaction with = 0.276 kJ/mol, Werder et al. (2003) (b) reduced LJ interaction with = 0.055 kJ/mol, and (c) purely repulsive Weeks-Chandler-Anderson (WCA) truncation of the full LJ potential  Weeks et al. (1971). For convenience, these three types of plates are denoted as (full) LJ plates, reduced LJ plates and WCA plates in the present work. Solute-solvent interactions are given by the geometric mean of the respective water and solute parameters.

In order to calculate the relative friction coefficient, two parallel plates are placed perpendicular to the z-axis and a series of simulations are run with the plates fixed at various separations, and all internal degrees of freedom are frozen. The inter-plate separation ranges from 0.4 to 2 nm (for the full and reduced LJ plates) or to 2.2 nm (for the WCA plates). The solvent-induced mean force is computed from such fixed configurations, although a finer grid spacing in is employed. In order to obtain the potential of mean force (PMF), this quantity is then integrated in combination with the direct plate-plate interactions.

In the simulations of the association process, two plates are initially placed nm apart. The lower plate is fixed to its position and a biased force along the direction is applied to the upper plate. Harmonic potentials are utilized to treat the intramolecular stretches ( kJ/mol nm and nm) and bends ( kJ/mol rad and ). Snapshots from the association of reduced LJ plates and representative trajectories plotted along the reaction coordinate are shown in Fig. 1. The results from these two sets of calculations are discussed in Section IV.2 and Section IV.3, respectively. In addition, we determine the friction coefficient on a single plate, in a system containing only one plate, by pulling it at a constant speed perpendicular to its plane and computing the drag force along its normal direction.

All systems are solvated by TIP4P water. Jorgensen and Madura (1985) The size of the solvation box of the full LJ and reduced LJ system is nm, and for the WCA system the box size is increased to nm. This is done in order to accommodate the large volume fluctuation that arising from the strong dewetting transition.

In the calculation of the friction, each system is equilibrated with a 2 ns NPT simulation (P= 1 atm, T = 300 K). The Berendsen method  Berendsen et al. (1984) and the stochastic velocity rescaling method Bussi et al. (2007) are chosen for the barostat and thermostat, respectively. For each system, there are 12-18 NVE production runs of 4 ns to compute the friction coefficient. In order to facilitate energy conservation, double precision routines are utilized for all NVE runs. To study the dragging force of the single plate, we choose four pulling rates ranging from 0.5 to 3 nm/ns. There are at least 12 NPT runs for each pulling rate that are utilized to obtain an accurate estimate of dragging force. We perform at least 55 runs with various initial configurations for each system to fully characterize the association process of two plates and estimate the distribution of mean passage times. All simulations are performed using GROMACS 4.5.3 and the Particle-Mesh Ewald technique is utilized to treat long range electrostatics. Darden et al. (1993); Hess et al. (2008)

Figure 1: (Upper left panel) The top view and side view of two parallel graphene-like plates. (Upper right panel) Representative association trajectories of two plates with the “reduced” LJ interaction. (Lower panels) Configurations taken from the association of the reduced LJ plates, corresponding to an inter-plate separation of (from left to right): = 1.7 nm (initial configuration), 1.0 nm (where the association stagnates), 0.7 nm (during the driving induced collapse) and 0.35 nm (when the plates come into contact). The plates and water molecules in the inter-plate region are depicted by a space filling representation, and the snapshots are rendered with VMD. Humphrey et al. (1996)

Iii Brownian dynamics with hydrodynamic interactions

Within the framework of Brownian dynamics, hydrodynamic interactions are encoded in the frictional force. The stochastic equation for the association process contains contributions from the mean force, the frictional force, and the Gaussian random force. The potential of mean force, , includes the contributions from the direct plate-plate interaction and the solvent-induced interaction. The separation between two plates, , is treated as the reaction coordinate. Hydrodynamic interactions give rise to the spatial dependence of the relative friction coefficient, and through this the spatially dependent diffusion coefficient .

The stochastic BD equation for plate association can be therefore written as: Tough et al. (1986)


where and the diffusion coefficient is related to the “random force” through the fluctuation dissipation theorem,


This equation of motion is an accurate description of the dynamics in the high friction limit for timescales greater than the momentum relaxation time, and in which the solvent time-scale is fast compared to the time scale for the motion of the heavy bodies. The dynamics can alternatively be expressed as a Smoluchowski equation for the probability distribution of finding the particle at at time :


This equation or BD can be used to determine mean first passage times for diffusion controlled reactions once and are known.

In order to apply BD we must determine the PMF and spatial dependence of the friction coefficient along the reaction path. As already pointed out, the continuum treatment of the hydrodynamic interaction breaks down at small separations. Thomas and McGaughey (2009); Bocquet and Charlaix (2010); Morrone et al. (2012) In this work, we calculate the friction coefficient on the molecular length scale from explicit molecular dynamics simulations. The friction coefficient can then be determined from the two-body friction tensor . In the calculation, the plates are fixed as required in the Brownian limit. The pair friction tensor can be expressed as a time integral of the correlation functions of the fluctuating force ,


where by symmetry and Here we focus on the direction parallel to the inter-plate separation and therefore study the force components along the plate normal vector. Eqn. 4 is only valid in the limit where the number of solvent molecules, , approaches infinity. Bocquet et al. (1997) In the present case of finite systems, it is still possible to relate the force autocorrelation function to the friction tensor, following the procedure of Bocquet et al. Bocquet et al. (1997) The friction along the relative coordinate is given by and depicted in Section IV.2.

Iv Results

iv.1 Friction and diffusion coefficients of a single plate

The friction coefficient perpendicular to the face of a single plate immersed in solvent may be extracted from the force autocorrelation function by means of techniques outlined in Ref. Bocquet et al., 1994. The resultant calculation yields values of , , kg/s, when the solvent-solute interaction is the full LJ, Reduced LJ, and the WCA potential, respectively. In agreement with prior work, Morrone et al. (2012) the friction experienced by the body decreases with increasing surface hydrophobicity. Each calculation is performed in a periodic box of dimensions given in Section II.

The friction coefficient may also be extracted from a set of non-equilibrium simulations where the plate is dragged through the solvent at a given velocity. As depicted in Fig. 2, the dragging force at four velocities is determined from MD simulations in explicit solvent, and grows with increasing velocity. The resultant plot of dragging force versus pulling rate exhibits a linear relationship with the slope equal to the friction coefficient. The extracted values of the friction are given in Table 1, alongside the results computed from the force autocorrelation function. The two estimates of the friction are in close agreement. The consistency between the results garnered from two different techniques serves to validate the calculation.

Next, we compare these results to the predictions of continuum hydrodynamics. The friction coefficient for a rigid body described by a set of spherical sites can be evaluated by means of the Rotne-Prager tensor. Rotne and Prager (50) The translational friction coefficient may be extracted from the full mobility tensor, where is the number of sites. Carrasco and de la Torre (1999) In the present calculation, periodic boundary conditions must be taken into account and are known to have a significant impact on hydrodynamic properties. Hasimoto (1959); Beenakker (1986); Yeh and Hummer (2004) Periodic effects may be treated by means of replacing the standard Rotne-Prager tensor with a form that accounts for periodicity via an Ewald summation. Beenakker (1986) This approach is discussed further in the Appendix.

As discussed above, different solvation boxes are utilized in the WCA and Lennard-Jones systems, and the continuum calculation must be undertaken for both sizes. The values of the friction coefficient with stick boundary conditions for the two box sizes are given in Table 1. The full LJ value is somewhat larger than the continuum result whereas, in line with expectations, the increasingly hydrophobic plates begin to show larger deviations from the stick boundary condition.

Figure 2: The dragging force of single plate moving at constant velocity. The results for the plate with full Lennard-Jones(LJ) interaction (LJ plate), the plate with reduced LJ interaction (reduced LJ plate), purely repulsive Weeks-Chandler-Anderson truncation of the full LJ potential (WCA plate) are depicted in black, red, and green, respectively. The dragging forces increase linearly as the pulling rate increases, and the fitted slopes are equal to the friction coefficient.
interaction box (nm) ( kg/s) ( kg/s) ( kg/s)
WCA 1.03 1.27
REDUCED 1.33 1.36 1.53
FULL 1.63 1.69 1.53
Table 1: The friction on a single plate in the -direction, as estimated from the force autocorrelation function, pulling simulations, and continuum hydrodynamics.

iv.2 Hydrodynamic and thermodynamic profiles

Friction coefficients, , and the potential of mean force, , are depicted as a function of inter-plate separation, , in Figs. 3, 4, and 5 (for the full LJ, reduced LJ, and WCA plate, respectively). The frictional profiles exhibit non-monotonic behavior as the two plates approach each other. The spatially dependent features of the molecular scale hydrodynamic interactions display similar trends to those found in the free energy profile. This finding is in agreement with our previous work, Morrone et al. (2012) and the recent work of Mittal and Hummer. Mittal and Hummer (2012) In the case of the full LJ plate, the friction peaks and the free energy increases as a layer of water is “squeezed out” Zangi (2011) at nm. The friction coefficient subsequently decreases at the minimum of the PMF. The final solvent layer is expelled as the separation decreases to 0.62 nm, and both the free energy and the friction coefficient increase. It is important to note that the friction profile peaks at 0.88 and 0.62 nm which are at the same positions as the PMF barriers. The peak heights at these two separations are and kg/s, respectively. Standard errors of the mean value are given in parenthesis.

Figure 3: Spatial dependence of relative friction coefficient (A) and the potential of mean force (B) for the association of two LJ plates. The friction peaks at = 0.62 nm and two free energy wells at = 0.68 and 0.97 nm are depicted in the insets of (A) and (B), respectively. (C) The relative density of water in the inter-plate region. (D) The ratio of the variance to the average of the number of water molecules (black) and the solvent relaxation time in the inter-plate region (red). For both curves in (D), the peaks at 0.9 nm are depicted in the inset. The values corresponding to the friction peaks are shown as filled symbols in panels (C) and (D).

For the reduced LJ plates, there is a low barrier at nm in the PMF, whereas the WCA plates exhibit barrierless assembly along the chosen reaction coordinate. In these cases, the corresponding friction profiles also display non-monotonic behavior. The friction profile peaks at = 0.92 nm and = 1.35 nm for the reduced LJ and WCA plates, respectively. As discussed below, these distances are in the region of the critical separation for dewetting. The corresponding peak heights of the friction coefficients are (reduced LJ) and kg/s (WCA). The somewhat large error bar results from the sizable force fluctuations that are present in the vicinity of .

The friction profiles converge to (LJ), (reduced LJ) and (WCA) kg/s at large separations. The standard error in the mean for these values is . These numbers are somewhat larger than expected upon comparison with the results given in Table 1. This difference may largely result from the finite box size and the impact of neighboring periodic images. However, the effect of periodic boundary conditions should not greatly impact the observed spatial dependence of the short-range hydrodynamics interactions.

In order to better understand the nature of molecular-scale hydrodynamic interactions, we analyze the density and the static and dynamic fluctuations in number of water molecules in the inter-plate region. The static fluctuations are measured by the ratio of the number variance to the average . An estimate of the solvent relaxation time can be computed from the integral of the normalized autocorrelation function of these fluctuations, . The resultant plots of these quantities as a function of inter-plate separation are depicted in Figs. 3-5. Additionally, the values of the solvent density and fluctuations at the separation that corresponds to the maximum friction are marked on the curves. For the LJ plates, the solvent density decreases when the plate separation decreases from to nm. The ratio of the variance to the average of water density peaks at nm, and the relaxation time peaks at nm. Both peaks occur at separations close to where the friction coefficient peaks. Moreover, when the separation decreases from to nm the solvent density sharply decreases. The static fluctuations and relaxation time grow in the process of water expulsion. Both the static variance of water density and the relaxation time peak at the same separation where the friction coefficient peaks ( = 0.62 nm).

Figure 4: Spatial dependence of relative friction coefficient (A) and the potential of mean force (B) for the association of two reduced LJ plates. (C) The relative density of water in the inter-plate region. (D) The ratio of the variance to the average of the number of water molecules (black) and the solvent relaxation time in the inter-plate region (red). The value corresponding to the friction peak is shown as filled symbols in panels (C) and (D).

In both the reduced LJ and WCA systems, the solvent density in the inter-plate region begins to decrease at the critical separation of the dewetting transition. For the reduced LJ plate, the solvent density dramatically decreases when the separation decreases from 0.95 to 0.92 nm. Both the variance of water density and the solvent relaxation time peak at nm. The density of water between the WCA plates decreases as the separation decreases from to nm. The ratio of the variance to the average of water density and the solvent relaxation time peak in the same region. At the critical separation for the dewetting transition, the inter-plate region fluctuates between wet and dry. The value of the critical distance between surfaces is on the order of one nanometer for the present-sized plates and decreases with decreasing hydrophobicity, in agreement with macroscopic thermodynamic analysis. Lum and Luzar (1997); Huang et al. (2003); Sharma and Debenedetti (2012) This characterization of the density and the fluctuation of solvent is consistent with the results of previous studies. Huang et al. (2003); Morrone et al. (2012)

In the three systems presently studied, it can be clearly seen that the friction coefficient increases where the solvent fluctuations become large and slow. Taken together, both static and dynamic solvent behavior engender the large frictions at small inter-plate separations. In agreement with our prior work, Morrone et al. (2012) molecular-scale hydrodynamic interactions largely result from such fluctuations when, in the case of hydrophobic bodies, the drying transition occurs or, for less hydrophobic species, when water molecules are expelled from the inter-plate region due to steric repulsion.

Figure 5: Spatial dependence of relative friction coefficient (A) and the potential of mean force (B) for the association of two WCA plates. (C) The relative density of water in the inter-plate region. (D) The ratio of the variance to the average of the number of water molecules (black) and the solvent relaxation time in the inter-plate region (red). The value corresponding to the friction peak is shown as filled symbols in panels (C) and (D).

iv.3 Comparison of Brownian dynamics with molecular dynamics

In order to further elucidate the impact of molecular-scale hydrodynamics on the kinetics of self assembly, we perform direct molecular dynamics simulations of this process. The two plates are initially placed perpendicular to the z-direction at a separation of 1.8 nm. A constant loading force is added to the upper plate and the lower plate is fixed in its initial position. The loading force utilized is pN (full LJ) or pN (reduced LJ and WCA). There are at least fifty-five association simulations performed for each system.

The assembly of both the WCA plates and reduced LJ plates slows down around the critical separation for dewetting. This process is illustrated for the reduced LJ case by the sample trajectories depicted in Fig. 1Humphrey et al. (1996) One can see there is a large dwell time at nm. In the case of the full LJ system, the upper plate initially rapidly approaches the lower plate and the inter-plate distance decreases to nm, corresponding to two inter-plate water layers. The association process then stagnates. After some period a water layer is expelled and the plate separation decreases to nm where a single water layer separates the plates. There is another dwell time before the last water layer is expelled and the two plates finally come into contact. During the association process the upper plate can rock. This is particularly pronounced around nm where the amplitude of the rocking can be as much as nm. In order to provide the best comparison with the (one-dimensional) BD scheme outlined below, the plates must be kept as parallel as possible. This is facilitated by the addition of a harmonic restraining potential on the plate’s internal degrees of freedom. It is important to note that outside the present context, “rocking” could be a viable degree of freedom, and we observe that assembly occurs significantly more rapidly when such effects are included. In the case of less attractive or purely repulsive solvent-solute interactions, the effect is much less prominent, as the displacement generated by the rocking mode is small with respect to the critical separation for dewetting that drives assembly ( nm).

As we will show, the mean first passage times observed in the MD association trajectories cannot be predicted by solely considering the free energy profile with constant friction. In order to evaluate the role of hydrodynamic interactions in assembly, we utilize a one-dimensional Brownian dynamics (BD) scheme as described in Section III, where the system evolves according to Eqn. 1 and can be integrated as described by Ermak and McCammon. Ermak and McCammon (1978) The dynamical degree of freedom is taken to be the inter-plate distance, , along which the friction and potential of mean force have been computed (see Figs. 3-5).

type (pN) (ps) (ps) (ps) (nm) (nm)
WCA 20 515 516 134 1.60 0.40
reduced LJ 20 1469 2390 200 1.20 0.40
LJ (1st barrier) 440 573 1460 85 1.20 0.80
LJ (2nd barrier) 440 32600 23078 243 0.75 0.40
Table 2: Mean first passage time, , of the plate association process as described in Section IV.3

This scheme may be utilized to generate predictions for the mean first passage time (mfpt) and distribution of first passage times (dfpt) from BD simulations with the loading force. Such a comparison serves to evaluate the degree to which the Brownian framework can yield an adequate description of the kinetics of assembly. The constant force applied to the molecular dynamics simulations is accounted in the BD by a means of the modified potential of mean force ,


where is the PMF determined from MD in the absence of a loading force. Because the loading force in Eqn. 5 is independent of the solution degrees of freedom, Eqn. 5 is rigorous. Moreover, because the loading potential is linear in we assume it will not alter the spatial dependence of the friction. However, for very large loading forces we expect that non-equilibrium effects will be significant and the BD picture will fail. For each system, we have generated 10,000 BD trajectories given initial state , and an absorbing boundary at . For each system, the values of these parameters and of the mean first passage time are given in Table 2. The corresponding first passage time distributions are shown in Fig. 6. The mean first passage time may also be computed directly from the Smoluchowski equation (Eqn. 3) by means of the following expression: Zwanzig (2001)


where and is the positon of the reflecting boundary. Due to the the driving force, the reflecting boundary can be taken to large values in our calculation. This equation and an ensemble of BD trajectories will yield the same result within numerical error.

The mean first passage time obtained from MD simulation for the association of the reduced LJ plates, is about 60% smaller than the result from BD simulation with hydrodynamic interaction (BD with HI), but 7 times larger than the value estimated from BD when the spatial dependence of the friction is not included and the single plate friction (see Table 1) is utilized instead (BD without HI). The fpt distribution obtained from MD is also close to, albeit narrower than, that obtained from BD with HI. In contrast, the distribution garnered from BD without HI is much more strongly peaked. In the case of the WCA plates, the mfpt and the fpt distribution from both MD and BD with HI are very close to each other, but both are approximately 4 times larger than the BD without HI result.

The results of MD simulation are in reasonable agreement with those obtained from BD with HI for both the WCA and the reduced LJ plate systems. The association process slows down in the region around the critical separation for dewetting transition where the friction peaks. As discussed above, the behavior of friction profile at the critical separation largely results from the solvent fluctuation due to the dewetting transition. The molecular-scale effect of hydrodynamic interaction evidently contributes to the slowing down of the association process near the critical separation, and if only barriers present in the PMF are considered (as in BD without HI), association occurs far too rapidly. These results indicate that a kinetic barrier along the reaction coordinate is present at the drying transition.

Figure 6: The distribution of first passage time (dfpt) for the association process of LJ plates (panel (A), first barrier; panel (B), second barrier), reduced LJ plates (C) and WCA plates (D), obtained from the MD simulation (black), BD simulation with and without the consideration of hydrodynamic interactions (HI) (red and blue, respectively). The dfpts obtained from MD are in reasonable agreement with the result of BD with HI, and the probability of having a fpt smaller than 500 ps is much larger in the results garnered BD without HI.

The values of the mfpt of the first (second) barrier of the full LJ plate system as obtained from BD with HI are 1460 (23078) ps, much larger than the results of 85 (243) ps for BD without HI. Hydrodynamic interactions strikingly slow down the first passage time over the second barrier by about two orders of magnitude. The mfpt obtained from the MD simulation is 573 and 32600 ps, for the first and second barrier, respectively. For the first barrier, the mfpt obtained from MD simulations is smaller than the result from BD with HI, while still being much larger than the result from BD without HI. Meanwhile, for the second barrier, the mean first passage time from MD simulations is larger than the result from BD simulations with HI. The distribution of the first passage times corresponding to the second barrier obtained from MD simulations is similar to the results from BD simulations with HI. For passage over the first barrier, the distributions exhibit greater deviation.

In general, the average value and the distribution of mean first passage times of the three types of plates calculated from MD simulations is consistent with the results of BD simulations with HI. In prior work, Morrone et al. (2012) we found that hydrodynamic interactions contributed approximately 40% to the assembly of two fullerenes. In the present set of calculations, hydrodynamic interactions contribute a much larger share. Comparison of spheres and plates indicate that the contribution of the hydrodynamic interaction is enhanced as the shape is flattened. More water molecules are confined by plates than spheres of the same surface area so that, at small separations, the degree of confinement and the length scale of dewetting is increased.

V Conclusion

In this work we study the impact of hydrodynamic effects on the kinetics of assembly of two plates of varying hydrophobicity. To this end, the potential of mean force and spatially dependent friction coefficient are determined along the inter-plate separation. The results show that there is a correspondence between peaks in the PMF and peaks in the frictional profile. High values of the friction are related to large and slow solvent fluctuations in the inter-plate region. Both solvent confinement and drying phenomena can play a critical role in the kinetics of assembly. In this way, molecular scale effects shape the hydrodynamic interactions, and give rise to their deviation from continuum theory, which predicts that the friction coefficient diverges as two bodies come into contact. Wolynes and Deutch (1976); Jeffrey and Onishi (1984)

The kinetics of assembly studied by means of molecular dynamics simulations in explicit solvent with a constant loading force applied along the reaction coordinate is compared with the predictions of Brownian dynamics with the same loading force and with the potential of mean force and hydrodynamic profile extracted from the data presented in Section IV.2. Brownian dynamics is a widely used technique owing to its computational efficiency as solvent degrees of freedom are not treated explicitly. There is reasonable agreement between the mean first passage times that are obtained from the two schemes. Indeed, we find that hydrodynamic interactions are essential to produce a reasonable description of the process and neglect of spatial dependence of the friction has a large impact on the kinetics. The HI give rise kinetic bottlenecks along the association pathway, Morrone et al. (2012) over and above the barriers in the potential of mean force. Interestingly, other degrees of freedom such as plate “rocking” can in some cases significantly increase the rate of assembly, probably by facilitating waters entry and exit from the inter-plate gap. In order to describe the effect of solvent relaxation on the rocking mode within the Brownian framework, it would require a determination of the friction coefficient experienced on this degree of freedom. The exploration of all possible pathways to plate assembly is beyond the scope of this work, and presently we concentrate on the approach of two parallel plates.

It has been suggested Huang et al. (2003); Willard and Chandler (2008) that including solvent degrees of freedom in the reaction coordinate is necessary in order to properly characterize the association process. Presently, a reaction coordinate is utilized that is only a function of the plate degrees of freedom. We have shown that such a choice is reasonable provided that (molecular-scale) hydrodynamic interactions are considered along the pathway. The inclusion of hydrodynamic effects captures the effect of solvent in an average sense, and the Brownian framework would be an exact treatment in the limit where the solvent time scales are much faster than those associated with the heavy bodies. Leading from the last point, the discrepancies between the MD and BD results can be largely attributed to either the impact of the pulling force, that is the friction experienced along the force-biased surface significantly differs from that extracted from the equilibrium calculations, or to non-Markovian effects. In the case of the passage over the first barrier of the LJ system, at least some of the deviation is likely due to the large pulling force, although non-Markovian effects should have an impact both in this case and in the other system studied. It has, for instance, been shown that dynamical caging effects can be exhibited along degrees of freedom that are associated with slowly varying memory functions. Tuckerman (2010) Such phenomena can also serve to explain the long dwell times exhibited in the association trajectories.

Unfortunately, given the broad distributions involved and limited number of MD trajectories that can be harvested, it is difficult to fully gauge the impact of non-Markovian effects based on the present work. Certainly, it is observed that, in the case of strongly hydrophobic plates, assembly is almost always preceded by the “first” drying transition, that is the system does not sample successive wet and dry states near the critical separation. One should expect such a process to be intrinsically non-Markovian, although the BD with HI approach that is presently employed appears to at least capture some aspect of this “waiting period” for drying in an average way. A more precise understanding of this phenomena will be the subject of future work.

This research was supported by a grant to B.J.B. from the National Science Foundation via Grant No. NSF-CHE- 0910943. J.L. was supported by the CAS Hundred Elite Program and MOST 973 program No.2012CB932504.

Appendix A Continuum hydrodynamic treatment of the friction on a rigid body

Consider a rigid body made up of spherical sites where the origin is taken to be the center of mass. A mobility tensor may be defined where:


and where and are dimensional vectors representing the velocities and external forces on the sites that comprise the body. In the case of stick boundary conditions, the elements of are approximated by the Rotne-Prager tensor Rotne and Prager (50) where the elements, , are given as:


where is the site radius, and is the shear viscosity. For TIP4P water the value is and has been taken from the literature. Gonzalez and Abascale (2010) The matrix is given by the following expression:


where and is the vector direct product of the unit vector of displacement. If the plate is treated as a body that is centro- and axi- symmetric, then the translational friction tensor, , is given by the following expression: Carrasco and de la Torre (1999)


where is a diagonal matrix. Presently, we are interested in the friction coefficient associated with motion perpendicular to the face of the plates, and this is what is reported in Section IV.1. For periodic systems, the elements of are given by


where is the periodic image vector, and the same as in Eqn. 8 except it is evaluated over periodic images. This expression may be evaluated by means of an Ewald summation, as was shown by Beenakker. Beenakker (1986) In the case of overlapping spheres, reciprocal space contributions are excluded from the Ewald sum, as the long-range portion of the hydrodynamic interaction only includes terms for which .


  • Berne et al. (2009) Berne, B. J.; Weeks, J. D.; Zhou, R. Annu. Rev. Phys. Chem. 2009, 60, 85–103.
  • Wallqvist and Berne (1995) Wallqvist, A.; Berne, B. J. J. Phys. Chem. 1995, 99, 2893–2899.
  • Chandler (2005) Chandler, D. Nature 2005, 437, 640–647.
  • Lum et al. (1999) Lum, K.; Chandler, D.; Weeks, J. J. Phys. Chem. B 1999, 103, 4570–4577.
  • Jamadagni et al. (2011) Jamadagni, S. N.; Godawat, R.; Garde, S. Annu. Rev. Chem. Biomol. Eng. 2011, 2, 147–171.
  • Huang et al. (2003) Huang, X.; Margulis, C. J.; Berne, B. J. Proc. Natl. Acad. Sci. USA 2003, 100, 11953–11958.
  • Willard and Chandler (2008) Willard, A. P.; Chandler, D. J. Phys. Chem. B 2008, 112, 6187–6192.
  • Patel et al. (2010) Patel, A. J.; Varilly, P.; Chandler, D. J. Phys. Chem. B 2010, 114, 1632–1637.
  • Zangi (2011) Zangi, R. J. Phys. Chem. B 2011, 115, 2303–2311.
  • Zhou et al. (2004) Zhou, R.; Huang, X.; Margulis, C. J.; Berne, B. J. Science 2004, 305, 1605–1609.
  • Liu et al. (2005) Liu, P.; Huang, X.; Zhou, R.; Berne, B. J. Nature 2005, 437, 159–162.
  • Bolhuis and Chandler (2000) Bolhuis, P. G.; Chandler, D. J. Chem. Phys. 2000, 113, 8154–8160.
  • Leung et al. (2003) Leung, K.; Luzar, A.; Bratko, D. Phys. Rev. Lett. 2003, 90, 065502.
  • Luzar (2004) Luzar, A. J. Phys. Chem. B 2004, 108, 19859–19866.
  • Sharma and Debenedetti (2012) Sharma, S.; Debenedetti, P. G. Proc. Natl. Acad. Sci. USA 2012, 109, 4365–4370.
  • Ermak and McCammon (1978) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352–1360.
  • Ando and Skolnick (2010) Ando, T.; Skolnick, J. Proc. Natl. Acad. Sci. 2010, 107, 18457–18462.
  • Brady (1988) Brady, J. Annu. Rev. Fluid Mech 1988, 20, 111–157.
  • Brady et al. (1988) Brady, J. F.; Philips, R. J.; Lester, J. C.; Bossis, G. J. Fluid Mech. 1988, 195, 257–280.
  • Thomas and McGaughey (2009) Thomas, J. A.; McGaughey, A. J. H. Phys. Rev. Lett. 2009, 102, 184502.
  • Bocquet and Charlaix (2010) Bocquet, L.; Charlaix, E. Chem. Soc. Rev. 2010, 39, 1073–1095.
  • Morrone et al. (2012) Morrone, J. A.; Li, J.; Berne, B. J. J. Phys. Chem. B 2012, 116, 378–389.
  • Giovambattista et al. (2007) Giovambattista, N.; Debenedetti, P. G.; Rossky, P. J. J. Phys. Chem. C 2007, 111, 1323–1332.
  • Rasaiah et al. (2008) Rasaiah, J. C.; Garde, S.; Hummer, G. Annu. Rev. Phys. Chem. 2008, 59, 713–740.
  • Baron et al. (2010) Baron, R.; Setny, P.; McCammon, J. A. J. Am. Chem. Soc. 2010, 132, 12091–12097.
  • Li et al. (2005) Li, J.; Liu, T.; Li, X.; Ye, L.; Chen, H.; Fang, H.; Wu, Z.; Zhou, R. J. Phys. Chem. B 2005, 109, 13639–13648.
  • Harada et al. (2012) Harada, R.; Sugita, Y.; Feig, M. J. Am. Chem. Soc. 2012, 134, 4842–4849.
  • Bocquet et al. (1997) Bocquet, L.; Hansen, J.; Piasecki, J. J. Stat. Phys. 1997, 89, 321–346.
  • Werder et al. (2003) Werder, T.; Walther, J.; Jaffe, R.; Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2003, 107, 1345–1352.
  • Weeks et al. (1971) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237–5247.
  • Jorgensen and Madura (1985) Jorgensen, W. L.; Madura, J. D. J. Chem. Phys. 1985, 56, 1381–1392.
  • Berendsen et al. (1984) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690.
  • Bussi et al. (2007) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126, 014101.
  • Darden et al. (1993) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092.
  • Hess et al. (2008) Hess, B.; Kutzner, C.; van der Spoel, D. J. Chem. Theor. Comput. 2008, 4, 435–447.
  • Tough et al. (1986) Tough, R. J. A.; Pusey, P. N.; Lekkerkerker, H. N. W.; Broeck, C. V. D. Mol. Phys. 1986, 59, 595–619.
  • Bocquet et al. (1994) Bocquet, L.; Hansen, J.; Piasecki, J. J. Stat. Phys. 1994, 76, 527–548.
  • Rotne and Prager (50) Rotne, J.; Prager, S. J. Chem. Phys. 50, 4831–4837.
  • Carrasco and de la Torre (1999) Carrasco, B.; de la Torre, J. G. Biophys. J. 1999, 75, 3044–3057.
  • Hasimoto (1959) Hasimoto, H. J. Fluid Mech. 1959, 5, 317–328.
  • Beenakker (1986) Beenakker, C. J. Chem. Phys. 1986, 85, 1581–1582.
  • Yeh and Hummer (2004) Yeh, I. C.; Hummer, G. J. Phys. Chem. B 2004, 108, 15873–15879.
  • Mittal and Hummer (2012) Mittal, J.; Hummer, G. J. Chem. Phys. 2012, 137, 034110.
  • Lum and Luzar (1997) Lum, K.; Luzar, A. Phys. Rev. E 1997, 56, R6283–R6286.
  • Humphrey et al. (1996) Humphrey, W.; Dalke, A.; Schulten, K. J. Molec. Graphics 1996, 14, 33–38.
  • Zwanzig (2001) Zwanzig, R. Nonequilibrium Statistical Mechanics; Oxford University Press: New York, 2001.
  • Wolynes and Deutch (1976) Wolynes, P. G.; Deutch, J. M. J. Chem. Phys. 1976, 65, 450–454.
  • Jeffrey and Onishi (1984) Jeffrey, D. J.; Onishi, Y. J. Fluid Mech. 1984, 139, 261–290.
  • Tuckerman (2010) Tuckerman, M. E. Statisical Mechanics: Theory and Molecular Simulation; Oxford University Press: New York, 2010.
  • Gonzalez and Abascale (2010) Gonzalez, M. A.; Abascale, J. F. J. Chem. Phys. 2010, 132, 096101.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description