Coarse-grained molecular models of the surface of hair

We present a coarse-grained molecular model of the surface of human hair, which consists of a lipid monolayer, in the MARTINI framework. Using molecular dynamics simulations, we identify a lipid grafting distance that yields a monolayer thickness consistent with atomistic simulations and experimental measurements of hair surfaces. Coarse-grained models for fully-functionalised, partially damaged, and fully damaged hair surfaces are created by randomly replacing neutral thioesters with anionic sulfonate groups. This mimics the progressive removal of fatty acids from the hair surface by bleaching. We study the structure of the lipid monolayer at different degrees of damage using molecular dynamics simulations in vacuum as well as in polar (water) and non-polar ( n -hexadecane) solvents. We also compare the wetting behaviour of water and n -hexadecane on the hair surfaces through contact angle measurements conducted using molecular dynamics simulations and experiments. Our model captures the experimentally-observed transition of the hair surface from hydrophobic (and oleophilic) to hydrophilic (and oleophobic) as the level of bleaching damage increases. By using surfaces with different damage ratios, we obtain contact angles from the simulations that are in good agreement with experiments for both solvents on virgin and bleached human hairs. In both the molecular dynamics simulations and further experiments using biomimetic surfaces, the cosine of the water contact angle increases linearly with the sulfonate group surface coverage. We expect that the proposed systems will be useful for future molecular dynamics simulations of the adsorption and tribological behaviour of hair surfaces.


Introduction
A detailed understanding of the chemical composition and structure of the surface of human hair is important for the development of more effective and sustainable hair care products. 1The human hair consists of a keratinaeous core, the cortex, and a medulla in thicker hair. 2 In virgin hair, the core is covered by thin (∼ 0.5 µm) cuticle layers, which are arranged in an overlapping structure from the root to the tip of the hair. 3The cuticles themselves are constituent of another set of characteristic layers, which are described in detail elsewhere. 2,4 fatty acid monolayer (F-layer) covers the cuticles of virgin hair and governs their surface properties. 5,618methyleicosanoic acid (18-MEA) has been identified as the dominant lipid species within the F-layer. 7,818-MEA is covalently bonded to the protein layers below, 9 mostly via thioester bonds to cysteine groups (Cys-18-MEA). 10owever, chemical or mechanical damage can lead to partial or complete removal of the F-layer. 1,8Even a single bleach treatment can cause significant damage to the F-layer, mostly through oxidation of the thioester groups in Cys-18-MEA to cysteic acid. 11There are spatial variations in the F-layer within individual hairs.Experiments have shown that the amount of 18-MEA decreases from the root of the hair to the tip. 12There are also higher 18-MEA concentrations in the cuticle centre than at its edges, which are more more prone to mechanical damage. 13Moreover, hair from older individuals contains lower 18-MEA levels than younger ones. 14he outer layers of hair are of particular interest since adsorption of formulation components, as well as lubrication between hair fibres, will largely depend on its surface properties. 1 Experimental surface characterization of virgin and damaged human hair has included the measurement of wetting by means of liquid contact angles, [15][16][17][18][19][20] surface charge density [21][22][23] as well as surface topography and friction using the atomic force microscope (AFM) 8,20,24,25 and high-load nanotribometers. 26These studies have shown that compared to virgin hair, damaged hair usually displays increased surface charge density, hydrophilicity, and friction.These effects can be attributed to removal of the 18-MEA monolayer and formation of cysteic acid (sulfonate) groups on the surface. 27,28n addition to experiments, a range of molecular modelling techniques have been used to study the mechanical and surface properties of hair.Akkermans and Warren 29 developed a multiscale modelling approach to the mechanics of human hair fibre.They used molecular dynamics (MD) simulations with different levels of coarse-graining to study the mechanics of keratin dimers and filaments. 29All-atom (AA) models of keratin dimers have also been developed and used to study the effect of disulfide cross-linker on mechanical properties. 30This atomistic model was used to parameterise a coarse-grained model with a bottom-up approach. 31his enabled the mechanical properties of keratin fibrils to be modelled at larger scales, approaching the size of entire hairs. 31cMullen and Kelty 32 used MD simulations with a united-atom (UA) force field to study the conformational and dynamical properties of 18-MEA and eicosanoic acid (EA) monolayers at the air-water interface.Natarajan and Robbins 33 used MD simulations to estimate the thickness of the 18-MEA layer on an ultra-high-sulfur protein surface.Their results suggested that the monolayer thickness was ∼1 nm, 33 which agreed well with previous X-ray photoelectron spectroscopy (XPS) experiments. 34Cheong and et al. 35 used MD simulations with an AA force field to investigate the optimum separation distance between 18-MEA chains grafted on a graphene surface, which was found to be in the range of d graft = 0.49 − 0.65 nm.These results were consistent with the epicuticle model of Negri et al., 36 but their optimal separation distance was much smaller than the frequently cited value of 0.94 nm. 37At the optimal molecular separation distance, the monolayer thickness was between 2.0-2.6 nm, which was consistent with experimental measurements using transmission electron microscopy (TEM), 38 but much thicker than previous results from MD simulations 33 and XPS experiments. 34n addition to structural and mechanical investigations, density functional theory calculations have been used to study the adsorption of explosive compounds onto 18-MEA monolayers. 391][42] A recent study 43 used coarse-grained MD simulations to investigate the adsorption of polymers onto heterogeneous surfaces obtained from the mapping of experimental AFM data of damaged hair surfaces.They investigated the polymer adsorption as a function of chain length and concentration.They also probed the stability of the adsorbed structures under shear flow. 43ere, we present new coarse-grained molecular models for the F-layer found on the outer layer of the hair cuticle in the fully-functionalised state and at several degrees of damage.The coarse-grained MARTINI force field 44,45 is used throughout.For fully-functionalised hair, we construct a three-dimensional surface model of fatty acids grafted to a graphene substrate in a hexagonal arrangement.This model is validated against previous AA-MD simulations.To represent damaged hair, we randomly replace an increasing proportion (0-100 %) of the fatty acid molecules with anionic sulfonate groups.We study how the surface structure changes with increased damage in vacuum, aqueous, and n-hexadecane environments.Using contact angle measurements, we confirm that the MARTINI model can capture the wetting behaviour of fully-functionalised and damaged hair in both polar (water) and non-polar (n-hexadecane) solvents.We anticipate that the proposed coarse-grained molecules will facilitate MD simulations of the adsorption of surfactants and polymers found in hair care formulations onto the surface of hair.The significant reduction in computational cost compared to atomistic force fields will enable the simulation of large macromolecules that are used to improve the feel and appearance of hair. 1 The  46 The EA molecule is represented by C 1 and N a beads. 47The sulfonate groups are represented by negatively charged Q a beads. 48b) Coarse-grained 4:1 mapping of the polarisable water beads 49 , n-hexadecane (C 1 ) and sodium cations with first hydration shell (Q d ). 45c) Nominal grafting distances relative to CG graphene.1][52] Ghost bead positions overlap with graphene beads on the graphene lattice where positions coincide.d) Different degrees of damage (0-100 % fatty acids removed) are considered.
coarse-grained models proposed here will also enable non-equilibrium MD simulations of hair friction in dry and aqueous environments.

