Determination and evaluation of the nonadditivity in wetting of molecularly heterogeneous surfaces

Significance Every folded protein presents an interface with water that is composed of domains of varying hydrophilicity/-phobicity. Many simulation studies have highlighted the nonadditivity in the wetting of such nanostructured surfaces in contrast with the accepted theoretical formula that is additive. We present here an experimental study on surfaces of identical composition but different organization of hydrophobic and hydrophilic domains. We prove that the interfacial energy of such surfaces differs by ∼20% and that a significant difference in the interfacial water H-bonding structure can be measured. As a result, in combination with molecular-dynamics simulations, we propose a model that captures the wetting of molecularly heterogeneous surfaces, showing the importance of local structure (first-nearest neighbors) in determining the wetting properties.


Effect of deuteration on the ligand shell nanoparticles
In the current work, most of the data i.e. SANS, SAXS, FTIR, TEM, SFG, AFM, and Contact angle, were measured on the same batches of nanoparticles, i.e. NPs coated with dMPA and OT. For the measurements of NMR and TGA, we used NPs coated with MPA and OT due to the requirements of the techniques. All the nanoparticles were synthesized using a ligand exchange reaction from the same batch of OT homoligand NPs. Therefore, they should have the same core size and size distribution. For the ligand shell composition, we have previously shown that the deuteration does not affect the ligand shell composition for the same(1) and similar(2) type of ligand coated NPs. Furthermore, while for the SANS modelling, only 50% penalty was imposed on the volume fraction of the ligands, i.e. the ligand shell ratio was not fixed. The volume ratio of the outcome model is MPA: OT = 1: 2.12. While the ligand ratio from NMR measurements is MPA: OT= 55%: 45%, which corresponds to approximately MPA: OT= 1: 2.18, agreeing well with the resulted 3D model. The good quality of the fits also confirms that the ligand ratios are consistent with NMR data.

Ab initio model fitting of SANS data
As described in a recent report, the construction of 3D models from SANS data was performed using MONSA software (3). Briefly, a spherical search volume with a diameter of 7.2 nm was generated according to the pair distribution function calculated using GNOM(4) program. The search volume composed of closely packed beads with a radius of 2 Å. Each bead could be assigned either to the gold, MPA, OT or solvent phase. To facilitate the calculation, the beads within 2 nm radius are fixed to be gold, while beads between 2.7 nm 3.6 nm radius could be assigned to either the two types of ligands or solvent. The ratio of volume fraction of the two types of ligands, calculated from the NMR and TGA value (OT: MPA= 2.0: 1), was used as a fitting penalty during the simulated annealing process. Starting from random configurations, simulated annealing was conducted in order to minimize the discrepancy between the experimental data and the form factor of the bead model. The discrepancy is defined as: Where, the k is the index of scattering curves; Nk is the number of points in the experimental scattering curves; Iexp(q) and Icalc(q) are the intensities from the experiments and the calculation of the bead model; ck is the scale factor and σ(q) is the experimental errors at each q.
The model fitting process was repeated multiple times always starting from random configurations and the resulting models all present similar ligand shell organization features.

Interpretation of the C-H region of SFG spectra
Since the morphological evolution during the annealing is driven by the gain in conformational entropy of the OT ligands, one would expect for the stripe-like NPs a higher gauche to trans ratio for the ligands compared to the patchy one. To measure such conformational changes, we performed SFG measurements also in the C-H stretch region (2800 -3000 cm -1 ) of the vibrational spectrum of NPs as shown in Figure S5. The ligand MPA was deuterated so that only the conformation of the OT was probed. Using SFG the alkyl chain conformations can be determined empirically using the amplitude ratio (referred to as d + /r + ratio) of the symmetric methylene (~2850 cm -1 , d + ) and the symmetric methyl (~2880 cm -1 , r + ) stretch vibrational modes. A value of d + /r + « 1 is associated with a stretched all-trans alkyl chain conformation, whereas a value of d + /r + > 1 indicates that gauche defects dominate the measured vibrational spectrum(5-7). The SFG spectra of both NPs display a different relative amplitude of the d + and r + mode indicating that the average ligand conformation on the NPs is different. The stripelike nanoparticle displays a larger d + /r + indicating less stretched alkyl chains.

