External Potential Modifies Friction of Molecular Solutes in Water

Stokes’s law for the friction of a sphere in water has been argued to work down to molecular scales, provided the effective hydrodynamic radius includes the hydration layer. In interpretations of experiments and in theoretical models, it is tacitly assumed that the solvent friction experienced by a solute does not depend on whether an external confinement potential acts on the solute. Using a novel method to extract the friction memory function from molecular dynamics simulations, we show that the solvent friction of a strongly harmonically confined methane molecule in water increases by 60% compared to its free-solution value, which is caused by an amplification of the slowest component of the memory function. The friction enhancement occurs for potential strengths typical of physical and chemical bonds and is accompanied by a significant slowing-down of the hydration water dynamics. Thus, the solvent friction acting on molecular solutes is not determined by solvent properties and solute-solvent interactions alone but results from the coupling between solute and solvent dynamics and thereby can be tuned by an external potential acting on the solute. This also explains why simulations of positionally constrained solutes do not reproduce free-solution diffusivities. Dynamic scaling arguments suggest similar effects also for macromolecular solutes provided the solution viscosity is sufficiently enhanced.


I. INTRODUCTION
Friction sets the fundamental time scale for all processes that occur in a solvent, ranging from molecular diffusion [1], macromolecular conformational changes [2,3], and chemical reactions [4], to protein folding [5,6]. For a sphere with radius R in a continuous solvent with viscosity η, the no-slip friction coefficient γ follows Stokes's law γ ¼ 6πηR, so that the friction force F f is for small enough velocity v given by F f ¼ γv. Modifications of Stokes's law due to water discreteness and water-solute interactions have been amply discussed [1,7]. In this paper we discuss a different modification of Stokes's law, demonstrated by molecular dynamics (MD) simulations of a single methane in water that is subject to a harmonic confinement potential of strength K. In fact, we find γ to depend sensitively on K in the range 10 2 < K < 10 4 kJ=ðmol nm 2 Þ so that γ for larger K is increased by about 60% compared to the value of γ for small (or vanishing) K. This dramatic slowing-down of the methane diffusivity in confinement is mirrored by an increase of the escape time from the first to the second hydration shell from τ esc1 ≈ 8 ps for small K to τ esc1 ≈ 18 ps for large K. The intimate coupling of solute and hydration shell dynamics had been suggested based on NMR studies that showed a solute to increase the viscosity in its hydration layer, which in turn slows down the solute diffusion, an effect that has been called secondary dynamic solvent effect [8]. Our results demonstrate a direct consequence of this coupling between solute and hydration shell dynamics: By a detailed analysis of the solute friction memory function, which we extract from our simulation trajectories using a novel method, we show that solute diffusivity and hydration shell kinetics are coupled and both influenced by the inherent time scale of solute motion, in our simulation model set by the external potential strength K.
There are various consequences of our findings: In simulation studies, it is common practice to constrain the position of a solute in order to determine spatially dependent solvation properties in inhomogeneous systems, for example, in hydrated lipid bilayer systems [9]. While this is unproblematic for static properties such as free-energy profiles, our results show that this procedure potentially perturbs kinetic properties such as the diffusivity profile. In fact, there are various alternative techniques that allow us to extract free-energy and diffusivity profiles from simulation trajectories of unconstrained solutes [10,11] and which therefore do not suffer from this potential complication.
The direct experimental test of our predictions is principally possible with optical traps, which are used to confine micron-sized plastic or metal beads in water [12]. Recent technological advances allow for the laser trapping of gold nanoparticles with a radius as small as 9 nm [13], which, however, is still much larger than the radius of a methane molecule employed in our simulations. Anti-Brownian electrokinetic trapping techniques have been demonstrated to efficiently trap single molecules with subnanometer hydrodynamic radii [14], but the achievable confinement potential strengths are rather weak. The potential-induced friction enhancement we observe for methane in water is expected to disappear for large solutes for which the inertial time scale τ m ≡ m=γ (where m denotes the solute mass) is larger than the longest solvent relaxation time. Based on time scaling arguments, we argue in Sec. V that the confinement-potential-induced modification of the solute friction coefficient should be observable experimentally also for larger solutes, if the solvent relaxation time is suitably increased by using high-viscosity solvents.
A less direct but nevertheless relevant experimental consequence concerns the coupling between the dynamics of molecular probes and the surrounding hydration shell dynamics, which is studied by various experimental techniques. It is known that the hydration shell dynamics slows down considerably upon transferring a molecular probe from bulk water to the surface of a macromolecule, in line with our results. While in the experimental system the anchoring of a molecular probe not only confines the molecular probe but also changes its environment due to the presence of a linker group and an anchoring scaffold, it is clear that the confinement-induced mechanism we demonstrate in this paper contributes also in these more complex experimental scenarios, as we explain in Sec. V.