Coarse-Grained Hair Model
Coarse-grained MD simulations are used to study both fully-functionalised and damaged hair surfaces.The MARTINI force field (version 2.0) 44,45 was selected for this purpose.The MARTINI model generally represents groups of four atoms into one single coarse-grained bead (4:1 mapping), although smaller bead sizes have also been developed. 53,54This model has been used extensively to study lipids in aqueous environments. 45Of particular relevance to the current study, the MARTINI model has been successfully employed in MD simulations of supported lipid monolayers 55 and bilayers. 47,56chematics of the model structures for fully-functionalised and damaged hair are depicted in Fig. 1.The fully-functionalised hair surface model includes a complete monolayer of 18-MEA groups, which are covalently bound to cysteine amino acids grafted to a rigid graphene substrate.This approach is similar to that used in the AA-MD study by Cheong et al. 35 The 18-MEA molecules consist of C 20 chains with a methyl branch on the eighteenth C atom and are modeled here as a chain of five C 1 (alkyl) beads.In a subset of MD simulations, a chain of four C 1 beads and a terminal C 2 bead is used to represent stronger alkyl-water interactions, 57 due to conformational disorder 58 caused by the methyl branch-containing alkyl group 59 (see Fig. 1 in the ESI † ).The C 1 beads are connected to a C 5 and P 5 bead, which are representative of the thioester and amine/carboxyl groups in cysteine, respectively.This structure is an extension of the palmitoyl cysteine (C 16 chain) MARTINI parametrisation suggested by Atsmon-Raz and Tieleman 46 by one additional C 1 /C 2 bead.The P 5 beads, which represent the outer protein layer of the epicuticle, are grafted to the graphene layer below.The graphene layer prevents molecules from escaping through the monolayer surface and can be used to control pressure and sliding in future non-equilibrium MD simulations.The graphene sheets also uses a 3:1 coarse-grained mapping [50][51][52] Unlike a 4:1 mapping, this retains a hexagonal structure, with a lattice parameter of factor twice that of atomistic graphene.The grafting positions of the EA, Cys-18-MEA, and cysteic acid groups are independent from the graphene bead positions to allow for intermediate grafting distances, similar to those used in the AA-MD simulations by Cheong et al. 35 .These intermediate values would otherwise not be possible using the 3:1 mapping for the graphene.Ghost beads were introduced in the plane of the graphene sheet and solely interact with the P 5 bead of the lipid chain through a bonded potential.For the comparisons with the AA-MD simulations by Cheong et al. 35 , lipid grafting distances of 0.49, 0.65, and 0.85 nm were considered, which correspond to surface coverages of 4.2, 2.4, and 1.4 molecules nm −2 , respectively.From these comparisons, the intermediate grafting distance of 0.65 nm was selected for the contact angle simulations.
In the case of damaged hair, the 18-MEA molecules are randomly removed by breaking the thioester bond (C 1 to C 5 ), while the neutral P 5 bead is replaced by a negatively charged Q a bead, which represents a sulfonate group. 48This is representative of the oxidation of Cys-18-MEA to cysteic acid, which has been observed in several experimental studies of bleached hair. 8,11,28Charge neutrality is maintained by adding the appropriate number of non-polarisable sodium cation (Na + ) beads Q d , 44 which also represent the first hydration shell in the MARTINI model. 45Since they are expected to be easily removed through washing, free fatty acid molecules are not considered in the current MD simulations.We randomly remove chains at a number corresponding to the desired damage ratio defined as: where N rem is the number of fatty acid chains to be removed and N FF is the number of grafted chains in the fully-functionalised system.Discrete damage ratios of χ N = 0 (fully-functionalised), χ N = 0.25, 0.5, 0.75, 0.85, 0.92 (partially damaged) and χ N = 1 (fully damaged) are considered here.Previous experiments have suggested that the surface of bleached hair mostly consists of negatively charged sulfonate groups, with a mean surface coverage of approximately 2.2 molecules nm −2 (2.8 molecules nm −2 in the damaged regions). 28For the chosen grafting distance (0.65 nm), this coverage of sulfonate groups is approximately reproduced by using a damage ratio of 0.85.
The projected surface area ratio between lipid coverage and exposed charged groups: is used as an additional measure for quantifying monolayer damage.Macroscopic wetting theories such as the Cassie-Baxter model 60 rely on the exposed surface area.Here, S MEA and S total are the projected area occupied by the Cys-18-MEA chains and the total projected surface area respectively.The Cys-18-MEA surface area S MEA is obtained from projecting all remaining Cys-18-MEA beads with an assumed radius of r b = 0.235 nm onto the xy plane of the periodic simulation box.This approach accounts for the effects of collapsing Cys-18-MEA chains onto the charged layers below.