Work of adhesion measurement from AFM
Images were collected working in the small-amplitude regime in order to map the work (8). The formula proposed in literature to estimate the work of adhesion from the imaging parameters is the following: where A is the working amplitude, A0 is the free amplitude, is the phase shift, is the molecular diameter of the liquid and is the exponential decay length of the density of work of adhesion at the solid-liquid interfaces.
The method does not provide an absolute measurement, but only a relative one that needs to be calibrated through the measurement of contact angle on reference samples. In fact, since the energy dissipation is calculated with respect to the oscillation of the cantilever in the bulk liquid, the intercept of the calibration line is always given by the surface tension of the imaging liquid, therefore only one reference sample is needed for the calibration. We used a film of MPA protected gold nanoparticles as a calibration system for the measurement of the WSL for the patchy NPs. The value of the work of adhesion calculated for these particles was then used as a calibration system for the stripe-like NPs.
To get a reliable comparison of different samples (and hence a good calibration), the exact same A0 was kept while scanning different samples (differences ≤ 0.05 nm where considered negligible), meanwhile A was kept as constant as possible. Furthermore, images were compared only if collected in the same day and with the same tip.
Each experiment was reproduced twice. A minimum of five nanoparticles per sample were analyzed. Analysis was performed on single particles, extracted from images similar to the ones shown in Figure S7. Data analysis was performed by calculating the average value (and the standard deviation) of ( , @ , ) over all the pixels of a single particle masked from an image. This value was then used to calculate the work of adhesion (after the calibration of the system was performed. In the calibration we used 72 mN/m for the surface tension of water).
For every experiment the average value of WSL on different particles and its standard deviation was calculated. 5 mN/m were added to each standard deviation in order to take into account the error associated with the calibration ( * = + 5 ⁄ ). The average value for different experiments has been calculated and the final standard deviation has been obtained by propagating the * of the single experiments.

Simulation details
All simulations were carried out using the GROMACS/2018.4 software package (9). Interligand interactions were treated using the OPLS/AA forcefield(10), the SPC/E model(11) was used for water, and the gold surface was described with the GolP model (12). For the patchy and stripe-like surfaces, we started with a box size of ~6 nm x 6 nm in the xy plane, with a height of 30 nm. A gold (111) surface with height 5 nm was placed at the center, and ligands attached to the surface (the positions of the sulfur atom of each ligand were constrained throughout the simulation, as were all gold atoms other than the image-charge sites (12) ). Slabs of water with height 5 nm were placed above and below the ligand surfaces and allowed to condense onto the surface. For the trench surfaces, a similar procedure was followed but with the dimensions of the box 12.1 nm x 2.6 nm in the xy plane. This gives 24 columns of ligands, from which we created trenches containing 0, 1, 2, 4, 8, 12, 16, 20, 22, 23 and 24 columns of MPA ligands (an example of the simulation setup for the trenches is shown in figure S9).
Surfaces with random motifs were simulated in boxes with the same dimension as the patchy and striped surfaces, with 20 different sets of random ligands generated, making 40 surfaces.
24 surfaces with tailored motifs were also designed to add to the training set if required. Figure   3A shows an example of the simulation setup for a trench surface after allowing the water slabs to condense.
For every surface, four independent runs were carried out with the initial orientation of ligands about the z axis generated randomly. Simulations were run in the NVT ensemble at 300 K with a timestep of 1 fs and a velocity-rescaling thermostat used to control the temperature (13). After allowing the water slabs to condense onto the ligand surface, the system was equilibrated for 10 ns and statistics were collected over 50 ns. The water-vacuum boundaries produced by expanding the simulation box meant that the pressure was held constant at 0 bar. Long-ranged electrostatics were dealt with using the particle-mesh Ewald (PME) method (14).
The structure of the solvent around the nanoparticle surface was quantified using two properties: the density ρ, obtained by taking a histogram of coordinates of the oxygen atoms in water molecules, and the dipole orientation density µ = <cosθ>ρ, where <cosθ> is the average orientation of the water dipole moment with the surface normal pointing away from the gold surface. The surface excesses were obtained by subtracting the predicted histogram from the one observed from simulations, where the prediction was obtained by taking the profiles calculated using surfaces with either full-MPA or full-OT coverage and linearly combining them: where fOT(z) is the profile from the 0-MPA trench and fMPA(z) is the profile from the 24-MPA trench. xOT and xMPA are the fraction of OT and of MPA ligands respectively.