II. FRICTION CONSTANT IN HARMONIC POTENTIAL
For an unconfined solute the diffusion constant and the friction coefficient follow from the long-time limit of the mean-square displacement. Obviously, for a solute that is confined by an external potential this standard protocol is not applicable. Instead, we have to use a more indirect route and extract the friction coefficient γ from the symmetric memory function ΓðtÞ ¼ Γð−tÞ, defined by the generalized Langevin equation [15,16], where UðxÞ ¼ 1 2 Kx 2 is a harmonic potential and the random force F R ðtÞ obeys hF R ð0ÞF R ðtÞi ¼ k B TΓðtÞ according to the fluctuation-dissipation theorem. We include the inertial term proportional to the solute mass m in order to correctly account for the short-time behavior where the solute dynamics is ballistic instead of diffusive. In some publications the upper integration boundary of the memory term in Eq. (1) is t instead of ∞. The two formulations are for positive time t > 0 related by a shift of the random force by ΔF R ðtÞ ¼ R ∞ 0 dt 0 Γðt þ t 0 Þ_ xð−t 0 Þ. Our formulation simplifies the Fourier analysis and does not require us to specify initial conditions; for a detailed discussion we refer to Ref. [17]. Our novel method to extract the memory function ΓðtÞ from simulation trajectories, which can also be used for experimental trajectories, is based on the solvent force experienced by the particle, The autocorrelation function, is after Fourier transformationC sol where we introduce the single-sided memory function Γ þ ðtÞ ≡ ΓðtÞ for t ≥ 0 and Γ þ ðtÞ ≡ 0 for t < 0 [see the Supplemental Material (SM) [18] for the full derivation].
Since the friction coefficient γ is given by γ ≡ R ∞ 0 dtΓðtÞ, we see from Eq. (4) that for K ≠ 0, Note that in the limit K ¼ 0 the integral over the force autocorrelation function C sol FF ðtÞ vanishes and therefore cannot be used to extract γ, as amply discussed in the literature [19,20].

III. SIMULATION SETUP
We perform MD simulations of a single water-solvated methane modeled as a Lennard-Jones (LJ) particle in a harmonic potential of strength K; see Fig. 1(a) for a schematic simulation setup. We use the GROMACS 5.1 [21,22] simulation package. The LJ parameters corresponding to methane are taken from the GROMOS 53A6 [23] force field, for water we use extended simple point charge (SPC/E) [24] parameters. We perform NVT simulations at fixed water number N, fixed volume V and fixed temperature T with a timestep of 2 fs at a temperature of T ¼ 300 K, controlled by the velocity rescaling [25] thermostat coupled with a time constant of 0.5 ps to water only, for three different cubic box sizes L ¼ 3 nm (894 H 2 O), L ¼ 4.5 nm, (3008 H 2 O) and L ¼ 6 nm (7159 H 2 O). In the SM [18] we compare these results to simulations in the NVE ensemble (fixed water number N, fixed volume V and constant energy E) and thereby demonstrate that the ensemble and the thermostat have no significant influence. We use simulation lengths of roughly 500 ns per parameter combination. Before the production runs, all systems are equilibrated for 5 ns in the NPT ensemble (fixed water number N, fixed pressure P and fixed temperature T) with P ¼ 1 bar using a Berendsen [26] barostat to determine the box size. The solvent force F sol ðtÞ acting on the methane is calculated from the difference of the total force, which is due to interactions with all water molecules, and the force generated by the harmonic potential; in the case K ¼ ∞, the methane is frozen at the potential minimum and the solvent force equals the total force. Position, velocity, and total force acting on the methane are saved at every time step. Fast Fourier transforms are used for the calculation of autocorrelation functions and integrals over the solvent force autocorrelation function are computed by trapezoidal integration.

IV. RESULTS
Typical trajectories of the 1D methane position xðtÞ and the conjugated total force FðtÞ ¼ mẍðtÞ for two different values of K are shown in Fig. 1(c). As expected, the methane oscillates in the confining potential, with frequency and amplitude depending on the force constant K. In Fig. 1(b), radial distribution functions (RDF) of the methane-water oxygen separation are depicted for different K, including the frozen limit K ¼ ∞. All RDFs superimpose perfectly, reflecting that the equilibrium properties of the methane hydration shell do not depend on K.
We show simulated solvent force autocorrelation functions C sol FF ðtÞ for different K as colored lines in Fig. 2(a); mild oscillatory behavior is seen, particularly for intermediate K. Note that position, velocity, and force autocorrelation functions show more pronounced oscillations, as discussed in Appendix D. In Fig. 2(b), we show running integrals over the autocorrelations, defined by as colored lines. According to Eq. (5), I sol FF ðt → ∞Þ= ðk B TÞ ¼ γ, and thus the plateau values observed in Fig. 2(b) for large t reflect the friction coefficient γ. Most importantly, we see that the simulation prediction for γ significantly depends on K, which is somewhat unexpected, since the hydration structure around methane is independent of K, as witnessed by the RDF in Fig. 1