Coarse-Grained Liquid-Phase Model
The original MARTINI water model was single-site and non-polarisable. 44Previous CG-MD simulations of supported lipid bilayers have suggested that the polarisable MARTINI water model gives more accurate results. 56n this study, we employ the original polarisable MARTINI water model developed by Yesylevskyy et al. 49 .Other polarisable water models compatible with the MARTINI framework, such as the refPol, 61 BMW-MARTINI, 62 and refined BMW-MARTINI 63 models were also tested in slab geometry simulations of the water-vapor and n-hexadecane-water interfaces. 64The underestimation of the water-vapour surface tension is a well known limitation of MARTINI water models. 45The refPol model gives significantly lower water-vapor and water-hexadecane surface tensions than both the original polarisable model and particularly experiments.It was therefore discarded after preliminary analysis.BMW-MARTINI was originally optimised to reproduce a range of properties, including the water-vapor surface tension.However, both investigated versions suffer from freezing artefacts when moving from a simple slab geometry of water confined between vacuum to a water droplet deposited on a surface in a vacuum environment.Moreover, previous studies 57,63 found that the water-alkane surface tension was significantly over-estimated with the original BMW-MARTINI model.This is an unfavourable side-effect of stronger LJ interactions used to increase the surface tension of water.The BMW-MARTINI model was therefore deemed unsuitable for the current study because the systems contain both water and alkane components.We selected the original polarisable MARTINI water model due to Yesylevskyy et al. 49 for the remainder of the CG-MD simulations.Although the water-vapor surface tension is underestimated by more than 50 % with this model 49 (similar to the non-polarisable MARTINI model 44 ), the n-hexadecane-water and n-hexadecane-water surface tensions are within 15 % and 30 % of the experimental values, respectively. 63,65Further details on the comparisons between the water models can be found in the ESI † .

