{'fraction_non_alphanumeric': 0.032418234807325585, 'fraction_numerical': 0.012485158133342928, 'mean_word_length': 4.97078410311493, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse Sampling of Water Density Fluctuations in Interfacial Environments

The free energetics of water density fluctuations near a surface, and the rare low-density fluctuations in particular, serve as reliable indicators of surface hydrophobicity; the easier it is to displace the interfacial waters, the more hydrophobic the underlying surface. However, characterizing the free energetics of such rare fluctuations requires computationally expensive, non-Boltzmann sampling methods like umbrella sampling. This inherent computational expense associated with umbrella sampling makes it challenging to investigate the role of polarizability or electronic structure effects in influencing interfacial fluctuations. Importantly, it also limits the size of the volume, which can be used to probe interfacial fluctuations. The latter can be particularly important in characterizing the hydrophobicity of large surfaces with molecular-level heterogeneities, such as those presented by proteins. To overcome these challenges, here we present a method for the sparse sampling of water density fluctuations, which is roughly two orders of magnitude more efficient than umbrella sampling. We employ thermodynamic integration to estimate the free energy differences between biased ensembles, thereby circumventing the umbrella sampling requirement of overlap between adjacent biased distributions. Further, a judicious choice of the biasing potential allows such free energy differences to be estimated using short simulations, so that the free energetics of water density fluctuations are obtained using only a few, short simulations. Leveraging the efficiency of the method, we characterize water density fluctuations in the entire hydration shell of the protein, ubiquitin; a large volume containing an average of more than six hundred waters.


Introduction
An understanding of density fluctuations in bulk water and at interfaces, has played a central role in the description of hydrophobic effects, 1,1-12,14-20 which drive biomolecular [21][22][23][24][25][26] and other aqueous assemblies. [27][28][29][30][31] For example, the fact that P v (N), the probability of observing N water molecules in a small observation volume (v 1 nm 3 ), obeys Gaussian statistics, 3,32 has provided molecular-level insights into the the pressure-induced denaturation of proteins, [33][34][35] as well as the convergence of protein unfolding entropies at a particular temperature. 19,32 Similarly, the Lum-Chandler-Weeks (LCW) theory 5 prediction that P v (N) should develop fat low-N tails for large volumes (v 1 nm 3 ) in bulk water, and near hydrophobic surfaces, 5,15,36 and its verification by simulations, 1,14,18,37 has clarified that water near hydrophobic surfaces is sensitive to perturbations, 18,[38][39][40][41][42][43] and led to the insight that extended hydrophobic surfaces could generically serve as catalysts for the assembly and disassembly of small hydrophobic solutes. 14,44 Importantly, the fact that low-N fluctuations are enhanced near hydrophobic surfaces, also makes them a suitable metric for quantifying the hydrophobicity of a surface, or the strength of its interactions with water. 1,16,18,20 Fluctuations are particularly useful in characterizing the hydrophobicity of complex surfaces with molecular-level heterogeneities, wherein conventional macroscopic measures, such as the water droplet contact angle, break down. 45 Indeed, for the rugged, heterogenous surfaces of proteins, hydrophobicity depends not only on the chemistry of the underlying residues, 12,18,20,46 but also on the particular topography [47][48][49] and chemical pat-tern 12,[50][51][52][53] presented by the protein, and can depend non-trivially on the specific combination of the two. 39,40,[54][55][56] Fluctuations have previously been used to characterize the hydrophobicity of such complex surfaces using small (e.g., methane-sized) v, 12,46 wherein the fluctuations of interest can be estimated using equilibrium molecular simulations. 3,57 However, the likelihood of low-N fluctuations decreases roughly exponentially with the size of v, and estimating P v (N) for larger v requires computationally expensive biased sampling techniques, such as umbrella sampling. 1,2,58 In particular, the larger the volume of interest, the larger the number of biased simulations required.
As a result, using large volumes, such as the hydration shells of entire proteins, to characterize surface hydrophobicity would be very expensive. Further, using umbrella sampling in conjunction with more computationally expensive treatments, which incorporate the polarizability of interfacial waters or electronic structure effects, 60,61 to characterize density fluctuations in large volumes, would also be prohibitively expensive.
Here we present a method that enables estimation of P v (N) at a number of sparsely separated N-values, albeit with a roughly two orders of magnitude increase in computational efficiency, as compared with conventional techniques, such as umbrella sampling or free energy perturbation.
To facilitate sparse sampling, we employ thermodynamic integration to estimate free energy differences between biased ensembles, thereby circumventing the umbrella sampling requirement that adjacent biased distributions overlap. Building upon recent work by Patel and Garde, 20 we employ a linear biasing potential, which enables efficient estimation of free energy differences between biased ensembles using averages that converge rapidly and require only short simulations. In the following section, we first derive the central equations underlying our method, followed by details pertaining to the systems studied and the simulations employed. In the subsequent section, we demonstrate how the method works using a small volume in bulk water. We then apply the method to characterize the hydrophobicity of interfaces; first, uniform self-assembled monolayer surfaces, and then the entire hydration shell of the protein, ubiquitin. We then discuss the underpinnings of the method's efficiency and its limitations, and conclude with a discussion of scenarios where the method may find broader applicability.

Derivation of the Central Equations
Consider an observation volume, v, described by its size, shape, and location. We are interested in characterizing the statistics of water density fluctuations in v, as quantified by the probability, to enabling sparse sampling of P v (Ñ), the method we present below, can also be used to obtain To estimate P v (Ñ) over the entire range ofÑ-values, including those that are highly improbable, standard approaches such as umbrella sampling, prescribe the application of a biasing potential that is a function ofÑ v . 58,62 Consider the biasing potential, U φ (Ñ v ), parametrized by φ , so that the Hamiltonian of the biased system becomes H φ = H 0 + U φ (Ñ v ); let Q φ be the partition function associated with H φ . The unbiased probability, P v (Ñ), of observingÑ coarse-grained waters in v can be readily related to the corresponding probability in the biased ensemble, P φ v (Ñ), as Here is the free energy difference between the biased and unbiased ensembles. The first term in Equation 2 can be readily estimated from biased simulations, and the second term is trivial to evaluate, because the functional form of the biasing potential is known. In contrast, the third term, F φ , is not known a priori; umbrella sampling relies on overlap in adjacent biased distributions to estimate F φ , usually with the help of algorithms such as WHAM or MBAR. [63][64][65] Such a prescription provides accurate estimates of F v (Ñ) over the entire range of sampledÑ v -values.
Nevertheless, obtaining a continuous F v (Ñ) distribution is not always necessary; instead, a knowledge of its functional form atÑ-values that are well separated can often be sufficient. However, the umbrella sampling requirement of overlap between adjacent biased distributions is incommensurate with the spirit of such sparse sampling. Importantly, the overlap requirement also necessitates a large number of biased simulations, and can be computationally expensive, in particular for large v. The proposed method aims to circumvent this overlap requirement by estimating F φ using thermodynamic integration. We will demonstrate that such an approach can be both accurate and highly efficient for a particular choice of the biasing potential, U φ (Ñ v ). In particular, we build upon recent work by Patel and Garde, 20 which employed a linear biasing potential to estimate the free energy of emptying v, as quantified by F v (Ñ = 0), with a roughly 2 orders of magnitude decrease in computational effort. Here we show how the entire P v (Ñ) distribution can be sparsely sampled with a similar increase in computational efficiency, through the use of such a linear potential, U φ (Ñ v ) = φÑ v , whose strength is determined by the parameter, φ . In the presence of the linear potential, U φ , Equation 2 becomes As discussed above, the first two terms of the above equation are readily determined. Here we determine F φ using thermodynamic integration, that is, by integrating ∂ F φ /∂ φ = Ñ v φ as: The integration is performed over a range of φ -values from zero (i.e., the unbiased system) to the particular φ -ensemble of interest. Increasing φ penalizes the presence of water molecules in v; To illustrate how the method works, consider employing Equation Thus, from estimates of the average values, Ñ v φ k , obtained from a series of biased simulations with potentials of strength, {φ k }, the corresponding estimates of , this method can also be used to obtain estimates of F v (Ñ) in the vicinity ofÑ = Ñ v φ k ; in particular, using Equations 3 and 4, F v (Ñ) can be estimated for eachÑ-value sampled by any of the biased simulations.
We note that while we have chosen to focus on the fluctuations inÑ v , the method presented here can readily be generalized to chararcterize fluctuations in any order parameter. Further, as mentioned above and shown in the Supporting Information, this approach can also be used to indirectly sample discrete order parameters; in particular, P v (N), the probability of observing N water molecules in v, can be readily estimated. Finally, we note that while there are certain similarities between our method and the Umbrella Integration (UI) method, proposed by Kästner and Thiel,3 there are important fundamental and practical differences between the two methods; a detailed discussion of these issues is included in the Supplementary Information.

Simulation Details
We illustrate the application of Equation 5 to compute the F v (Ñ) for 3 different classes of systems: (1) a small spherical volume of radius R v = 0.5 nm in bulk water; (2) a cylindrical disk of width, w = 0.3 nm and radius, R v = 2 nm, in the hydration shell of self-assembled monolayers (SAMs) of alkyl chains with two different head groups: -OH and -CH 3 , which correspond to hydrophilic and hydrophobic SAM-water interfaces respectively; and (3)  To study each of these systems, we use all-atom MD simulations using the GROMACS package, 68 suitably modified to incorporate the biasing potentials of interest. We perform biased simulations at different φ -values, and use our central Equation 5 to estimate F v (Ñ) at sparsely separated N-values. In all cases, a simulation box with periodic boundary conditions was employed, and the leap frog algorithm 62 with a 2 fs time-step was used to integrate the equations of motion. For the bulk water and SAM systems, the SPC/E model of water was used, 69 whereas the protein was hydrated using TIP4P water. 70 The short-range Lennard-Jones and electrostatic interactions were truncated at 1 nm, whereas the long range electrostatic interactions were treated using the Particle Mesh Ewald algorithm. 71 Bonds involving hydrogen atoms were constrained; the SHAKE algorithm 72 was employed for the bulk water system, and the LINCS algorithm 73 was used for the SAM and protein systems. In all cases, a constant temperature of 300 K was maintained by using the canonical velocity-rescaling thermostat. 74 For the bulk water system, the pressure was additionally maintained at 1 bar using the Parrinello-Rahman barostat. 75 While the SAM and protein systems are simulated in the canonical ensemble, a buffering water-vapor interface was employed to ensure that the system remains at coexistence pressure. 1,2,76 Additional system specific details including those pertaining to the placement of observation volumes are given below.
Bulk Water: The spherical observation volume is situated at the center of a cubic simulation box with a side of length 6 nm.

SAM Surfaces:
The simulation setup is similar to that used in ref. 14. The disk-shaped observation volumes are placed adjacent to the SAM surfaces, with their axes perpendicular of the surface. For the CH 3 -terminated SAM surface, the observation volume is centered at the first peak of the water density distribution in the direction perpendicular of the SAM surface. For the OH-terminated SAM surface, the proximal edge of the observation volume is placed at a location where the averaged density of the SAM heavy atoms drops below 2% of the bulk water density.
Protein Hydration Shell: We used the CHARMM27 forcefield 77 to simulate ubiquitin, which was hydrated by roughly 13,000 TIP4P water molecules. 70  ps is discarded as equilibration in response to the biasing potential, and the subsequent 100 ps is used to estimate N v φ , whereas the subsequent 2.9 ns were used to estimate δ N 2 v φ and to obtain accurate estimates of P v (N) using umbrella sampling in conjunction with WHAM. To obtain P v (N) using umbrella sampling, several additional biased simulations are also needed to ensure overlap between adjacent biased distributions.

Results and Discussion
Bulk Water readily obtained atÑ = Ñ v φ , and are shown in Figure 1d. The biased distribution free energies, , are small because typical water number distributions are peaked at their means; resulting from the addition of the three terms in Figure 1d, displays a fat low-Ñ v tail as shown in Figure 1e, and is in excellent agreement with the exact F v (Ñ) obtained by umbrella sampling.

SAM -water interfaces
The low-Ñ tail of water density fluctuations, estimated in the vicinity of a surface, is expected to depend strongly on surface hydrophobicity, or the strength of surface-water interactions. 1,5,[9][10][11][12][14][15][16]18,20,37,56 To illustrate that the sparse sampling method can capture such marked differences in the respective F v (Ñ) adjacent to the hydrophilic SAM surface is parabolic (black line) to a good approximation, consistent with the underlying density fluctuations being Gaussian. In contrast, a marked fat tail in water density fluctuations is observed adjacent to the hydrophobic SAM surface, in agreement with previous findings. 1,15,18 Such a difference has previously been demonstrated to arise from the fact that water near hydrophobic surfaces is situated at the edge of a dewetting transition. 18,78 This difference is also reflected in the sensitive response of interfacial water to perturbations, that is, in the sigmoidal decrease in Ñ v φ with increasing φ (Figure 3b), and in a peak in the corresponding susceptibility, (Figure 3d). In contrast, the Gaussian fluctuation adjacent to the hydrophilic surface are associated with linear decrease in Ñ v φ ( Figure 3b) and a constant susceptibility (Figure 3d).
It is clear that there is a correspondence between the functional forms of the free energetics of (c) , also shows distinct features for the two SAM surfaces. While the susceptibility of the OH-terminated SAM surface is roughly constant, a sharp peak is observed in the susceptibility of the CH 3 -terminated SAM surface due to the corresponding sigmoidal response of Ñ v φ to φ .
water number fluctuations, F v (Ñ), and the corresponding response of the average water number, Ñ v φ , to the strength of the potential, φ . In particular, if F v (Ñ) is known over the entire range ofÑ-values of interest, Ñ v φ can be readily obtained at all values of φ through reweighting. 20 Indeed, the curves shown in Figures 3b and 3d were obtained in that manner. In contrast, the sparse sampling method introduced above not only allows us to perform the inverse operation, it does so with only a select few values of Ñ v φ ; that is, given well-separated Ñ v φ estimates at a select few φ -values, the sparse sampling method allows us to estimate F v (Ñ) at thoseÑ = Ñ v φvalues.

Entire Protein Hydration Shell
Water density fluctuations have previously been used to characterize the hydrophobicity of regions within protein hydration shells, in volumes containing a few to roughly ten waters on average.
Here we leverage the efficiency of the sparse sampling method to estimate water density fluctuations in the entire hydration shell of the protein, ubiquitin (Figure 4a), which is significantly larger, and contains roughly 660 waters on average. In contrast with the uniform SAM surfaces studied in the previous section, the protein-water interface is chemically and topographically heterogeneous; such complexity influences the corresponding water density fluctuations in a non-trivial manner. 12,16 Using our central Equation 5, the response, Ñ v φ , to the biasing potential, φÑ v , shown in Figure 4b, can readily be transformed into the free energy, F v (Ñ). As shown in Figure 4c, percolation transitions on partially hydrated protein surfaces. It has been shown that when proteins are hydrated with insufficient waters to form a complete protein hydration shell, an inter-connected network of waters appears abruptly over a narrow range of hydration levels. [80][81][82] When a protein is partially hydrated, its hydration waters are essentially in vacuum; in contrast, the hydration shell waters in our simulations are surrounded by and interact with other waters. However, in both cases, a collective (percolation or dewetting) transition is facilitated by a competition between protein-water and water-water interactions, with the protein providing a heterogeneous surface that displays a wide range of surface chemistries and thereby protein-water interaction strengths.

Efficiency and Limitations of the Method
We note that F v (Ñ), shown in Figure 4c, was obtained using 10 simulations run for 0.
We choose β φ 0 = 1, and to facilitate comparison with the ubiquitin F v (Ñ) (Figure 4), we choose β κ 0 = 10 −8 so that the two distributions have roughly the same mean and variance. As shown in Figure 5a, Ñ v φ for the analytical model features a sharp cliff, wherein the value of Ñ v φ decreases dramatically over a narrow range of φ -values. The corresponding F v (Ñ) obtained using the sparse sampling method is shown in Figure 5b, and highlights the limitations of the method. Due to the sharp cliff in Ñ v φ , and the associated inability to sample intermediateÑ-values, estimates of F v (Ñ) can not be obtained for a large rangeÑ-values. In addition, an inability to capture the precise location of the cliff, leads to substantive errors in F φ for higher φ -values, and thereby in F v (Ñ) for the lowerÑ-values. A sharply decreasing response function akin to that for the analytical model is thus a strong indicator of bistability in the underlying free energetics; in such cases, the sparse sampling method should not be used in conjunction with linear potentials as proposed here.
Strategies for overcoming such limitations, which are likely to be encountered when water is at or close to coexistence with its vapor in the volume of interest (e.g., in hydrophobic confinement), will be the subject of a future publication.

Outlook
By circumventing the umbrella sampling requirement of overlap between adjacent biased distributions, and instead using thermodynamic integration to estimate free energy differences between the biased and unbiased ensembles, the method presented here enables sparse sampling of the free energetics of water density fluctuations, F v (Ñ). Furthermore, a judicious choice of the functional form of the biasing potential, that is, one which is linear in the order parameter of interest,Ñ v , enables estimation of F v (Ñ) in a computationally efficient manner. The low-Ñ behavior of F v (Ñ) serves to characterize the hydrophobicity of complex surfaces with nanoscale heterogeneities, such as those presented by proteins; here, we use the method to characterize F v (Ñ) in a large volume, constituting the entire hydration shell of the protein, ubiquitin. Such a characterization of protein hydration shells not only provides an overall measure of protein hydrophobicity, but also quantifies the free energy required to displace water molecules from the protein hydration shell, which could inform its propensity to bind hydrophobic ligands. 83 In addition to facilitating characterization of F v (Ñ) in large volumes, the efficiency of our method may also enable estimation of F v (Ñ), using more detailed treatments of water in the bulk or at interfaces, which are inherently expensive form a computational standpoint. Such treatments include force fields that explicitly account for molecular polarizability, 60,84 or ab initio molecular dynamics simulations, which incorporate electronic structure effects. 61,85 The method presented here is fairly general, and can be readily generalized to other order parameters besidesÑ v as well as to higher dimensions. A straightforward generalization of the method could facilitate characterization of the free energetics of concentrations (as opposed to number density) fluctuations in multi-component aqueous solutions and mixtures. 86 Biasing potentials that couple linearly to an order parameter of interest are also employed in a variety of other contexts; examples include constant electrostatic potential simulations used to study charge fluctuations in capacitors, 87 simulations that bias trajectory space using a dynamical order parameters such as activity, 88 and alchemical methods for estimating binding free energies. 89,90 While free energy perturbation is typically employed in alchemical calculations to estimate the free energy differences between the ensembles of interest, the method introduced here could additionally provide the statistics of the energy differences between those ensembles; an understanding of such statistics may further inform optimal strategies for accurately and efficiently estimating the corresponding free energies. Such an understanding may also lead to the development of analytical expressions for the estimation of free energies. Due to the dramatic increase in computational efficiency that the method provides, we believe that it will also be well-suited for the characteri-

Derivation of the Umbrella Sampling Equation
Here we derive the umbrella sampling equation used in the main text, for a a biased Hailtonian, . The probability, P v (Ñ), of observingÑ coarse-grained waters in the unbiased ensemble is given by: where Q 0 ≡ e −β H 0 is the partition function of the unbiased ensemble. Employing a straightforward reweighting strategy, we get where Q φ ≡ e −β H φ is the partition function of the biased ensemble. Finally, recognizing that the delta function allows us to pull the e βU φ (Ñ v ) term outside the integral, we get Taking the logarithm of both sides and recognizing that δ (Ñ v −Ñ) φ is simply the probability, P φ v (Ñ), of observingÑ coarse-grained waters in v in the biased ensemble, we get the central result of umbrella sampling as Getting P v (N) using the sparse sampling method In the main text, we develop a sparse sampling method to estimate P v (Ñ), the probability of ob-servingÑ coarse-grained waters in a volume, v, of interest. Here we show that using the same framework, a closely related discrete order parameter, such as the number of waters in v, can also be sparse sampled, that is, the probability distribution, P v (N), of observing N waters in v can also be estimated. As in the previous section, the joint probability distribution, P v (N,Ñ), of observing N waters andÑ corse-grained waters in the unbiased ensemble is given by: The central idea is to bias the system using the coarse-grained order parameter,Ñ, but collect both N and N data from the biased simulations to obtain P φ v (N,Ñ), the biased joint distribution, which is the first term on the right side of Equation S11. The biasing potential (the second term) is known analytically and the free energy differences (the third term) between biased and unbiased ensemble can be readily estimated using thermodynamic integration, as discussed in the main text. Finally, P v (N) can be estimated by an integrating the unbiased joint distribution overÑ, To illustrate Equation S11 and Equation S12, we apply this strategy to the hydration shell of the protein, ubiquitin, which was discussed in the main text. An example of the biased and corresponding unbiased joint distributions at an N-value (in this case, N = N φ with β φ = 4) of interest, are shown in Figures S1a and S1b respectively. Figure S1c shows that the sparse sampled free ener- The corresponding unbiased joint distribution was obtained using equation S11. (c) Each biased simulation was used to sparsely sample F v (N) by integrating the respective unbiased joint distributions using equation S12; the integral of the joint distribution from the β φ = 4.0 ensemble is shown in purple. The resulting F v (N) agrees well with the exact results obtained using the Indirect Umbrella Sampling method. 1,2 Comparison between the Sparse Sampling Method and Umbrella Integration The Umbrella Integration (UI) method, proposed by Kästner and Thiel, 3 utilizes principles of both umbrella sampling and thermodynamic integration to estimate free energies, akin to the sparse sampling method that we present here. However, there are important distinctions between the two methods, which we discuss below.
1) While we explicitly estimate the free energy differences between the biased and unbiased ensembles using thermodynamic integration, UI prescribes eliminating them from the analysis altogether by differentiating Equation 2 of the main text. ∂ F v /∂Ñ v -values are then estimated using biased simulations, and integrated numerically to yield estimates of F v (Ñ v ).
2) The authors focus on the widely-used harmonic form of the biasing potential, and show that in the limiting case of a stiff spring, UI becomes equivalent to thermodynamic integration. In contrast, we use a linear biasing potential, which contributes significantly to the efficiency of our method; in particular, the linear potential ensures that Ñ v φ decreases monotonically with φ (because ∂ Ñ v φ /∂ φ = −β δÑ 2 v φ < 0), and enables estimation of F φ using only a few biased simulations.
3) Finally, as recognized by its authors, UI does not, in principle, require overlap between adjacent windows. In practice, however, the functional form of ∂ F v /∂Ñ v is not known a priori, so accurate estimates of F v (Ñ v ) rely on estimates of ∂ F v /∂Ñ v at finely spacedÑ v -values. Indeed, the authors use UI not to perform sparse sampling, but to illustrate that it leads to smaller statistical errors relative to WHAM (Weighted Histogram Analysis Method).

Fluctuations in the Hydration Shell of the Protein, Hydrophobin
Here, we use the sparse sampling method to estimate P v (Ñ) in the hydration shell of the protein, hydrophobin II (PDB ID: 2B97), which is known to have a large hydrophobic patch on the protein surface. As shown in Figure S2, the protein displays characteristics that are quite similar to that of the protein, ubiquitin, discussed in the main text.