Results on trench configurations
Based on the expectation that the boundary between OT and MPA regions is responsible for the non-additive behavior observed, we ran simulations, as described in the previous section, on a series of trench geometries, in which the length of the boundary between the two regions is kept constant, but the area of each region differs. If the boundary region alone is responsible for non-additivity, we would expect that the excess properties of interest are independent of the size of the trench, so long as it is large enough that the two boundaries do not "interact" with each other. Figure S11 shows that, indeed, while the surface excesses differ between small and large trench sizes, for large enough trenches the excesses are relatively independent of trench size. This is further borne out by Figure 3A in the main text, in which the local water density is shown in a cut across the trench, and is similarly independent of trench size for large trenches.
When investigating the local density of water, we found that the structure of water in contact with MPA ligands follows one of two regimes, where the two structures are related (in the case of a bare surface) by a 60 o rotation, as shown in Figure S10. To obtain the results in Figure 3A, we separated each surface into one or the other type of surface-water structuring, in order to avoid spuriously averaging over two different types of structure and washing-out the local water density.
These trench configurations give a qualitative idea of how the water structure depends on the size of regions of ligand, but for a quantitative picture a wider variety of surface morphologies must be considered.

Nearest-Neighbor Model for Surface Excesses
In order to predict the surface excesses on arbitrary ligands, we developed a model in which an arbitrary excess quantity Δf(z) (either Δρ(z) or Δµ(z)) is given by the linear combination,  (15,16), with surface excess functions used in our case rather than, e.g., energetics.
We chose to classify ligands based on their nearest neighbors, with the remainder of their surroundings providing a mean-field background. This approach was motivated by the fact that simulations containing trenches of ligands showed the structure of water to be quite independent of the trench size above the smallest trenches, so that long-ranged effects on the local structure were minimal. The nearest-neighbor environment of a ligand is described by the six ligands surrounding it, meaning that there are 2 7 =128 possible ligand plus environment combinations. Of these, many are related by rotational symmetry, and only 26 must be considered, of which the two cases where all six surrounding ligands are identical to the central one were assumed not to contribute to the excess. The remaining 24 possible motifs are shown in Figure S12, along with the fitting functions Δρi(z) and Δµi(z) corresponding to each.
with < a regularization parameter used to avoid overfitting. Figure S12 shows that the contribution of each motif to a given surface excess is essentially of the same order of magnitude: no particular motif is a priori more important than the others in determining the excess. However, the striped surfaces contain both a larger number of environments that give an excess and a larger variety of environments (with several motifs appearing almost exclusively on striped surfaces). These factors lead to more pronounced surface excesses for the striped than for the patchy surfaces.     on different days, we choose to report our data separately. We calculated the peak to peak ratio error of the sum frequency scattering setup on 10 spectra of different samples of sodium dodecylsulfanate stabilized hexadecane nanoemulsions. and the deviation from spectrum to spectrum is less than 10 %. This value is used as an error for the values reported on the above graph. We have never observed an increase in the ratio during an experiment. The fact that during one experiment (C) we did not observe any difference between the two sets of nanoparticles can be explained by the fact that annealing may not have worked perfectly on that occasion. Since in the formation of the film it is likely that similar nanoparticles will  Peak ratio assemble together, it is possible that in this case the laser was focused on an area composed of nanoparticles that were not affected by annealing.