MD Simulation Details
The molecular dynamics simulations are performed in LAMMPS 66 using the velocity-Verlet 67 integration scheme with a timestep of 5 fs.Non-bonded interactions between beads are represented by the Lennard-Jones (LJ) potential for uncharged beads and the LJ and Coulomb potentials for charged beads.LJ interactions between beads are shifted to zero between the cut-off radius r LJ,cut = 0.9 and r LJ,shift = 1.2 nm. 44Bonds and angles are treated using weak harmonic potentials as in the original MARTINI framework. 44Water unit bond lengths are constrained using the SHAKE 68 algorithm.The full set of force field parameters used in this study can be found in the ESI † .Long-range electrostatic interactions are treated using a slab implementation 69 of the particleparticle particle-mesh (PPPM) method 70 with a Fourier grid spacing of at least ∆x = 0.2 nm, a Coulombic switching radius of r C,cut = 0.12 nm and a relative energy tolerance of 10 −5 .
The systems are periodic in the directions parallel to the surface (x and y).Reflective boundary conditions are imposed in the direction perpendicular to the surface (z) to prevent water molecules from escaping the simulation cell.Interactions of vapor phase beads with these boundaries are rare due to sufficient separation distance of surfaces and droplets from the surface.A nominal system box size for water droplet wetting of L x = 35.7336nm, L y = 30.9462nm and L z = 25 nm was found sufficient to prevent periodic interactions between droplets.For n-hexadecane, which has lower contact angles at moderate damage levels, the system size in x and y was increased to L x = 71.4672nm, L y = 61.8924nm, respectively, to avoid interaction effects between the periodic images of the droplets.
The thickness of a human hair is typically 50-100 µm and each cuticle cell is approximately 5-10 µm long. 27This is much larger than can be directly simulated in the MARTINI framework.The surfaces used in the current CG-MD simulations are therefore representative of a small region of the epicuticle surface.For the damaged surfaces, the outer F-layer is partially removed, which occurs during a single bleach treatment. 11ore severe damage, such as that caused by repeated bleaching, 71 can lead to cuticle removal.This can change the microscale topography 72 and porosity 73 of the hair surface, which cannot be accounted for in the current CG-MD simulation framework.
The initial conditions for the monolayer-vacuum configurations consist of fully extended lipid chains normal to the graphene sheet and counter-ions placed at random positions in a thin sheet (∆z = 0.4 nm) above the monolayer for damaged hair surfaces.Other initial configurations for the counter-ions were tested but led to similar distributions after equilibration.Energy minimization is performed for the monolayer surfaces and counterions in vacuum.The model surfaces are then further equilibrated for 10 ns in the canonical (NVT) ensemble using a Nosé-Hoover thermostat 74,75 set at T = 300 K, prior to extracting any statistics during the last 10 ns of the production run.The graphene layer and corresponding ghost beads are excluded from thermostatting.Convergence to equilibrium values was monitored by means of the global kinetic and potential energies, monolayer thickness and tilt angles as well as first and second moments of the monolayer bead positions.
Coarse-grained droplet wetting simulations were run in the NVT ensemble for at least 20 ns.Statistics were only extracted over the final 10 ns.The initial hemispherical droplet configuration is obtained from bulk simulations in the isoenthalpic-isobaric (NPT) ensemble at atmospheric pressure and a temperature of 300 K. 76 The total number of polarisable water units in the base configuration is N = 18, 121, leading to an initial droplet radius, r ≈ 10 nm.For n-hexadecane, N = 4, 036 molecules are used for the initial droplets, which resulted in a similar initial droplet radius.The hemispherical droplets are placed at the centre of the fully equilibrated monolayer surfaces and integrated in the NVT ensemble along with the lipid molecules.
The possible effect of the droplet size on the contact angle was investigated by testing droplets with approximately half (N = 8, 994), two times (N = 36, 464), and four times (N = 73, 268) the number of beads in the initial hemisphere at t 0 . 53All damage ratios where χ N < 1 were considered for this comparison (see Fig. 2 in the ESI † ).A significant increase in contact angle with droplet size is observed between N = 8, 994 and N = 18, 121, indicating that the former droplet size is too small to be representative of macroscale contact angles. 77The contact angle is more sensitive to droplet size on surfaces with greater the damage ratio.The changes in contact angles for larger droplets (N ≥ 18, 121) is sufficiently small (within the statistical uncertainty intervals) to be disregarded here.The change in contact angle with droplet size can be fit using the modified Cassie-Baxter approach for fuzzy interfaces (see Fig. 3 in the ESI † ). 78However, accurate extrapolation to microscale droplets, as used in the accompanying experiments, would require further simulations at larger droplet sizes, which are beyond the scope of this current study.
Contact angles for water and n-hexadecane are obtained based on the procedure described by de Ruijter et al. 79 and Werder et al. 80 and successfully applied to MARTINI water droplets on graphene surfaces by Sergi et al. 53 .Discretization sizes of ∆z = 0.05 nm and ∆z = 0.025 nm in the surface-normal direction are used for high contact angle droplets (χ N ≤ 0.75) and low contact angles (χ N > 0.75) respectively.Circle fits are obtained from averaging over five individual snapshots within a bin of ∆t bin = 0.25 ns each.One LJ cut-off length normal to the hair surface is excluded at both the base and top of the droplet.As noted in de Ruijter et al. 79 , the droplet base is prone to deviations from the spherical shape as non-bonded interactions with the surface are strongest in that region.This is expected particularly for our heterogeneous surfaces in the case of damaged hair.We furthermore exclude the top of the droplet due to the high angular offset with respect to the droplet surface which leads to poor fitting results for the density profiles in the Cartesian coordinates.The vertical position of measurement for the contact angle is determined as the surface-normal position z corresponding to a threshold of 99 % of the cumulative monolayer mass distribution function: The heterogeneous nature of the damaged surfaces and nanoscale droplet sizes mean that pinning is possible in the wetting simulations. 81To account for this variability, two different initial water droplet positions shifted diagonally across the periodic box (∆x = ±138 nm, ∆y = ±90 nm) and two different random monolayer damage seeds with a water droplet at the centre position are considered.This leads to five trials per damage level which are used to obtain wetting statistics.Similar position and random seed investigations were omitted for n-hexadecane because these results show less variability.
New MD simulations with AA force fields were also conducted to validate the structure of the fully-functionalised CG hair model.For this purpose, we used the optimised potentials for liquid simulations (OPLS) 82,83 force field for the 18-MEA monolayers and the extended simple point charge (SPC/E) 84 force field for water, were employed.Two lipid grafting distances were considered (d graft = 0.49 nm, 0.65 nm).The simulation protocol of fully-functionalised monolayers in vacuum follows the approach by Cheong et al. 35 All bonds containing H atoms were constrained using the SHAKE algorithm. 68A water droplet consisting of 14,642 molecules (droplet radius ∼ 5 nm) was introduced on the equilibrated surface and contact angles were obtained following the same procedure as for the CG model. 77The surface tension of the water-vapour and n-decane-vapour interfaces was recently evaluated for several force fields, including SPC/E and OPLS. 85The authors found good agreement with experimental data for both the water-vapour (underestimated by 15 %) and n-decane-vapour (underestimated by 5 %) interfaces, which underlines the suitability of the AA-MD as reference cases for comparison with the corresponding CG-MD configurations.
3 Experimental setup

Hair contact angle
Dynamic contact angle measurements on hair surfaces were conducted using the Wilhelmy force balance method 87 .The experimental procedure closely follows the approach described in the literature for dynamic contact angle measurements of single hairs. 15,17A ∼2 cm segment of single fibre hair is mounted on a Kruss mount (Hamburg, Germany -Part number: SH0801) with double sided sticky tape tabs.A straight 1.5 cm segment is extended from the mount parallel to the axis of the mount to ensure perpendicular hair contact with the solvent.The diameter of the hair was measured optically.A bath of 99+% n-hexadecane or distilled water is raised at 6 mm min −1 until the lower end of the fibre is wetted.The segment is then immersed 0.2 mm to avoid the cut tip.The hair is further immersed at 3 mm min −1 , stopping every 0.2 mm for one second to measure the wetting force on the hair.This is continued until 6 mm length of hair is measured.
To validate the coarse-grained hair models at different levels of damage, we investigate two types of hair samples experimentally: virgin (pristine) and medium bleached (single bleach treatment).Five hairs from 120 virgin and 80 medium bleached hair swatches were analysed at the tip, middle and root positions.The results given below correspond to the mean, minimum, and maximum of these measured values.