A. Memory function
To gain insight into this puzzling finding, we extract the memory function from simulation data. For this we introduce a variant of existing methods [4,16,[27][28][29] that gives robust and reliable results for kernels and friction coefficients over a wide range of confinement potential strengths. We first observe that in the frozen limit, K ¼ ∞, Eq. (4) predicts that C sol FF ðtÞ ¼ k B TΓðtÞ; i.e., in this limit, the solvent force autocorrelation function equals the memory function. We fit C sol FF ðtÞ in Fig. 2(a) for K ¼ ∞ by a sum of stretched exponentials, with n ¼ 2, shown as a black dotted line. To extract ΓðtÞ also for K ≠ ∞, we calculate C sol FF ðtÞ and I sol FF ðtÞ numerically based on Eqs. (4) and (6) using the functional form Eq. (7) and extract the parameters by simultaneous fits (black broken lines) to the simulation data (colored lines) in Figs. 2(a) and 2(b); the agreement between fits and data is perfect.
The resultant memory kernels are presented in Fig. 3(a). They show a fast decay at 50 fs and a long-time tail extending to about 5 ps. In between, a pronounced and quite abrupt shoulder at a decay time of about 100 fs is present for intermediate K values, qualitatively similar to previous results for the diatomic kernel spectrum in Lennard-Jones fluids [27]. The fitted decay times τ i are shown in Fig. 3(b); note that for intermediate K values we use n ¼ 3 stretched exponentials for the fit (see SM [18] for details on the fit procedure, analysis of the robustness of the fit results, and a comparison with alternative methods). The shortest decay time τ 1 reflects water-methane repulsive interactions: from a fit to the short-distance part of the RDFs in Fig. 1(b) we estimate the fastest relaxation time as τ 1 ¼ 58 fs for K ¼ 0 and τ 1 ¼ 85 fs for K ¼ ∞ (see Appendix A for details), in good agreement with the τ 1 data and indicated in Fig. 3(b) by red arrows. The intermediate K-dependent decay time τ 2 agrees quantitatively with the harmonic oscillation period τ 0 ¼ 2π ffiffiffiffiffiffiffiffiffiffi m=K p , shown as a broken line. The longest decay time τ 3 is rather independent of K and similar to the hydrogen-bond breakage time τ HB ≈ 1.4 ps [30]; we thus associate this time with an intrinsic water relaxation time. We see in Fig. 3(a) that the memory kernel ΓðtÞ changes with K only in the finite K range for which τ 1 < τ 2 ≈ τ 0 < τ 3 , i.e., when the harmonic oscillation period is between the shortest and longest memory relaxation times. The friction kernel thus results from the dynamic interplay of the methane harmonic oscillations, characterized by the relaxation time τ 2 and governed by the external potential strength K, with the solvent relaxation modes in the relaxation time window between τ 1 and τ 3 . Based on this insight we develop a dynamic scaling argument for the behavior of large solutes in generic viscous solvents; see Sec. V.
While the fitted exponents α 1 and α 3 are not too different from unity and thus can be thought of as representing single-exponential relaxation modes, the exponent α 2 is of the order of 10, which indicates strongly nonlinear relaxation (see SM [18] for details). It is this high value of α 2 which leads to the abrupt drop of the memory function at the time scale τ 2 visible in Fig. 3(a). Note that the hydrodynamic power-law tail predicted from continuum hydrodynamics [31], recently observed in optical trap experiments [12] and MD simulations of a supercritical LJ fluid [32], makes a negligible contribution to ΓðtÞ in the friction-relevant ps time scale; see SM [18].

B. Friction coefficient
The results for γ, obtained by the integral over the fitted ΓðtÞ, are presented in Fig. 4 for the three simulation box  (7) as a function of K. Statistical errors are smaller than the symbol size. The two time scales for the frozen limit K ¼ ∞ are denoted by horizontal colored lines on the right. The harmonic oscillation period τ 0 ¼ 2π ffiffiffiffiffiffiffiffiffiffi m=K p is shown as a dashed line, the relaxation times in the K ¼ 0 and K ¼ ∞ limits estimated from methane-water repulsive interactions are denoted by red arrows. sizes considered and agree nicely with the limiting results for K ¼ ∞ (obtained from simulations of a completely frozen methane) and K ¼ 0 (obtained from mean-square displacements for a freely diffusing methane), as indicated by colored bars. The deviation of the K ¼ 0 value from the experimental result γ ¼ 2.2 × 10 −12 kg=s [33] reflects that the SPC/E water viscosity is smaller than the experimental value. From the close agreement between the data for L ¼ 4.5 nm and L ¼ 6.0 nm we conclude that hydrodynamic finite-size effects [34] are negligible for L ≳ 4.5 nm. Figure 4 demonstrates that γ significantly depends on K and reaches at strong confinement (K → ∞) a value about 1.6 times the free-diffusion limit (K ¼ 0). These results are in qualitative agreement with previous simulations that demonstrated the friction of frozen ions to be larger than their free-diffusion values [35]. Note, however, that earlier simulation studies suggested frozen ions to be characterized by a smaller friction coefficient compared to free ions, a disagreement that was never discussed or settled [36,37]. Together with our simulation results for a hydrated water molecule in an external potential, which are shown in Appendix B and exhibit similar effects as for methane, we conclude that friction modification in confinement is thus a general phenomenon that applies to ions, nonpolar as well as polar uncharged molecules. It transpires that simulations of frozen or confined molecules cannot be used to estimate the friction coefficient or the memory function of free molecules, since-as we show here-the memory function depends crucially on the precise strength of the confinement potential. The increase of γ with K, at first sight surprising in light of the decrease of τ 2 with K in Fig. 3(b), is solely caused by the prefactor A 3 of the slowest kernel contribution, as shown in Appendix C. Together with our interpretation that the slowest time scale τ 3 is related to an intrinsic water relaxation time that reflects the hydrogen-bond breakage dynamics, this means that the external confinement potential mostly modifies the hydration-shell contribution to the total friction.
We include typical values of K for van der Waals (vdW), hydrogen, ion, and covalent bonds (estimated in the SM [18]) at the top of Fig. 4, and we see that the most drastic change of γ occurs in the range K ≈ ð10 2 -10 4 Þ kJ=ðmol nm 2 Þ, which matches the strength of typical noncovalent bonds. We also include in Fig. 4 an empiric fit, with the fit parameters γ ∞ ¼ 2.67 × 10 −12 kg=s, γ 0 ¼ 1.73 × 10 −12 kg=s, K 0 ¼ 1 234 kJ=ðmol nm 2 Þ, and c ¼ 2=3 for the system with L ¼ 4.5 nm. This simple function will be useful for comparison with other simulations and experimental data. The parameters K 0 and c correspond to the midpoint and the width of the potential strength range over which the friction changes with K. The small value of the exponent c ¼ 2=3 reflects a rather broad range over which γ changes with K. Interestingly, the value K 0 corresponds to a methane mean-square displacement ffiffiffiffiffiffiffiffi comparable to the width of the first hydration shell in Fig. 1(b).

C. Hydration-shell dynamics
As mentioned in Sec. I, based on NMR experiments [7,8] the friction coefficient of molecular solutes and the hydration-shell dynamics have been shown to be coupled. In this section, we explore the influence of an external potential acting on a solute on the hydration-shell dynamics and in particular check whether the change of the friction coefficient with K is paralleled by a change of the hydration water relaxation dynamics. In Fig. 5(a), we plot the water escape time τ esc1 [38], defined as the mean first passage time to reach the second hydration shell at water-methane separation r ¼ 0.65 nm starting from the first hydration shell at r ¼ 0.37 nm [see Fig. 1(b) for a graphical definition] as a function of K. The averages are calculated from typically 50 000 escape events per parameter combination, from which we estimate the relative error to be less than 1%. τ esc1 has the limiting values τ esc1 ¼ 8 ps for K ¼ 0 and Fig. 4 is paralleled by a slow-down of the escape dynamics of water molecules from the first hydration shell.
Water mean escape times and water translational relaxation times around molecular probes can experimentally be probed by Overhauser dynamic nuclear polarization techniques [39,40]. Obviously, the increase of τ esc1 with K has a trivial contribution due to the fact that the relative distance coordinate relevant for the escape dynamics is governed by the sum of the methane and water translational diffusion constants [39]. To see this, we neglect for a moment that the methane position and the methane-water separation are dynamically coupled, and approximate the relative diffusion constant by the sum of the bulk translational diffusion constants D trans rel ¼ D trans CH 4 þ D trans H 2 O . In the completely confined case, corresponding to K ¼ ∞, the methane diffusion constant D trans CH 4 is obviously zero and thus D trans rel ¼ D trans H 2 O . Approximating the relative diffusion constant to be independent of the methane-water separation (which we know from simulation studies on the relative motion of two water molecules not to be entirely true [41]), the reciprocal relationship between reaction times and diffusion constant, τ esc1 ∼ 1=D trans rel , leads to the prediction From the simulation values D trans CH 4 ¼ ð2.38 AE 0.01Þ × 10 −5 cm 2 =s (see SM [18]) and D trans H 2 O ¼ ð2.70 AE 0.02Þ × 10 −5 cm 2 =s (calculated from the mean-square displacement of a single water molecule, taken from a 250-ns-long simulation of pure water in a L ¼ 4.5 nm box) we thus obtain the estimate τ esc1 ðK ¼ ∞Þ=τ esc1 ðK ¼ 0Þ ¼ 1.88, which is considerably smaller than the ratio τ esc1 ðK ¼ ∞Þ= τ esc1 ðK ¼ 0Þ ¼ 18 ps=8 ps ¼ 2.25 extracted from the simulation data in Fig. 5(a). Is the deviation of the simulated escape time ratio in Fig. 5(a) from the simple estimate associated with a change of the hydration water viscosity (and thus a change of D trans H 2 O ) in the first hydration shell with rising confinement potential strength K?
To look into this, in Fig. 5(b) we show the mean escape time τ esc2 for water molecules starting from the first hydration shell at r ¼ 0.37 nm and reaching the much higher distance r ¼ 1.37 nm, for which we obtain values from τ esc2 ¼ 60 ps for K ¼ 0 to τ esc2 ¼ 110 ps for K → ∞. The mean escape times are higher by about an order of magnitude compared to the results for τ esc1 in Fig. 5(a) and the change of τ esc2 with K is shifted to considerably lower values of K. The ratio of escape times from the completely immobilized and the free methane is now τ esc2 ðK ¼ ∞Þ=τ esc2 ðK ¼ 0Þ ¼ 1.83 and thus considerably smaller than the ratio obtained in the first hydration shell and quite close to our estimate based on the bulk diffusion constants ðD trans , as expected: The methane confinement does not seem to significantly influence the water dynamics beyond the first solvation shell. We conclude that while the main contribution to the translational slow-down of water in the first hydration shell around a methane with rising K comes from the trivial shift of the relative diffusion constant, a significant contribution to this slow-down comes from the modification of the water viscosity in the first hydration shell that accompanies the increase of the methane friction coefficient with rising K. This relates well with our finding that the change of γ with increasing K is predominantly caused by a variation of the long-time contribution to the memory function, as shown in Appendix C, which in turn is related to the intrinsic water dynamics.
The influence of the confinement potential on the hydration-shell dynamics is more directly reflected by the water orientational dynamics around methane, since here we do not need to consider a relative coordinate. In Fig. 6(a), we show the orientation correlation function.
of water molecules with an initial separation below 0.4 nm around a free and a frozen methane molecule. Following previous studies [42], uðtÞ denotes the normalized vector along each OH bond of a water molecule and P 2 ðxÞ ¼ ð3x 2 − 1Þ=2 is the second Legendre polynomial. We include exponential fits C rot ðtÞ ∝ expð−t=τ rot CH 4 Þ for t > 2 ps as dashed lines. The orientational relaxation times are τ rot CH 4 ðK ¼ 0Þ ¼ 3.1 ps and τ rot CH 4 ðK ¼ ∞Þ ¼ 3.6 ps for free and frozen methane and thus differ by a factor of τ rot CH 4 ðK ¼ ∞Þ=τ rot CH 4 ðK ¼ 0Þ ¼ 1. 16. The orientational relaxation time is inversely related to the rotational diffusion constant, τ rot CH 4 ∼ 1=D rot , in the overdamped limit. We thus conclude that the rotational diffusion of water in the first solvation shell around a methane molecule slows down considerably as the translational methane motion is progressively inhibited by the increasing confinement potential. In Fig. 6(b), we show orientation correlation functions of water molecules in the first hydration shell around a free and a frozen water molecule (note that in the frozen case both position and orientation of the central water molecule are fixed). The fitted relaxation times are τ rot H 2 O ðK ¼ 0Þ ¼ 2.5 ps (in agreement with literature values [42]) and τ rot H 2 O ðK ¼ ∞Þ ¼ 3.8 ps and thus differ by a factor τ rot Not surprisingly, due to the formation of directional hydrogen bonds, a frozen water molecule slows down the orientational hydration water dynamics much more than a frozen methane molecule.

V. DISCUSSION AND CONCLUSION
In conclusion, the friction coefficient γ of a methane molecule is shown to significantly depend on the confinement potential strength K, which constitutes a generic and unexpected modification of Stokes's law γ ¼ 6πηR. This reflects, on the one hand, that friction coefficients of fixed solutes differ from free solutes, as suggested previously [43]; on the other hand, it means that freesolution friction coefficients and memory functions cannot be obtained from confined or frozen simulations, contrary to common practice. More generally, our results demonstrate that friction and diffusion on the molecular scale result from the intricate dynamic coupling of solute and hydration degrees of freedom.
Interestingly, the maximal variation of the friction coefficient γ with K we find is similar to the change in γ of a sphere as one goes from a slip (characterized by γ ¼ 4πηR) to a stick situation (in which case one recovers the standard result γ ¼ 6πηR). Clearly, the similarity of this variation of γ to our results is purely coincidental.
For more complicated potentials involving multiple local minima and barriers, for example, for proteins in a suitably defined folding landscape, we speculate that in analogy to our results obtained for methane in a harmonic potential, a local free-energy minimum would produce a local increase in the conjugate friction landscape and, conversely, a free-energy barrier would tend to reduce the local friction. Indeed, this might explain certain universal features seen in diffusivity landscapes extracted from water-explicit simulation trajectories of simple proteins [10]. Along the same lines, the diffusivity profile of a water molecule as a function of the distance from a planar wall indicates a reduced diffusion at local minima of the free-energy profile [44].
The limit of an infinitely strong confinement potential K → ∞ is equivalent to a solute with an infinite mass m → ∞, which directly follows from the fact that the solute is at rest in both cases. By analogy with our results for the friction coefficient γ in Fig. 4, where a continuous increase with rising K is observed, we would expect that γ also goes up continuously as m increases. This is indeed confirmed by experimental measurements of the diffusion coefficient of different isotopes of molecules (such as 13 CO 2 = 12 CO 2 ) and atoms (such as 3 He= 4 He) in water [45].
The methane molecule in our simulations is represented as a simple Lennard-Jones sphere, so the friction increase and the hydration-shell dynamics slow-down with rising K presumably is a rather universal effect that should hold for other molecular solutes as well. Indeed, the effect has been seen for ions in previous simulations [35][36][37] and in Appendix B is demonstrated for a hydrated water molecule in an external potential, which serves as a simple model for a confined polar molecule.
Our simulation setup closely resembles optical trap experiments, where particles that are dispersed in aqueous solution are confined in laser-light-induced harmonic potentials [12]. The lower size limit of the trapped particles has reached the 10 nm scale [13], which is still substantially larger than the size of a methane molecule used in our simulations. Anti-Brownian electrokinetic trapping techniques allow us to trap nanometer-sized molecules [14], but the confinement potential strengths are rather weak. So the legitimate question is whether the effects we predict survive for larger particles, for which a time scale discussion is needed.
The time scale above which a particle in a viscous solvent starts to feel frictional effects is on the scaling level given by the inertial time scale τ m ≡ m=γ (where m denotes the solute mass), at which the particle behavior crosses over from ballistic to diffusive [46]. Using the particle density ρ, particle radius R, and the Stokes friction γ ¼ 6πηR, we obtain τ m ∼ ρR 2 =η, from which we see that the inertial time scale increases strongly with particle radius R. An external harmonic potential of large strength K gives rise to an oscillation time τ 0 ≃ 2π ffiffiffiffiffiffiffiffiffiffi m=K p ; for small K, in the overdamped regime, the time scale is rather τ γ ≃ 2πγ=K. The crossover between the oscillatory inertial and the overdamped regimes occurs at a critical potential strength K Ã ≃ γ 2 =m ∼ η 2 =ρR, where one has τ 0 ∼ τ γ ∼ τ m . Using the notation we introduced when discussing our simulation results, we assume that the memory function is characterized by the longest and shortest relaxation times τ 3 and τ 1 , respectively; in between these times there could in principle be a whole spectrum of intermediate relaxation modes. For water-solvated methane, we have argued that τ 3 is an intrinsic relaxation time of the solvent and thus independent of solute properties (in principle, τ 3 could also depend on internal solute relaxation modes for more complex, large solutes, which however does not change our argumentation and thus is not explicitly considered). The time scale τ 1 stems from the fast relaxation of solutesolvent interactions, i.e., the Lennard-Jones repulsion. Besides a trivial dependence via the relative mass of a solute-water pair (see Appendix A for details), also τ 1 will be rather insensitive to the solute size. In Fig. 7(a), we show the potential time scales τ 0 and τ γ as a function of K for a large particle, for which the inertial time scale τ m is larger than the longest solvent relaxation time τ 3 . In this case the particle oscillates in an inertial, frictionless fashion over the entire solvent relaxation time range τ 1 < t < τ 3 and we do not expect any dynamic coupling between the particle and the solvent for any value of K; in other words, an external potential does not influence the particle friction and all the effects we discuss in this paper are absent.
In Fig. 7(b), we show the opposite situation where the inertial time scale τ m is smaller than the longest solvent relaxation time τ 3 . In the example shown in Fig. 7(b), the particle dynamics changes from inertial to diffusive within the solvent relaxation time range and in this case we do expect the particle friction γ to depend on the potential strength K. This is the situation we encountered in our simulations of methane in water and for which the memory kernel relaxation times are depicted in Fig. 3(b). Note that we draw the potential time scale in Fig. 7(b) in the time range τ m < t < τ 3 as a dotted line, since it is not clear whether in this time range it follows the oscillatory or the overdamped prediction; in fact, the simulation results in Fig. 3(b) seem to follow the oscillatory prediction τ 0 ≃ 2π ffiffiffiffiffiffiffiffiffiffi m=K p in the entire range τ 1 < t < τ 3 . To put in explicit numbers, with the methane mass m CH 4 ¼ 16u ¼ 2.7 × 10 −26 kg and the methane diffusion constant D trans CH 4 ¼ 2.38 × 10 −5 cm 2 =s, we obtain for the inertial time scale τ m ¼ m CH 4 D trans CH 4 =k B T ¼ 15 fs, which is even smaller than the fastest water relaxation time τ 1 ¼ 58 fs for K ¼ 0 or τ 1 ¼ 85 fs for K ¼ ∞.
It transpires that the coupling scenario depicted in Fig. 7(b) can also for a large particle be obtained if the longest solvent relaxation time τ 3 is increased accordingly. This can be achieved by the use of highly viscous solvents, such as glycerol or polymer solutions. For entangled polymer solutions, the longest relaxation time scales as a power law with the polymer length and can thus be increased straightforwardly [47]. In essence, in order to be able to observe the K-dependent friction effects we describe in this paper also with large particles, which can be easily confined in optical traps, one has to sufficiently increase the solvent viscosity which will lift the upper solvent relaxation time τ 3 and at the same time bring down the inertial time scale τ m ≡ m=γ.
A fundamentally different experimental consequence of our results concerns the hydration-shell dynamics around FIG. 7. Dynamic scaling diagrams for (a) a large solute particle and for (b) a small particle. In (a) the inertial time scale τ m lies above the longest water relaxation time τ 3 and the particle relaxation time in the potential (shown as a solid line) changes from the harmonic oscillator prediction τ 0 ≃ 2π ffiffiffiffiffiffiffiffiffiffi m=K p (for K > K Ã ) to the overdamped prediction τ γ ≃ 2πγ=K (for K < K Ã ) above τ 3 . In (b) τ m is smaller than τ 3 and the potential relaxation time becomes overdamped within the relaxation time range. In this case we expect the external potential of strength K to modify the particle friction. probe molecules, which is measured by a number of different experimental methods such as combined 2 H-17 O nuclear relaxation [7,8], nuclear Overhauser effect [48], dynamic Stokes shift [49], Overhauser dynamic nuclear polarization [39], fluorescence [50], 2D infrared spectroscopy [51], quasielastic neutron scattering [52], and THz absorption [53,54]. The water hydration dynamics around methane slows down significantly above potential strength values that correspond to covalent bonds, this is seen for the translational water motion in Fig. 5, although here the trivial shift of the effective relative diffusivity contributes (as discussed in Sec. IV C), as well as for the orientational water dynamics in Fig. 6. These results thus suggest that hydration water relaxation times not only depend on the type of probe molecule, i.e., polar versus nonpolar [55], but also on how tightly the probe molecule is confined or bound to a macromolecule.
The translational diffusion of hydration water relative to probes attached to proteins or lipid assemblies is experimentally observed to be about 3.5-5 times slower than around probes that freely diffuse in bulk water [39,52,56], for probes attached to DNA segments a slow-down by only a factor of about 2 was reported [40], which seems to be the lower bound for this type of measurements. A similar slowdown of the orientational correlation time of hydration water at the protein water interface by a factor of about 2 compared to bulk was reported based on experiments [57] and simulations [42]. Clearly, when a probe molecule is anchored at a larger molecule, not only is its motion confined by an effective potential whose stiffness depends on the elasticity of the anchoring group, also the environment around the probe molecule is modified due to the presence of the linker groups and the scaffold to which the probe is attached, which for typical systems has been shown to be the main factor determining the hydration water dynamics around the probe [42,58,59]. Nevertheless, we argue that the confinement-induced friction modification we demonstrate in this paper will certainly contribute to the water hydration slow-down-among other factorsand thus is a noteworthy mechanism.
In fact, in two experiments the flexibility of the scaffold onto which a probe molecule was anchored was varied without a drastic change in the environment of the probe. In a site-specific femtosecond-resolved fluorescence study using 16 tryptophan labeled myoglobin mutants, the water relaxation time was shown to be correlated with the local protein structural rigidity [50]. Based on THz absorption studies, the hydration-shell relaxation dynamics was argued to become faster with increasing protein flexibility; in these experiments, the protein flexibility was modified by suitably chosen mutations [53]. The range over which a protein perturbs the hydration-shell dynamics according to the analysis of THz experiments extends to several hydration layers, whereas in our analysis of translational water motion we only see an effect of the confining potential in the first hydration layer. Nevertheless, we conclude that these two experimental studies suggest that water hydration dynamics is coupled to the local rigidity of the probe molecule, even when the probe molecule is attached to a large protein. This is fully confirmed by simulations that show a significant slow-down of hydration shell dynamics around a frozen protein compared to a flexible protein of the same structure [60].

ACKNOWLEDGMENTS
Financial support from the Deutsche Forschungsgemeinschaft (DFG) via Grant No. SFB 1114 is acknowledged. We are pleased to acknowledge helpful discussions with Bertil Halle and Bill Eaton.

APPENDIX A: ESTIMATE OF THE CONTRIBUTION OF THE REPULSIVE INTERACTION BETWEEN METHANE AND WATER TO THE MEMORY KERNEL
Here we show that repulsive methane-water interactions give rise to a relaxation time that matches the shortest memory time scale τ 1 . The force autocorrelation function of an undamped harmonic oscillator is where the characteristic frequency is and denotes the relative mass of the molecule pair. We determine the force constant K rel from the free energy F ðrÞ as a function of the methane-water distance r, shown by a black line in Fig. 8, which is obtained from the radial distribution function gðrÞ shown in Fig. 1 Note that F ðrÞ is very anharmonic even on energy scales of the order of the thermal energy k B T, so a harmonic fit to the entire function F ðrÞ is not useful. To extract the fastest time scale of particle motion in this free-energy profile, we fit a harmonic potential, to the repulsive part in the range 0.32 < r < 0.36 nm only, which is denoted by the red line in Fig. 8. The result is K rel ¼ 3; 548 kJ=ðmol nm 2 Þ. To obtain the fastest time scale τ 1 of the memory kernel defined in Eq. (7), we demand that the stretched exponential function in Eq. (7) and the oscillatory function in Eq. (A1) have both decayed to 1=e, which is equivalent to from which we obtain with our estimate for K rel and m rel given by Eq. (A3) the time scale which holds for unconfined methane. If the methane molecule is frozen, which corresponds to m CH 4 → ∞, we use m rel ¼ m H 2 O and obtain These times are included in Fig. 3(b) as red arrows and match the simulated times quite nicely. For comparison, if we repeat the analysis but use instead of the fit to the repulsive part of the free energy a harmonic fit around the minimum of the free energy F ðrÞ, denoted by the blue line in Fig. 8, we obtain the time scales τ 1 ¼ 119 fs for free methane and τ 1 ¼ 174 fs for frozen methane. These time scales are significantly larger than the fastest time scale of the memory function. We conclude that the fast initial decay of the memory function is indeed caused by repulsive interactions between the solute and the solvent molecules.

APPENDIX B: CONFINEMENT-DEPENDENT FRICTION FOR A HYDRATED WATER MOLECULE
To demonstrate that the confinement dependence of the friction coefficient is not limited to hydrophobic solutes, we simulate a confined SPC/E [24] water molecule solvated in a cubic SPC/E water box with L ¼ 4.5 nm for 250 ns for each K. We present in Fig. 9 the solvent force autocorrelation functions C sol FF ðtÞ and the running integrals I sol FF ðtÞ of the confined water molecule for different potential strengths K between 25 and 250 000 kJ=ðmol nm 2 Þ. As shown in Sec. II, the height of the plateau of I sol FF ðtÞ corresponds to the friction constant; i.e., We conclude from the K dependence of the heights of the plateaus of I sol FF ðtÞ shown in Fig. 9(b) that also the friction coefficient of water increases significantly with rising confinement potential strength.

APPENDIX C: DECOMPOSITION OF THE FRICTION COEFFICIENT
Here, we decompose the methane friction coefficient into contributions from different terms in the memory kernel according to γ ¼ P i γ i with the definition FIG. 8. The free energy F ðrÞ associated with the distance r between the methane and surrounding water molecules together with a harmonic fit to the repulsive part (red line) and a harmonic fit around the free energy minimum (blue line).  where A i , α i , and τ i are the parameters of the memory kernels defined in Eq. (7). The contributions γ i are shown as a function of K for the system with L ¼ 4.5 nm in Fig. 10 together with the total friction γ and the sum of the short-time contributions γ 1 þ γ 2 . The results demonstrate that the total friction coefficient is dominated by the longtime contribution γ 3 , whereas the sum of the short-time contributions γ 1 þ γ 2 is constant. Since the time scale τ 3 is rather constant as a function of K, as shown in Fig. 3(b), we conclude that the change of γ with K is solely due to an increase of the amplitude A 3 with increasing K.

APPENDIX D: CORRELATION FUNCTIONS
AND ALTERNATIVE METHODS

Correlation functions for methane
In Fig. 11, the (normalized) autocorrelation functions C FF ðtÞ ¼ hFð0ÞFðtÞi, C _ x _ x ðtÞ ¼ h_ xð0Þ_ xðtÞi, and C xx ðtÞ ¼ hxð0ÞxðtÞi are presented for the medium system size (L ¼ 4.5 nm) and five different harmonic spring constants K. We observe that with larger spring constant all functions are subject to pronounced oscillatory behavior with a frequency close to ω 0 ¼ ffiffiffiffiffiffiffiffiffiffi K=m p . The amplitude of these oscillations decays on time scales between 0.5 and 5 ps, force and velocity autocorrelation functions with small spring constants decay the fastest. Note that these oscillations are much more pronounced than the oscillations in C sol FF ðtÞ presented in Fig. 2(a).

Analytic expressions for correlation functions
Here, we give analytic expressions for the correlation functions C xx ðtÞ ¼ hxð0ÞxðtÞi, C _ x _ x ðtÞ ¼ h_ xð0Þ_ xðtÞi, and C FF ðtÞ ¼ hFð0ÞFðtÞi, which can be derived analogously to Eq. (4) from the generalized Langevin equation Eq. (1): As we haveC _ x _ x ð0Þ ¼ 0 andC FF ð0Þ ¼ 0 for K ≠ 0, clearly, the integrals of these two correlation functions will not be suited to extract the friction constant. Only the position autocorrelation function is nonzero at zero frequency, i.e., C xx ð0Þ ¼ k B TΓð0Þ=K 2 ¼ 2k B Tγ=K 2 for K ≠ 0, and thus in principle allows us to extract the friction coefficient. We remark that for vanishing confinement K ¼ 0, the velocity autocorrelation fulfills  x _ x ðtÞ=C _ x _ x ð0Þ, and C xx ðtÞ=C xx ð0Þ for the system with the size L ¼ 4.5 nm and five different spring constants K of the harmonic potential. and thus allows us to extract γ. Note, however, that all given correlation functions are subject to large amplitude oscillations for sufficiently high values of K, which limits the practical usefulness, as is demonstrated in Fig. 11.

Comparison to an extraction method based on the position autocorrelation function
Hummer et al. developed a method to calculate the friction coefficient from the position autocorrelation function in systems with a biasing potential [61]. Starting from the work of Woolf and Roux [62], they showed that the diffusion coefficient D of a particle in a harmonic potential UðxÞ ¼ Kx 2 =2 can be written as or, in terms of the friction coefficient, We note that this result can be directly derived from Eq. (D1) according to which is equivalent to Eq. (D6). We compare the (normalized) position autocorrelation function C xx ðtÞ ¼ hxð0ÞxðtÞi with the (normalized) solvent force autocorrelation function C sol FF ðtÞ ¼ hF sol ð0ÞF sol ðtÞi in Fig. 12. We conclude that for sufficiently small values of K, the two methods should work equally fine. However, as already pointed out in Ref. [61], the position autocorrelation function is dominated by oscillations, which render the integration numerically unstable for large values of K. Our method that is based on C sol FF ðtÞ is suitable to extract the friction coefficient and the memory kernel for all values of the confinement strength K.

Parametrization-free methods to extract the memory kernels
In principle, the memory functions can be extracted numerically in a parametrization-free way from correlation functions in the time domain [3,16,17]. Because of the instability of the inversion scheme [3], this approach is in practice quite unstable for systems with high values of K, which are subject to fast oscillatory behavior in the relevant correlation functions (see Fig. 11). We remark that the parametrization of the memory kernel as it is used in this work can be regarded as a regularization of the problem [3].