Crystalline Wax Esters Regulate the Evaporation Resistance of Tear Film Lipid Layers Associated with Dry Eye Syndrome

Dry eye syndrome (DES), one of the most common ophthalmological diseases, is typically caused by excessive evaporation of tear fluid from the ocular surface. Excessive evaporation is linked to impaired function of the tear film lipid layer (TFLL) that covers the aqueous tear film. The principles of the evaporation resistance of the TFLL have remained unknown, however. We combined atomistic simulations with Brewster angle microscopy and surface potential experiments to explore the organization and evaporation resistance of films composed of wax esters, one of the main components of the TFLL. The results provide evidence that the evaporation resistance of the TFLL is based on crystalline-state layers of wax esters and that the evaporation rate is determined by defects in the TFLL and its coverage on the ocular surface. On the basis of the results, uncovering the nonequilibrium spreading and crystallization of TFLL films has potential to reveal new means of treating DES.

. Representative WE topology (BO). Atoms with L-OPLS nonbonded parameters are colored black. Atoms with OPLS-AA non-bonded parameters are colored red. Bonds with OPLS-AA parameters are colored black, and bonds with L-OPLS torsion potential are blue. Initial configurations were created using the Packmol software. 9 The monolayer simulations were equilibrated for at least 50 ns or until the total energy of the system had converged, which took up to 300 ns in some systems. Analysis was performed over 50 ns following equilibration. The bulk WE simulations were equilibrated for 45-495 ns, until the density was constant, and then analyzed over the production simulations.

Bulk systems composed of wax esters
The employed simulation model was first validated via bulk WE simulations. To this end, initial configurations consisting of 200 WE molecules were generated with approximately half of the molecules in a disordered conformation, and half in a crystalline lattice, based on the experimental crystal structure of saturated wax esters. [10][11][12] Each configuration was simulated until a completely crystalline or disordered state was reached. Final densities as a function of temperature are shown in Figure S2. Bulk simulations were performed for palmityl palmitate, palmityl laurate, and BO. Good agreement was found between the simulated and experimental melting points of palmityl palmitate and palmityl laurate. The simulated melting point of BO was ~10°C higher than the experimentally determined melting point (38°C), 13,14 possibly due to its very long alkoxy chain (22 carbons), which is beyond the parametrization used for L-OPLS hydrocarbons. 2 Therefore, reported temperatures were scaled down by a factor of 1.032 in BO monolayer simulations to allow direct comparison with experiments. Figure S2. Simulated melting points of behenyl oleate (BO), palmityl palmitate and palmityl laurate. Shaded regions depict standard deviation. Dotted lines show experimentally determined melting points. 13,14 Wax ester monolayers For the WE monolayer simulations, we used a symmetric configuration with two monolayers containing 100 or 200 BO molecules on both sides of a 4.7-nm thick water slab. Since crystal nucleation from the liquid phase takes place in time scales beyond those achievable with atomistic MD simulations, 70% of BO molecules were placed in a tightly packed extended conformation, to act as an initial crystal nucleus. The rest of the BO molecules were initially placed in a hairpin conformation, with the ester group facing the water surface. The initial conformation states were based on existing experimental evidence of wax ester monolayer structure. 15,16 Water permeation simulations Film permeability for water was determined using the inhomogeneous solubility-diffusion model. 17,18 In the permeation simulations, umbrella sampling (US) was employed to calculate the potential of mean force (PMF), while the depth-dependent diffusion coefficients were estimated from the mean squared displacement (MSD). The distance from the water surface along the direction normal to the surface was employed as the reaction coordinate for the US simulations. Starting configurations for the permeation simulations were generated at 0.1 nm intervals along this reaction coordinate by pulling a single water molecule through both monolayers. A harmonic potential with a force constant of 1000 kJ mol -1 nm -2 was applied to the sampled water molecule along the z-coordinate. Each simulation window was equilibrated for 1 ns followed by a 4-ns production simulation for the disordered and crystal defect cases, and a 19-ns production simulation for the completely crystalline monolayer. Temperature was maintained at 37°C during the permeation simulations. Weighted histogram analysis method, implemented in the GROMACS tool g_wham, 19 was used to construct the PMF. To ensure the convergence of the PMFs, the obtained hydration free energy of water was compared to the value calculated from free-energy perturbation simulations as the latter is free of errors due to hysteresis. In the free-energy perturbation simulations one water molecule was decoupled from the water phase, effectively placing it in the vacuum phase. The free energies were calculated using the Bennett Acceptance Ratio method. 20 First, Coulombic interactions were decoupled linearly, followed by decoupling of the van der Waals interactions using a soft-core function 21 with parameters α = 0.5, σ = 0.3 and p = 1. The resulting hydration free energy for the OPC water was 31.5 ± 0.2 kJ mol -1 . This is higher than the value determined experimentally (26.4 kJ/mol), 22 but is consistent with other non-polarizable water models, such as SPC/E. 23 The commonly used positional autocorrelation method 24 for determining the diffusion coefficient along the reaction coordinate (normal to the water surface) was found to produce artificially low diffusion coefficients near the water surface ( Figure S3), since the slow motion of the interface contributes to the position and velocity autocorrelation functions. Therefore, the lateral diffusion coefficient, determined using the MSD, was used to approximate the diffusion coefficient along the reaction coordinate. This appears to be acceptable, since the MSD and autocorrelation methods agree in the regions not affected by surface fluctuations.

Methods
Lipids were added to the air-water interface of a KSV Mini trough (Biolin Scientific, Stockholm, Sweden) in chloroform solution. Chloroform was allowed to evaporate for 5 min. The film was compressed continuously at a rate of 3-12 Å 2 molecule -1 min -1 . Slower compression was used at small mean molecular areas. Surface potential was determined using KSV SPOT (Biolin Scientific, Stockholm, Sweden). Since the surface potential probe is sensitive to the distance from the water surface, a linear correction with respect to time was applied to each measurement to consider the lowering of the water surface due to evaporation of water ( Figure S4). Temperature was maintained within ± 1 °C of the target temperature using a Lauda ECO E4 thermostat (Lauda, Germany). A KSV NIMA microBAM (Biolin Scientific, Stockholm, Sweden) was used to record images of the films during compression. Due to the limited trough area, final results were combined from separate compressions ranging from 300 Å 2 molecule -1 to 80 Å 2 molecule -1 and from 90 Å 2 molecule -1 to 18 Å 2 molecule -1 . Figure S4. Surface potential calibration of a pure water surface used to remove the linear drift caused by evaporation of the subphase.

Hole healing in a solid BO monolayer
To demonstrate the healing time scale of a solid BO film, BO was initially spread to a small area of the trough (3 Å 2 /molecule), after which the barriers were rapidly opened to an area of 16 Å 2 /molecule. Temperature was maintained at 37 ± 1 °C, and the film was observed using microBAM. The solid film spread to roughly complete coverage within 1 h, after which evolution of small holes was followed ( Figure S5).    * L-OPLS force field parameters. † HC-CM-CT-CT torsion potential is undefined in the original OPLS-AA force field. However, in contrast to this study, HC-CM-CT-CT was replaced with HC-CM-CT-OH potential in the original L-OPLS force field.

Topology files
For behenyl oleate, palmityl palmitate, and palmityl laurate, the used topology files in a GROMACS-compatible format are included separately in the Supporting Information.