Biomimetic surface contact angle
Biomimetic hair surfaces were produced by functionalising atomically-smooth silicon wafers with eicosanoic acid (C 18 ) and/or sulfonate (SO - 3 ) groups.Virgin hair was represented by a pristine C 18 -terminated monolayer, while complete damage was represented by a SO - 3 -terminated monolayers.Intermediate levels of damage (e.g.medium bleached and heavily bleached) are represented by co-functionalised surfaces.The use of biomimetic surfaces eliminates any damage-induced change in microscale surface topography 72 or porosity, 73 from the experimental contact angle measurements.They therefore facilitate a more direct comparison to our MD simulation results.
][91] Both temperature (23.1 ± 0.3 • C) and humidity (50 ± 1%) were controlled.Approximately 3 µL droplets of DI water droplets were deposited on horizontally positioned silica wafers by dangling a drop about 200 µm from the surface and slowly growing the drop at 0.5 µL s −1 until it detached from the needle.Profiles of the droplets wetting and equilibrating on the surface were captured at 50 frames per second for 10 s.The First Ten Angstroms (Portsmouth, U.S.A.) software (Version 2.1, Build 381) was used to extract contact angles from the movie.Equilibrium contact angles were determined and measurements were repeated five times for each sample.

Monolayer Properties
Figure 2 depicts the variation in film thickness and tilt angle of the equilibrated fully-functionalised lipid monolayers at increasing nominal grafting distance.Comparisons to the AA-MD reference data by Cheong et al. 35 are also shown.A simplified CG model equivalent to their AA model of EA has been included to facilitate more direct comparison to our CG-MD results.This simplified CG EA model only differs from our final model (Cys-18-MEA) by means of the P 5 bead being omitted and a change in bead type from C 5 (thioester) to N a (ester) for the first bead in the lipid chain.The monolayer thickness is evaluated as the mean of the distance in the surface-normal direction z between the P 5 bead (N a for EA) and the bead with the highest position in z.The length of a fully extended EA and Cys-18-MEA chain in the MART INI framework is ∼3 nm.Note that the definition for the tilt angle used by Cheong et al. 35 is based on the difference in centre of mass between the top and bottom C atoms in the lipid chain within a finite set of regions.Such a definition effectively captures regional anisotropic lipid tilt, which is not what is of relevance in our CG-MD simulations.Instead, the tilt angle in our CG-MD simulations is computed as an average over all individual lipid chain tilt angles between the P 5 (N a for EA) and the terminal C 1 beads.
All of the molecules and models in Figure 2 show a general collapse of fatty acids with increasing grafting distance, as indicated by the decreasing monolayer thickness and increasing tilt angles.As the grafting distance increases, the molecules move from being aligned mostly perpendicular to the surface (tilt angle = 0 • ) to aligning more parallel to the surface.This is due to a reduction in inter-chain van der Waals interactions as the distance between them increases and a concurrent increase in van der Waals interactions between the chains and the graphene layer as more of the surface is exposed. 92There is excellent agreement between the AA-MD 35 and CG-MD results for the thickness and tilt angle of the EA monolayer at the three chosen grafting distances.In general, the CG monolayers are thicker than the AA ones, although the difference is always less than the van der Waals radius of a single CG bead (∼0.5 nm).The thickness of the Cys-18-MEA decreases faster than EA due to the presence of the P 5 bead, is more sensitive to packing distance and enables the chains to align more parallel with the surface to maximise the van der Waals interactions.The estimate of the optimum separation distance of the 18-MEA molecules of d graft = 0.49 nm to d graft = 0.65 nm presented by Cheong et al. 35 from their AA-MD simulations is supported by our CG-MD results.This grafting distance yields a thickness that is consistent with previous TEM experiments of virgin hair. 38Moreover, the cysteic acid molecular density presented by Korte et al. 28 for bleached hair closely agrees with a grafting distance of 0.65 nm.In the remainder of this study, we only consider a nominal d graft = 0.65 nm for both fully-functionalised and damaged hair.
We also investigated the effects of random damage to the lipid monolayer.Several methods for quantifying the degree of damage on the hair surface can be applied within our CG-MD simulation framework.The simplest approach is to directly use the damage ratios by means of the number of randomly removed fatty acids, χ N , and surface coverage, χ S , introduced in Eq. 1 and Eq. 2 respectively.The area damage ratio, χ S , is selected as the damage quantification parameter in the remainder of this study, since it is most representative of the degree of damage that can be measured experimentally (e.g. from AFM scans 43 ).The relationship between the various damage ratios in vacuum and the presence of wetting fluids is shown in Fig. 4. The morphology of the damaged hair surfaces will depend on the external influences causing such damage.Random removal of lipid chains is likely from chemical treatments such as bleaching.Other authors have observed roughly circular 43 and striped 28 damage patterns on hair surfaces.However, the size of these regions is comparable to the size of the surfaces employed in our CG-MD simulations.
Fig. 3 shows the variation in the film thickness and tilt angle of the monolayers at a nominal grafting distance of 0.65 nm with an increasing degree of damage (χ S ) in vacuum and in the presence of water and nhexadecane.In vacuum, a collapse of the lipid chains is observed with an increasing degree of damage, similar to that observed with increasing grafting distance.
Introducing wetting fluids to the surface leads to intermolecular interactions with the solvent beads as well as between the chains.The addition of n-hexadecane molecules on the surfaces does not have any significant effect on the thickness and tilt angles compared to those obtained in vacuum.This is because the van der Waals interactions between the n-hexadecane molecules and the grafted chains are similar in strength to those between proximal grafted chains.On the other hand, swelling is observed for the damaged surfaces as water penetrates into the monolayer.In polymer brushes, such swelling behavior is observed when the interactions between solvent molecules (i.e.water hydrogen bonding) are strong compared to the interactions between solvent and polymer (i.e.van der Waals). 93This leads to a a slight increase of the monolayer thickness with higher degrees of partial damage in water.The degree of swelling increases with increasing damage because a reduced number of lipid chains result in weaker inter-chain interactions and allow increased water penetration into the monolayer.This is consistent with experimental measurements that have shown an increase in swelling in less densely packed monolayers. 94The maximum degree of swelling relative to the vacuum case (∼0.5 nm) is consistent with previous experimental measurements of ionic surfactant monolayers formed on mica surfaces from water. 94Swelling is further confirmed by increased area damage ratios with water as opposed to the vacuum cases (see Fig. 4).
Figure 4 shows the evolution of the surface coverage of the hair monolayer (χ S ) for different degrees of damage (χ N ) in vacuum, water, and n-hexadecane environments.For all of the partially damaged systems, χ S is lower than χ N .This is because the flexible lipid molecules rearrange to cover the surface, increase the surface-molecule van der Waals interactions, and reduce the surface energy.The surface coverage decreases less for water than for vacuum and n-hexadecane environments.This is due to the swelling behaviour observed in Fig. 3

in water.
Surface coverage maps corresponding to the surface damage ratios in Fig. 4 are shown in Fig. 2 in the ESI † for vacuum, water and n-hexadecane cases.The coverage maps were obtained by projecting individual beads using their van der Waals radius onto a continuous planar surface.At 25 % substitution, separated charged damage islands are observed within an otherwise dense monolayer surface.When half of the lipid molecules are removed, the size of the damaged islands increases and they begin to merge.As a result of lipid reorganisation to maximise the chain-surface van der Waals interactions, less than half of the underlying charged surface is exposed.The surfaces at higher damage ratios then show a transition from charged damaged islands to isolated lipid islands.These islands of agglomerated lipid chains shrink in size as the damage ratio is increased further.At 92 % substitution, some islands constituent of single lipids are observed due to the mean distance to neighboring lipids being too large to be bridged even when fully collapsed onto the charged surface below.The shape of the damage patterns is qualitatively similar to that observed in previous AFM experiments. 23,28,43imilar coverage surface coverage projections are also observed for water and n-hexadecane (Fig. 2 in the ESI) † .For water, more of the charged surface is exposed at the same damage ratio due to swelling of the monolayer.Water penetrates into the lipid layer and reduces agglomeration of the lipid chains.For nhexadecane, the surface coverage is essentially identical to the vacuum case.This is because of the similar interaction strength between the grafted chains and the solvent in this case, which means that swelling is negligible (Fig. 3).As a result of these similar interactions, there is some coarsening of the coverage pattern underneath the n-hexadecane droplet, but this does not affect χ S .

Liquid Droplet Wetting
The CG hair surface model is further evaluated by placing nanoscopic water and n-hexadecane droplets on the surface and measuring their contact angle.Examples of the fits to the density profiles of water from the CG-MD simulations are shown for fully-functionalised and partially damaged surfaces in Fig. 5. Examples of the contact angle fits for CG-MD n-hexadecane and AA-MD water are shown in Fig. 4 in the ESI † .For the fully-functionalised system shown in Fig. 5 a), there is negligible penetration of water into the closely-packed lipid layer and the lipid-water interface is easily defined.The water penetration distance into the fully-functionalised monolayer is consistent between the AA-MD (1.29 nm) and CG-MD (0.85 nm) simulations (see Table 1 in the ESI † ).For the damaged systems shown in Fig. 5 b) and c), water molecules penetrate into the lipid monolayer and form layers in the z-direction (see Fig. 5 in the ESI † ).The penetration distance and solvent layering increase with increasing damage level for water.For n-hexadecane, the opposite trend is observed (see Fig. 6 in the ESI † ); the molecules penetrate into the lipid monolayer for the fully-functionalised system but the solvent and monolayer form separate layers for the damaged systems.The penetration of both fluids into the monolayers layering justifies the exclusion of the region close to the interface from the density fitting for the contact angle determination.There is good agreement between the mass density profiles from the AA-MD and CG-MD simulations for fullyfunctionalised surfaces (see Fig. 7 in the ESI † ).
The nominal contact angle data from the CG-MD model, AA-MD and experimental data of dynamic contact angle measurements for water and n-hexadecane on virgin hair are given in Table 1.Experimental values for advancing and receding contact angles for medium bleached hair are given in Table 2. Water shows much larger hysteresis between the advancing and receding contact angles, which is partially due to the swelling effects described above. 95The advancing contact angles are used for comparison with the MD simulations because the receding contact angles are particularly sensitive to microscale effects (e.g.cuticle edge trapping 17 ), which cannot be captured in the current simulations.The preferential use of the advancing contact angle is common practice. 95Previous experiments have suggested that the pseudo-static water contact angle (as used in the MD simulations) is similar to the advancing contact angle on both virgin hair (both 103 • ) and medium bleached hair (86 • vs. 79 • ). 18Good agreement is obtained between the AA-MD and CG-MD wetting configurations for the fully-functionalised hair surface.In both the experiments and CG-MD simulations, the fullyfunctionalised/virgin hair surfaces are hydrophobic (water contact angle > 90 • ) and oleophilic (n-hexadecane contact angle < 90 • ).The experimental advancing contact angles for virgin hair are slightly overestimated by both of the molecular modeling approaches.This suggests that even the virgin hair surfaces in the experiments are not completely covered by a lipid surface.The observed deviations from experiment are deemed acceptable given the imposed simplifications and assumptions of the models, specifically the underestimated water-vapor surface tension. 44The effect of terminal bead type (C 1 or C 2 ) on water and n-hexadecane contact angle for the different damage levels was investigated in further CG-MD simulations, as shown in Fig. 1 in the ESI † .As expected, the use of the C 2 bead, which has stronger LJ interactions, resulted in a lower contact angle for both water and n-hexadecane. 57Use of the C 2 bead improves the agreement with experimental results for the water contact angle on virgin hair compared to C 1 , but leads to poorer agreement for n-hexadecane (Table 2).Therefore, for the remainder of this study, we employ a C 1 terminal bead.surface fit.A detailed description of the uncertainty quantification process is described in detail in Appendix A. The fully-functionalised hair surface is characterised by large water droplet contact angles, which is indicative of a hydrophobic surface.The droplet shape in this case is well defined because the hair surface is homogeneous and there is clear separation between the water and lipid monolayer layers.This is reflected in relatively small uncertainties in the water contact angle at low damage ratios.Increasing the damage ratio leads to exposure of the charged SO - 3 groups and a reduction of water contact angles, suggesting a more hydrophilic surface.7][18]  moderate random damage, an approximately linear reduction of the water contact angle with damage ratio is observed.Full wetting occurs for water when the lipid F-layer is completely removed.Previous experiments have measured contact angles as low as 48 • for heavily bleached hair, 18 which from Fig. 6 suggests that around 85% of the lipid molecules have been removed.An increase in uncertainties of the contact angles is observed with lower contact angles resulting from both the increase in the droplet radius and increasing heterogeneity of the hair surface.The transition from a strongly hydrophobic to an increasingly hydrophilic character of the surface agrees with the experimental observations by means of water contact angles of virgin and bleached hair, as shown in Table 1 and Table 2, respectively.To match the experimental advancing contact angle of virgin hair (106±3 • ), a damage ratio, χ S ≈ 0.2 is required in the CG-MD simulations (see Fig. 8 in the ESI † ).A damage ratio larger than zero seems reasonable for virgin hair since it is subjected to weathering effects and is therefore unlikely to be completely covered by a pristine 18-MEA monolayer.The experimental advancing contact angle of medium bleached hair (72±11 • ) can be approximately matched by using χ S ≈ 0.8 (see Fig. 8 in the ESI † ).
For n-hexadecane, low contact angles are sustained over a wide range of area damage ratios.Electrostatic interactions between the wetting fluid and the surface are not relevant due to the uncharged, non-polar nature of the coarse-grained n-hexadecane molecules.The hair surface becomes more oleophobic as the damage ratio increases, as indicated by an increase in n-hexadecane contact angles.The experimental measurements of the advancing contact angle of n-hexadecane on virgin (36±7 • ) and medium bleached (72±8 • ) hair confirm this trend.These experimental contact angles can be approximately matched in the CG-MD simulations by using χ S ≈ 0.2 (χ N ≈ 0.25) for virgin hair and χ S ≈ 0.8 (χ N ≈ 0.85) for medium bleached hair, which is similar to the values required to match the experimental water contact angle (see Fig. 8 in the ESI † ).The observed χ N for both water and n-hexadecane are also consistent with chromatography measurements of the content of the 18-MEA on virgin hair (25 % reduced compared to the root content) and single-bleached hair (80 % content removed compared to the root). 11he dynamic contact angle experiments are further complemented by static contact angle measurements on biomimetic silica wafer surfaces with different degrees functionalisation of eicosanoic acid (C 18 ) and sulfonate (SO - 3 ) group coverage.The experiments on these model surfaces are not perturbed by hysteresis, curvature effects and mesoscale cuticle effects.The Cassie-Baxter 60 model is applied to the experimental contact angle data to estimate the surface area coverage of C 18 and SO - 3 groups respectively.Fig. 7 compares the cosines of contact angles of water from the CG-MD simulations with the corresponding experimental values.A good agreement in the slope of the contact angle as a function of damage is found between the experimental and numerical values for water at moderate damage (χ S < 0.7).The offset in absolute water contact angle values is a result of the underestimated water-vapor surface tension of the polarizable MART INI water model 44,49 (and consequent overestimation of the contact angle 53 ).Lines of best fit corresponding to the linear Cassie-Baxter 60 model have been included for the experimental and simulation data.The Cassie-Baxter line of best fit for the simulations lies within the uncertainty bounds for all of the data points.Reasonable agreement between the slopes of the experimental and simulation data fits is obtained suggesting that the Cassie-Baxter model (or its extensions 78 ) can serve as a simple approach to describe the hydrophobichydrophilic transition of the idealised surface of the epicuticle of virgin and bleached hair.

Conclusions
We have developed a coarse-grained molecular model of the lipid monolayer on the surface of the hair epicuticle within the MART INI framework for use in MD simulations.The 3D models presented accurately reproduce the surface properties of virgin and medium bleached human hair.
The surface coverage of lipids grafted to the underlying protein layer is determined from comparisons to atomistic MD simulations and previous TEM experiments.We replicate the effects of bleaching by randomly removing different fractions of fatty acid molecules from the surface and replacing them with anionic sulfonate groups.This leads to nanoscale heterogeneities due to the clustering of the remaining lipid chains to form damaged islands.The size of these islands increases with an increasing degree of damage.
The wetting characteristics of the models show good agreement with the trends observed in experiments and atomistic MD simulations.The highly hydrophobic surface in the fully-functionalised state becomes increasingly hydrophilic as the degree of damage increases, up to complete wetting when all of the fatty acids are removed.Conversely, the surfaces become increasingly oleophobic as the level of damage increases.Accounting for the uncertainties from the two methods, we can define number damage ratios (fatty acid/sulfonate) in the MD simulations that replicate the experimental water and n-hexadecane contact angles for both virgin hair (χ N = 0.25) and medium bleached hair (χ N = 0.85).The cosine of the water contact angle increases linearly with the area damage ratio, as predicted by the Cassie-Baxter equation.The slope of this increase is similar to that obtained with biomimetic hair surfaces functionalised with different alkyl/sulfonate ratios.
These results pave way for further MD simulations of adsorption on the hair surface and non-equilibrium MD simulations of the tribological behavior of hair.These could facilitate the virtual screening of surfactants and polymers for the development of improved hair care formulations.

Fig. 1
Fig.1Schematics of the coarse-grained MARTINI model44,45 of the hair surface.a) Coarse-grained mapping of the Cys-18-MEA, EA and cysteic acid molecules.The fully-functionalised 18-MEA molecules are represented by P 5 , C 5 , C 1 , and C 2 beads.46The EA molecule is represented by C 1 and N a beads.47The sulfonate groups are represented by negatively charged Q a beads. 48b) Coarse-grained 4:1 mapping of the polarisable water beads49 , n-hexadecane (C 1 ) and sodium cations with first hydration shell (Q d ). 45c) Nominal grafting distances relative to CG graphene.The graphene layer, as employed in previous atomistic MD studies of the hair surface,35 is represented by a 3:1 mapping to C 1 beads.[50][51][52]Ghost bead positions overlap with graphene beads on the graphene lattice where positions coincide.d) Different degrees of damage (0-100 % fatty acids removed) are considered.

Fig. 2
Fig. 2 Snapshots of Cys-18-MEA lipid monolayers at different grafting distances rendered using VMD 86 a).C 1 beads are shown in black, C 5 in blue, and P 5 in green.EA and Cys-18-MEA monolayer b) film thickness and c) tilt angles compared to AA-MD reference data for EA and 18-MEA from Cheong et al. 35 .Lines are guides for the eye.The thickness and tilt angles are measured from the centre of mass of the N a bead to the terminal C 1 bead for EA and the P 5 bead to the terminal C 1 /C 2 bead for Cys-18-MEA.

Fig. 3
Fig. 3 Snapshots of monolayers in vacuum from χ S = 0 (top) to χ S = 1 (bottom) rendered using VMD. 86a).C 1 beads are shown in black, C 5 in blue, and P 5 in green.Monolayer thickness b) and tilt angle c) at a nominal grafting distance of d graft = 0.65 nm and various degrees of random damage in vacuum environment and with water and n-hexadecane.Lines are guides for the eye.

Fig. 4
Fig. 4 Surface coverage ratio, χ S , as a function of the number damage ratio, χ N , for the CG hair surfaces in vacuum, water, and n-hexadecane.Snapshots show the projected lipid coverage for the partially damaged surfaces in vacuum.Details of the calculation of χ S and examples of the projected lipid coverage for water and n-hexadecane are shown in Fig. 2 in the ESI † .

Fig. 5
Fig. 5 Circular surface fits of a water droplets on the a) fully-functionalised, χ N = 0, and partially damaged, b) χ N = 0.5, c) χ N = 0.85, monolayers from the CG-MD simulations.Four temporal bin fits for the droplet surface (solid circles) are shown.Averaged contact angles (solely from the shown four temporal bins) are indicated by bold red lines.The surface-normal positionfor the contact angle measurement is shown (dashed red line).Individual surface density fits in the radial direction are omitted for the sake of clarity.Note that the CG-MD simulation domain is larger than the radial section shown here.

Fig. 6
Fig.6shows the change in the contact angles of water and n-hexadecane droplets on the hair surface as the level of damage is increased.The vertical bars consider contributions from temporal fluctuations, different initial droplet positions on the hair surface, different random damage seeds, and uncertainties of the droplet

Fig. 6
Fig. 6 Change in the water and n-hexadecane contact angles on model hair surfaces with an increasing degree of random damage a).Points show mean values from three different initial droplet positions and three different initial random damage seeds.Vertical bars show the minimum and maximum values observed across the temporal binning of the different position and seed trials.Solid lines are guides for the eye.Snapshots of water (blue) and n-hexadecane (green) droplets at low, medium, and high levels of damage rendered using VMD b). 86Surface coverage ratios χ S are computed independently for monolayers with water and n-hexadecane.The inset shows the change in the cosine of the angle, rather than the angle.

Fig. 7
Fig. 7 Cosine of the water contact angle from the CG-MD simulations at different degrees of damage and experiments using silica wafer surfaces (Exp.) at various degrees of C 18 /SO - 3 functionalisation.Simulation points show mean values from three different initial droplet positions, vertical bars represent minimum and maximum values observed across the temporal binning of the different position and seed trials.Dashed lines are Cassie-Baxter 60 fits (cos(θ ) = χ S cos(θ d ) + (1 − χ S )cos(θ ff )) predict contact angles for fully-functionalised and fully damaged monolayers of θ CG,ff = 121 • (θ exp,ff = 102 • ) and θ CG,d = 14 • (θ exp,d = 5 • ) respectively.Open symbols represent points where the surface coverage has been obtained from the Cassie-Baxter relation based on experimental contact angles.The inset shows the change in the angle, rather than the cosine of the angle.

Table 1
Contact angle data (in degrees) on fully-functionalised hair from CG-MD (MART INI) and reference AA-MD (OPLS-AA, SPC/E) and experimental advancing and receding contact angle data for virgin hair.Ranges shown are maximum and minimum values measured.

Table 2
Experimental advancing and receding contact angle data (in degrees) for medium bleached hair.Ranges shown are maximum and minimum values measured.

Table 3
Contact angle data from CG-MD for the C 1 type terminal lipid beads at various degrees of random lipid damage