[1]R.J. Braun

1]Department of Mathematical Sciences, University of Delaware, 210 South
College Avenue, Newark, 19716, DE, USA

# On Modeling Tear Breakup Dynamics with a Nematic Lipid Layer

M.J. Taranchuk<EMAIL_ADDRESS><EMAIL_ADDRESS>[

###### Abstract

One of the main roles of the lipid layer (LL) of the tear film (TF) is to help
prevent evaporation of the aqueous layer (AL). The LL thickness, composition,
and structure all contribute to its barrier function. It is believed that the
lipid layer is primarily nonpolar with a layer of polar lipids at the LL/AL
interface. There is evidence that the nonpolar region of the LL may have
liquid crystalline characteristics. We investigate the structure and function
of the LL via a model of the tear film with two layers, using extensional flow
of a nematic liquid crystal for the LL and shear-dominated flow of a Newtonian
AL. Evaporation is taken into account, and is affected by the LL thickness,
internal arrangement of its rod-like molecules, and external conditions. We
conduct a detailed parameter study with a focus on the evaporative resistance
parameter, the Marangoni number, and primary liquid crystal parameters
including the Leslie viscosities and director angle. This new model responds
similarly to previous Newtonian models in some respects; however,
incorporating internal structure via the orientation of the liquid crystal
molecules affects both evaporation and flow. As a result, we see new effects
on TF dynamics and breakup.

###### keywords:

Tear film, Nematic liquid crystal, Lipid layer, Marangoni effect

##  1 Introduction

With each eye blink, a thin liquid tear film is established that helps protect
and nourish the ocular surface, and helps provide a good smooth surface for
clear sight. The tear film (TF) is a multi-layer film, and each layer plays a
role in TF dynamics [1]. Despite the importance of the TF, there is still much
that is not understood. This paper develops a new model for how the TF may
fail, which is often called tear breakup (TBU). The purpose is to better
understand the function and structure of the outermost layer of the TF, called
the lipid layer (LL). A prime motivation for understanding TF dynamics and TBU
is that these phenomena are thought to contribute to the initiation and
progression of dry eye disease (DED), which affects millions of people [2, 3].
For comprehensive views of clinical aspects of DED and the tear film see Refs.
[4, 5, 6, 7, 8]. Reviews related to the tear film include Refs. [9, 1, 10]. In
this introduction, we now describe the overall structure of the LL. We then
introduce previous work on properties, proceed to mathematical models of the
tear film and its dynamics, and finally discuss details of the LL and our
goals for this paper.

###  1.1 Overall structure of the tear film

At the anterior interface with air is an oily LL that retards evaporation [11,
12], and helps to retain a smooth, well-functioning tear film [13]. The LL
also contains polar lipids that are surface active and lower the surface
tension between the LL and the aqueous layer (AL), which is posterior to the
LL [14, 15]. The AL consists mostly of water [16], but it also contains
various ions and proteins [11]. At the corneal surface, which is posterior to
the AL, there is a region with transmembrane mucins protruding from the
epithelial cells of the cornea or conjunctiva. This forest of glycosolated
mucins is called the glycocalyx. It is generally agreed that the presence of
the healthy hydrophilic glycocalyx prevents the tear film from dewetting the
ocular surface [17, 18, 19].

The overall thickness of the TF is a few microns [20], while the average
thickness of the LL is on the order of tens to 100 nanometers when the eye is
open [13, 21, 22, 23, 10]; the thickness of the glycocalyx is a few tenths of
a micron [19, 24]. The overall structure is reformed on the order of a second
after each blink in a healthy tear film [10, 14, 25, 26].

TBU occurs when a thinned region forms in the TF. Clinically, this is defined
as the first dark area that is observed in the fluorescent TF following
instillation of fluorescein dye [27]. Various mechanisms are thought to cause
different types of TBU, or to combine in various ways: evaporation of water
from the AL [28, 7, 29], Marangoni-driven tangential flow (from non-uniform
polar lipid surface concentrations) [30], or dewetting-driven tangential flow
(from uneven ocular surface conditions) [31, 32, 33]. Here, we use the phrase
tangential flow to mean flow along the eye. The evaporation of water from the
AL causes relatively slow thinning [34], while one may expect rapid thinning
to be explained by Marangoni-driven tangential flow [30] or dewetting-driven
tangential flow. Concentration gradients from surface active polar lipids in
the LL induce shear stress at the LL/AL interface and drive outward tangential
flow via the Marangoni effect [14, 15]. Dewetting from irregularities in a
corneal surface region has been hypothesized to generate outward tangential
flow from internal pressure gradients inside the AL due to van der Waals type
forces [35, 36, 31]. A related term is full-thickness tear film breakup (full
thickness-TBU), which is when the AL has thinned to the point where the LL and
glycocalyx touch [29, 37]. It is our belief that the strongest evidence to
date regarding mechanisms for TBU are evaporative and Marangoni-driven
thinning and their interaction [38, 39]; unambiguous evidence of van der
Waals-driven dewetting is not yet in hand in our opinion [29]. We now discuss
studies of these mechanisms in more detail.

The aqueous part of tear fluid is primarily supplied from the lacrimal gland
near the temporal canthus (corner) and the excess is drained through the
puncta near the nasal canthus [40]. Water lost from the tear film due to
evaporation into air is an important process [11, 41, 42]. We believe that
this is the primary mechanism by which the osmolarity (essentially, the total
concentration of ions on a molar basis) is increased in the tear film [10].
Some water is supplied from the ocular epithelia via osmosis when the
osmolarity increases above the isotonic level (the so-called hyperosmolar
state) [43, 44, 10, 41, 45]. According to the Dry Eye WorkShop (DEWS) reports,
the osmolarity is important because it plays an important etiological role in
the development of dry eye disease [28]. An accepted view in the DEWS reports
is that repeated exposure to sufficient hyperosmolarity in tears, as may occur
in TBU, can result in inflammation and ocular surface damage in a vicious
cycle that may lead to DED. There is still considerable ongoing research to
more fully understand the causes and etiology of DED; some mathematical
modeling results discussed in this paper help to quantify the relationship
between tear breakup and hyperosmolarity.

###  1.2 Measurements and properties

Both TF instability and hyperosmolarity are important to study because they
are proposed as etiological causes of dry eye disease [6, 7]. An osmolarity
difference between the corneal epithelium and AL induces osmotic flow from the
cornea to the TF [46, 10]. TF osmolarity may be measured in the inferior
meniscus clinically [47]; the healthy range is 296-302 mOsM [47, 48, 49]. Dry
eye measurements in the same location can reach 316-360 mOsM [50, 48, 51, 52]
but estimates for the TF over the cornea reach 900 mOsM or higher [53, 10, 46,
54]. High levels of TF osmolarity are associated with pain, inflammation and
cellular changes [55, 56, 53]. In support of these potentially high levels of
TF osmolarity over the cornea, mathematical models without spatial variation
have estimated peak osmolarities up to ten times the isotonic concentration
[43, 10]. The modeling work of Peng et al. [46] found that evaporation
elevates osmolarity in breakup regions.

TF thinning rates have been measured experimentally or estimated in many
studies. A few experimental methods include spectral interferometry [57, 42,
34], an open chamber [58], an evaporimeter [59], and averaging pure water and
functioning LL rates over the cornea obtained by heat transfer analysis and
thermal imaging [60]. Parameter estimation schemes were developed in Luke et
al. [54, 39] for fitting partial differential equation (PDE) models to
experimental fluorescent intensity distributions [61]. They found evaporative
thinning rates ranging from 36.9 to $-4.91$ $\mu$m/min (negative indicates
thickening) and overall TF thinning rates ranging from 23.5 to 1.85
$\mu$m/min. These estimates were quite similar to previous experimental spot-
thinning rates that did not specifically seek TBU instances [57].

Tiffany and Dart [62] measured values of the viscosity of human meibum (the
material expressed from the meibomian glands inside the lids) using a
capillary tube and they estimated values from 9.7 to 19.5 Pa s at 308 K. This
is a representative value for the temperature on the surface of the tear film
[63]. In a series of papers by Tiffany and coworkers, the human tear film has
been shown to have non-Newtonian properties [64, 65, 66]. Tiffany [66] found
that whole tears, that is, tear sampled directly from the eye, are shear
thinning. The shear thinning has been fit with empirical models, namely Cross
[66] and Ellis models [67]. Tiffany [64] also found weak elastic effects;
Tiffany and Nagyová [68] found that removing all lipids from tears caused tear
fluid to become Newtonian. The rheology of the LL depends on temperature [69];
in vitro observations showed temperature dependence of meibum on saline for
monolayers [70, 71] and for bulk samples [72]. The elastic contributions to
viscous response was small above 305 K for these measurements [70]. King-Smith
et al. [23] used high-resolution microscopy to observe LL dynamics, and they
saw fine structure of the LL during the interblink period that included
islands of lipids.

The surface tension of the tear/air interface was measured to be 46.2 mN/m by
Miller [73] in vivo using a special contact lens; this is about 2/3 that of an
air/pure-water interface. Pandit et al. [65] measured the surface tension of
whole tears with a capillary tube and obtained $43$ mN/m for unstimulated
tears and $46$ mN/m for stimulated tears; Miller [73] likely measured
stimulated tears. A value of $45$ mN/m is often used in theoretical studies.
The components responsible for lowering the surface tension from that of pure
water has been the subject of a number of papers [74, 68, 75, 76]. A challenge
for this measurement is the dependence of surface tension on surface active
components like phospholipids and other lipids [77, 78, 79] and still other
species [80]. Detailed equations of state for studying the Marangoni effect
from these species is not available for the tear film in vivo to our
knowledge, and this is an ongoing area of research.

###  1.3 Mathematical modeling

A variety of mathematical models have incorporated various important effects
of tear film dynamics as has been recently reviewed [43, 81, 10]. The most
common assumptions for these models are a Newtonian tear fluid and a flat
cornea [14, 82]. One category of models treats the flow over an open eye-
shaped domain; these models have used lubrication theory with two space
dimensions and time as independent variables [83, 84, 85, 86, 87]. These
papers were on a stationary domain, and emphasized the effects of supply and
drainage, osmolarity in various locations in the tear film, and the overall
appearance of fluorescein imaging of the tear film. Some initial studies have
examined flow on a blinking domain [88, 89, 81], but they have yet to
incorporate the same level of effects as in the static boundary papers. There
have also been useful compartment models developed for tear film dynamics [90,
44]. These models use compartments of the tear film as dependent variables in
systems of ordinary differential equations (ODE) and thus avoid the
computational difficulties of solving PDEs on 2D domains. The models can be
solved over many blinks and give insight into the osmolarity and other
variables at various operating conditions.

There are quite a few models that have been developed using a single variable
in space and time for independent variables. There are two main categories
here: (1) a spatial domain that extends across the eye opening vertically,
with or without a moving end, and (2) a shorter domain that is a local model
for TF dynamics. We first discuss category (1). First, a number of papers
considered a domain that extended across the open eye with fixed ends which
typically called post-blink drainage models [91, 92, 93, 94]. These models
typically imposed large menisci at the domain ends that drained fluid from the
relatively flat interior of the TF, and that drove the dynamics. Notable among
these papers is that evaporation was included to find that it could drive
thinning to TBU [94], and that a moving end, representing the upstroke of the
blink, was treated as a coating flow problem [92]. Both of those aspects would
find many applications in later papers. Heat transfer and evaporation were
considered in Li and Braun [95], and later treated together with thermal
imaging in Dursch et al. [60]. A counterpoint to these papers Berger and
Corrsin’s treatment of the Marangoni effect in post blink spreading [14],
which provided a basis for using this effect in many later papers.

Using a moving end on a one-dimensional domain to simulate the upper eyelid
during a blink has been employed in a number of papers. The first to do this
was by Wong et al. [92]; they treated the problem as a coating flow analysis
and developed a closed form equation for the thickness of the deposited film
as a function of end (lid) speed. Later papers extended this analysis by
computing the thickness and other TF variables numerically. These included a
combined TF formation and drainage model [96], and then adding surfactant
transport on the surface of the TF [97]. Additional papers using the eye
opening and interblink together included gravity and reflex tearing [98];
using optimization varying viscosity to simulate and design eyedrops [67]; a
detailed look at insoluble surfactant [99]; a (Newtonian) two-layer model of
the tear film with insoluble surfactant between them [100]; the effect of
(constant) globe curvature on tear film thickness [101]; and a shear-thinning
(power law) model for TF thickness on a curved substrate [102]. Some papers
also treated the moving end for complete blink cycles. Simple sinusoidal
motion was treated in Braun and King-Smith [103], while more realistic motion
was treated in Heryudono et al. [104], and insoluble surfactant was added in
Zubkov et al. [105]. Heat transfer was later considered with either full [106]
or partial [107] blink cycles.

We next discuss category (2) which are local one-dimensional models for TF
dynamics. These models are along a short line segment and are designed to be
local models of TF flow. The effects of evaporation and the Marangoni effect
on TBU have been studied in a number of papers. Evaporation drove all the
dynamics in the models in [46, 10, 61], though the treatments were different.
In [46], a stationary LL with variable thickness caused increased localized
thinning; also included were osmolarity transport in the AL and osmosis across
the AL/cornea interface. An important result from that work is that diffusion
of osmolarity out of the region of evaporation prevents osmosis from stopping
thinning as it would in spatially uniform models [43, 10]. Simpler evaporation
models were used in [10, 61], but they included transport of fluorescein so
that they could explain what was seen in TF visualization experiments. Zhong
et al. [30] developed a PDE model with one spatial variable that incorporated
both mechanisms. Including fluorescein dye transport and fluorescence enabled
fitting of models to in vivo fluorescence data within TBU regions to estimate
parameters that are not possible to directly measure there at the time of
writing [54, 39]. Recently, ODE models with no space dependence have been
successfully fit to fluorescence data in small TBU spots and streaks [108].
Those models have been coupled to a neural-network based data extraction
system to greatly expand the amount of TBU instances that may be studied
[109].

Other effects have been included in local models of TBU [31, 32, 46, 61, 110,
111, 30, 112, 113, 114]. Some models incorporated the effect of soluble mucins
in the AL necessitating departure from more standard lubrication models [112,
113]. Localized non-wettability of the ocular surface was modeled in [114].
Other models have provided fundamental analyses on evaporation [115, 116] and
instability [117, 118].

###  1.4 LL structure and the model in this paper

The composition and structure of the LL, and their relation to function, are
the subject of much ongoing research [77, 119, 120]. The complexity of the
composition affects the determination of structure and function, and has led
to a variety of perspectives in this area [80, 121]. It is generally thought
that the LL has an outermost nonpolar LL with a layer of polar lipids and
other surface active components at the LL/AL interface [122, 1]. Measuring all
but the most basic properties can be a challenge, however. Microscopy has
shown some example forms that the LL may take [23]; the shapes observed there
suggest non-Newtonian behavior of the LL [71, 72]. For a number of years, it
was very difficult to reproduce significant resistance to evaporation using in
vitro meibum or LL films [123, 124], even though resistance to evaporation has
been thought to be an important property in vivo [11, 34]. However, recent
progress on in vitro films with many components has been made [12, 125]. The
barrier function does depend on thickness [34, 126], but clearly there is an
effect from internal structure as well [127, 128].

Molecular dynamics simulations, together with Langmuir trough experiments,
have also given insight into the structure and function of the LL or model
systems for it. Coarse grained models were initially employed [129, 130] and
those models suggest that the LL may not be so as well organized as layers of
nonpolar lipids atop a polar lipid monolayer. Later, all-atom models have
suggested that the specific polar lipids [131, 132] and wax esters [133] may
form ordered states to slow evaporation of water from the AL. Molecular
dynamics simulations are likely to continue to give important insights into LL
structure and function, but they currently can only be solved at small space
and time scales. A link to bigger scales has been made via Langmuir trough
observations with Brewster angle microscopy, e.g. [131, 132] for equilibrium
states. Those insights need to be converted to larger scale dynamic models to
address other situations such as TBU. In this paper, we develop a two-layer
continuum model that attempts to incorporate some internal structure of the LL
that may be expected to produce evaporation resistance [38, 29] that depends
on that structure. We use a nematic liquid crystal structure as a model for
the LL, which can show birefringence as has been observed for meibum [134].

Adapting two-layer models of tear dynamics [100, 46] to dynamic models of
local TBU regions [135, 136] has made use of extensional flow for the LL and
shear dominated flow for the AL. These models have a common basis from the
approach of Matar et al. [137]. Of particular interest here is including
extensional flow of a nematical liquid crystal [138, 139] for the LL with a
Newtonian shear-dominated AL layer. That extensional flow modeling used the
Ericksen-Leslie equations and multiple scale perturbation methods to reduce
the equations to a small PDE system. A detailed derivation and preliminary
results for the combined model were developed in Taranchuk and Braun [136].
With this new model, they found dependence of TF dynamics on the director
angle $\theta_{B}$ with all other parameters fixed, and we continue this work
with a detailed parametric study. We find new ways the effects combine to
affect TF dynamics and TBU.

The paper is organized as follows. In Section 2 we describe the problem
formulation. In Section 3 we present our results. These include profiles of
the film thicknesses, velocities and pressures that lead to different outcomes
depending on the internal orientation of the liquid crystal molecules.
Finally, in Section 4 we discuss the results and outline our conclusions.

##  2 Model

Our two-layer model of the tear film lies on top of the corneal surface, which
we take to be flat. It is composed of a Newtonian AL containing various
solutes and a liquid crystal LL containing non polar lipids. At the interface
of these two fluid layers is a monolayer of surfactant comprising primarily
polar lipids. Water in the AL is lost through evaporation, $J_{e}$, and
supplied from the corneal surface through osmosis, $J_{o}$. Figure 1 shows a
sketch of the model. In the AL, we solve for the thickness $h_{1}$ and the
osmolarity concentration $c$. In the LL we solve for the following quantities:
thickness $h_{2}$, axial velocity $u_{2}$, angle of the director $\theta$, and
surface concentration of surfactant $\Gamma$.

Figure 1: Schematic of the two-layer tear film model.

We model the AL using the Navier-Stokes equations, and the LL using the
Ericksen-Leslie equations for nematic liquid crystals. We use the following
scalings, where primes denote dimensional quantities. The coordinates
$(x^{\prime},y^{\prime})$ and velocity components
$(u_{i}^{\prime},v_{i}^{\prime})$ correspond to the axial and transverse
directions respectively, $t^{\prime}$ is time, and $p^{\prime}_{i}$ is
pressure. The subscript $i$ indicates the fluid layer under consideration;
$i=1$ corresponds to the AL, and $i=2$ corresponds to the LL.

$\displaystyle x^{\prime}$ $\displaystyle=Lx,\qquad\qquad
h_{1}^{\prime}=\epsilon Lh_{1},\qquad\qquad h_{2}^{\prime}=\epsilon\delta
Lh_{2},$ $\displaystyle y^{\prime}$ $\displaystyle=\epsilon Ly,\hskip
36.0ptu_{1}^{\prime}=Uu_{1},\hskip 44.0ptu_{2}^{\prime}=Uu_{2},$
$\displaystyle h^{\prime}$ $\displaystyle=\epsilon Lh,\hskip
36.0ptv_{1}^{\prime}=\epsilon Uv_{1},\hskip 42.0ptv_{2}^{\prime}=\epsilon
Uv_{2},$ $\displaystyle t^{\prime}$ $\displaystyle=\frac{L}{U}t,\hskip
38.0ptp_{1}^{\prime}=\frac{\mu_{1}U}{\epsilon^{2}L}p_{1},\hskip
33.0ptp_{2}^{\prime}=\frac{\mu_{2}U}{L}p_{2},$ (1)
$\displaystyle\Gamma^{\prime}$ $\displaystyle=\Gamma_{0}\Gamma,\hskip
38.0ptc^{\prime}=Cc,\hskip 50.0pt$ $\displaystyle J_{o}^{\prime}$
$\displaystyle=\epsilon UJ_{o},\hskip
30.0ptJ_{e}^{\prime}=\epsilon\rho_{1}UJ_{e}.$

Here the total thickness of the tear film
$h^{\prime}=h^{\prime}_{1}+h^{\prime}_{2}\implies h=h_{1}+\delta h_{2}$. The
scalings used in the model, chosen with the dimension of the human eye in
mind, are defined in Table 1. The LL is the thinnest layer in the tear film,
and has an average thickness of tens to 100 nm during the interblink [13, 140,
141, 142]. The AL is much thicker, averaging about 3 $\mu$m [143]. The
glycocalyx is less than a micron thick [19, 24]; while it does not appear in
the model explicitly, we account for its height by defining TBU as an AL
thickness of 0.5 $\mu$m or less. We consider a domain of length 2 mm.

Table 1: Parameters used in dependent variable scalings Parameter | Value | Units | Description | Reference
---|---|---|---|---
$L=H_{1}(\gamma_{1}+\gamma_{2})/(\mu_{1}E_{0})^{1/4}$ | $3.1779\times 10^{-4}$ | m | Length scale | [135]
$H_{1}$ | $3.5\times 10^{-6}$ | m | AL thickness | [143]
$H_{2}$ | $5\times 10^{-8}$ | m | LL thickness | [140]
$V=E_{0}$ | $5.093\times 10^{-7}$ | m/s | Thinning rate | [135]
$U=V/\epsilon$ | $4.624\times 10^{-5}$ | m/s | Horizontal velocity | [135]
$C$ | $300$ | mOsM | Reference osmolarity | [46]
$\Gamma_{0}$ | $4\times 10^{-7}$ | mol/$\text{m}^{2}$ | Surfactant concentration | [99, 100]
|  |  |  |

The nondimensional quantities that result from scalings are given in Table 2,
and the physical parameters are listed in Table 3. We note that while
experiments show evidence of liquid crystal structure in the LL of the TF, the
parameter values for the LL have not been measured. For liquid crystals, we
return to a Newtonian fluid by setting the Leslie viscosities
$\alpha^{\prime}_{i}=0,\,i=1,\cdots,6,\,i\neq 4$, and
$\alpha^{\prime}_{4}=2\mu_{2}$. The dynamic viscosity of the LL, $\mu_{2}$,
has previously been taken to be 0.1 Pa s [100, 135]. Thus, to match this
viscosity in the Newtonian limit, as default values we choose our Leslie
viscosities to be three times the properties of 5CB, a liquid crystal that has
been relatively well studied near physiological temperatures [144]. We explore
a larger range of values in Section 3.1.3.

Table 2: Nondimensional Parameters Parameter | Formula | Value | Description
---|---|---|---
$\epsilon$ | $H_{1}/L$ | $0.0110$ | Aspect ratio of AL
$\delta$ | $H_{2}/H_{1}$ | $0.0143$ | Lipid to aqueous thickness ratio
$\Upsilon$ | $\epsilon^{2}\mu_{2}/\mu_{1}$ | $0.0091$ | Reduced viscosity ratio
$\mathcal{C}_{1}$ | $\epsilon^{3}(\gamma_{1}+\gamma_{2})/(\mu_{1}U)$ | $1.0001$ | Reduced aqueous capillary number
$\mathcal{C}_{2}$ | $\epsilon\gamma_{2}/(\mu_{2}U)$ | $43.8371$ | Reduced lipid capillary number
$\mathcal{M}$ | $\epsilon RT_{0}\Gamma_{0}/(\mu_{1}U)$ | $187.7687$ | Reduced Marangoni number
$\mathcal{E}$ | $E_{0}/(\epsilon\rho_{1}U)$ | $1.0001$ | Evaporation parameter
$\mathcal{R}_{0}$ | $k_{m}H_{2}/(Dk)$ | $17.68$ | Evaporative resistance parameter
$\mathcal{P}$ | $CP_{c}/(\epsilon U)$ | $0.1355$ | Osmosis parameter
$\textrm{Pe}_{1}$ | $UL/D_{1}$ | $9.1842$ | Péclet number for osmolarity diffusion
$\textrm{Pe}_{s}$ | $UL/D_{s}$ | $0.4898$ | Péclet number for surfactant diffusion
|  |  |
Table 3: Physical Parameters Constant | Value | Units | Description | Reference
---|---|---|---|---
$\rho_{1}$ | 1000 | kg/$\text{m}^{3}$ | Density of liquid water | [145]
$\rho_{2}$ | 900 | kg/$\text{m}^{3}$ | Density of lipid | [100]
$\mu_{1}$ | $1.3\times 10^{-3}$ | Pa s | Dynamic viscosity of water | [100, 66]
$\mu_{2}=\alpha_{4}^{\prime}/2$ | $0.0978$ | Pa s | Dynamic viscosity of lipid | [144]
$\gamma_{1}$ | $0.027$ | N/m | Surface tension of water | [146]
$\gamma_{2}$ | $0.018$ | N/m | Surface tension of lipid | [147]
$R$ | $8.3145$ | J/(K mol) | Ideal gas constant | [145]
$D_{1}$ | $1.6\times 10^{-9}$ |  $\text{m}^{2}$/s | Osmolarity diffusion constant | [148]
$D_{s}$ | $3\times 10^{-8}$ |  $\text{m}^{2}$/s | Surfactant diffusion constant | [100]
$k_{m}$ | $0.0182$ | $\text{m}/\text{s}$ | Mass transfer coefficient | [135]
$Dk$ | $5.136\times 10^{-11}$ | $\text{m}^{2}/\text{s}$ | Lipid permeability | [135]
$P_{c}$ | $2.3\times 10^{-10}$ | $\text{kg}/(\text{mOsM}\ \text{m}^{2}\ \text{s})$ | Corneal osmosis coefficient | [46]
$E_{0}$ | $5.0928\times 10^{-7}$ | $\text{m}/\text{s}$ | Evaporative thinning rate | [135]
$T_{0}$ | 308.15 | K | Eye surface temperature | [46]
$\alpha_{1}^{\prime}$ | $-0.0180$ | Pa s | Leslie viscosity | [144]
$\alpha_{2}^{\prime}$ | $-0.2436$ | Pa s | Leslie viscosity | [144]
$\alpha_{3}^{\prime}$ | $-0.0108$ | Pa s | Leslie viscosity | [144]
$\alpha_{4}^{\prime}$ | $0.1956$ | Pa s | Leslie viscosity | [144]
$\alpha_{5}^{\prime}$ | $0.1920$ | Pa s | Leslie viscosity | [144]
$\alpha_{6}^{\prime}$ | $-0.0624$ | Pa s | Leslie viscosity | [144]
$\theta_{B}$ |  $0$ \- $\pi/2$ | radians | Director angle | [144]
|  |  |  |

###  2.1 Governing equations

After nondimensionalizing, we asymptotically expand all of the dependent
variables in powers of the small parameter $\epsilon$. We then collect like
powers of $\epsilon$, and solve to find a closed system of leading order
equations. In the process, we find that the director angle is constant at
leading order, $\theta=\theta_{B}$. We also determine the leading order AL
pressure $p_{1}$, and the AL depth averaged axial velocity $\bar{u}_{1}$. This
process is described in more detail in the appendix A, and derived thoroughly
in Taranchuk and Braun [136]. We find the closed leading order system to be

$\displaystyle h_{10t}+(\bar{u}_{10}h_{10})_{x}=J_{o}-J_{e},$ (2)
$\displaystyle h_{20t}+\left(u_{20}h_{20}\right)_{x}=0,$ (3)
$\displaystyle-\frac{A(\theta_{B})}{B(\theta_{B})}\delta\Upsilon(u_{20x}h_{20})_{x}-\frac{u_{20}}{h_{10}}=\frac{h_{10}}{2}\bigg{(}-\mathcal{C}_{1}h_{10xxx}-\mathcal{C}_{2}\Upsilon
h_{0xxx}\bigg{)}$ $\displaystyle-\mathcal{C}_{2}\delta\Upsilon
h_{20}h_{0xxx}+\mathcal{M}\Gamma_{0x},$ (4)
$\displaystyle\Gamma_{0t}+(\Gamma_{0}u_{20})_{x}=\frac{1}{Pe_{s}}\Gamma_{0xx},$
(5)
$\displaystyle(c_{0}h_{10})_{t}=\frac{1}{Pe_{1}}(c_{0x}h_{10})_{x}-\left(c_{0}\bar{u}_{10}h_{10}\right)_{x},$
(6)

where

$\displaystyle h_{0}=h_{10}+\delta h_{20},$ (7) $\displaystyle
J_{o}=\mathcal{P}\left(c-\frac{c_{0}}{C}\right),$ (8) $\displaystyle
J_{e}=\frac{\mathcal{E}}{1+\mathcal{R}(\theta_{B})h_{20}},$ (9)
$\displaystyle\mathcal{R}(\theta_{B})=\mathcal{R}_{0}(0.1+0.9\sin\theta_{B}),$
(10) $\displaystyle
p_{10}=-\Upsilon\mathcal{C}_{2}h_{0xx}-\mathcal{C}_{1}h_{10xx},$ (11)
$\displaystyle\bar{u}_{10}=-p_{10x}\frac{h_{10}^{2}}{12}+\frac{u_{20}}{2},$
(12)

and

$\displaystyle A(\theta_{B})=$
$\displaystyle\;2\cos(2\theta_{B})(\alpha_{1}+\alpha_{5}+\alpha_{6}+4)(\alpha_{2}+\alpha_{3}-\alpha_{5}+\alpha_{6})$
$\displaystyle+\cos(4\theta_{B})(\alpha_{1}[\alpha_{2}-\alpha_{3}]+[\alpha_{2}+\alpha_{3}][\alpha_{5}-\alpha_{6}])+\alpha_{1}\alpha_{2}-\alpha_{1}\alpha_{3}-2\alpha_{1}\alpha_{5}$
$\displaystyle-2\alpha_{1}\alpha_{6}-8\alpha_{1}+\alpha_{2}\alpha_{5}+3\alpha_{2}\alpha_{6}+8\alpha_{2}-3\alpha_{3}\alpha_{5}-\alpha_{3}\alpha_{6}-8\alpha_{3}-2\alpha_{5}^{2}$
$\displaystyle-4\alpha_{5}\alpha_{6}-16\alpha_{5}-2\alpha_{6}^{2}-16\alpha_{6}-32,$
(13) $\displaystyle B(\theta_{B})=$
$\displaystyle-\alpha_{1}\cos(4\theta_{B})+\alpha_{1}-2\cos(2\theta_{B})(\alpha_{2}+\alpha_{3}-\alpha_{5}+\alpha_{6})$
$\displaystyle-2\alpha_{2}+2\alpha_{3}+2\alpha_{5}+2\alpha_{6}+8.$ (14)

Equation (2) and (3) represent conservation of mass (of solvent) in the AL and
LL, respectively. Equation (4) is the axial force balance in the LL. Equation
(5) conserves mass of surfactant, and Equation (6) conserves osmolarity in the
AL.

The evaporation function in Equation (9) depends on the evaporative resistance
in the LL which is described by Equation (10). We vary the evaporative
resistance function $\mathcal{R}(\theta_{B})$ by changing both
$\mathcal{R}_{0}$ and $\theta_{B}$. Various director angles are considered
throughout the paper, and Sections 3.1.2 and 3.2.2 focus on the effect of
$\mathcal{R}_{0}$ for two sets of initial conditions.

###  2.2 Initial and boundary conditions

####  2.2.1 Evaporation-driven flow

To illustrate flow driven by evaporation, we choose initial conditions in such
away that valley, or thin spot, in the LL drives all the dynamics of the tear
film toward TBU. We choose a Gaussian perturbation of the LL

$\displaystyle h_{2}(x,0)=1-0.9e^{-x^{2}/(2\sigma^{2})},$ (15)

where $2\sigma^{2}=1/9$. We take $h_{1}=\Gamma=c=1$ initially. We compute
$p_{1}(x,0)$ from Equation (21), and solve the discrete version of the axial
force balance Equation (22) using the other initial values to obtain
$u_{2}(x,y,0)$. This initial condition corresponds to the Gaussian disturbance
described in Stapf et al. [135]; we refer to it as a valley here.

####  2.2.2 Marangoni effect

To illustrate surfactant gradient driven flow, we choose initial conditions in
which all variables are uniform initially ($h_{1}=h_{2}=c=1$), except for a
local excess of surfactant in the center of the domain. We compute
$p_{1}(x,0)=0$ and $u_{2}(x,y,0)$ as described above. The initial condition
for $\Gamma$ is

$\displaystyle\Gamma(x,0)=1+s-s\left(\frac{1+0.5\tanh\left(\frac{x-x_{0}}{x_{w}}\right)+0.5\tanh\left(\frac{-x-x_{0}}{x_{w}}\right)}{n}\right)$
(16)

where the height of the region of excess surfactant is $s=1.8$, the inflection
point is $x_{0}=1$, and the width of the region is $x_{w}=0.36$, and
normalization factor $n=0.9915$.

####  2.2.3 Combining LL and surfactant perturbations

Having considered perturbations in the LL and surfactant surface concentration
individually, we now consider two scenarios which combine perturbations in
both. Namely, we consider a bump in $h_{2}$ combined with a local excess of
$\Gamma$, and a valley in $h_{2}$ combined with a local excess of $\Gamma$.

Imaging of the LL has shown instances of regions of thicker lipid [22]. It is
reasonable to assume that when the LL is thicker, there is more surfactant. To
illustrate this, we consider a thin LL with a central bump and central
increased $\Gamma$ in that area. We take $h_{1}(x,0)=c(x,0)=1$, and we
calculate $p_{1}(x,0)$ and $u_{2}(x,0)$. The initial condition for surfactant
is Equation (16), and for LL thickness

$\displaystyle
h_{2}(x,0)=1+s-s\left(\frac{1+0.5\tanh\left(\frac{x-x_{0}}{x_{w}}\right)+0.5\tanh\left(\frac{-x-x_{0}}{x_{w}}\right)}{n}\right)$
(17)

where the height of the bump is $s=1$, the inflection point is $x_{0}=1$, and
the width of the bump is $x_{w}=0.36$, and normalization factor $n=0.9915$.
This corresponds dimensionally to a LL with background thickness of 25 nm, and
a bump of lipid of height 50 nm in the center of the domain.

We also consider a pathological scenario in which a valley in the LL is
combined with a local excess of $\Gamma$. Keeping all other initial conditions
as described above, we modify only Equation (17) by changing $s=-0.5$,
creating a valley in the LL. This corresponds dimensionally to a LL with
thickness 50 nm, and a valley of depth 25 nm in the center of the domain.

####  2.2.4 Boundary conditions

At both boundaries we impose the following no flux boundary conditions,

$\displaystyle h_{1x}=h_{2x}=p_{1x}=\Gamma_{x}=c_{x}=0,$ (18)

and specify $u_{2}=0$ for the LL velocity.

###  2.3 Numerical solution

Before solving numerically, we make two change of variables. In order reduce
the order of the system, we define $p_{10}$ as a dependent variable (Equation
(21)) and substitute into Equation (4). We also use the substitution
$m_{0}=h_{10}c_{0}$ in Equation (6), where $m_{0}$ is the mass. This yields
the following closed system which we wish to solve

$\displaystyle h_{10t}+(\bar{u}_{10}h_{10})_{x}=J_{o}-J_{e},$ (19)
$\displaystyle h_{20t}+\left(u_{20}h_{20}\right)_{x}=0,$ (20) $\displaystyle
p_{10}=-\Upsilon\mathcal{C}_{2}h_{0xx}-\mathcal{C}_{1}h_{10xx},$ (21)
$\displaystyle-\frac{A(\theta_{B})}{B(\theta_{B})}\delta\Upsilon(u_{20x}h_{20})_{x}-\frac{u_{20}}{h_{10}}=\frac{h_{10}}{2}p_{10x}-\mathcal{C}_{2}\delta\Upsilon
h_{20}h_{0xxx}+\mathcal{M}\Gamma_{0x},$ (22)
$\displaystyle\Gamma_{0t}+(\Gamma_{0}u_{20})_{x}=\frac{1}{\mbox{Pe}_{s}}\Gamma_{0xx},$
(23) $\displaystyle
m_{0t}=\frac{1}{\mbox{Pe}_{1}}\left(m_{0x}-\frac{m_{0}h_{10x}}{h_{10}}\right)_{x}-\left(m_{0}\bar{u}_{10}\right)_{x},$
(24)

where

$\displaystyle h_{0}=h_{10}+\delta h_{20},$ (25) $\displaystyle
J_{0}=\mathcal{P}\left(\frac{m_{0}}{h_{10}}-1\right),$ (26) $\displaystyle
J_{e}=\frac{\mathcal{E}}{1+\mathcal{R}(\theta_{B})h_{20}},$ (27)
$\displaystyle\bar{u}_{10}=-p_{10x}\frac{h_{10}^{2}}{12}+\frac{u_{20}}{2}.$
(28)

We then use the method of lines to solve the system numerically subject to the
initial and boundary conditions of Section 2.2. Spatial derivatives are
approximated with second-order finite differences on a uniform grid of $512$
points. This results in a system of differential algebraic equations that we
solve forward in time in Matlab (MathWorks, Natick, MA, USA) using ode15s. We
found that the solutions were converged using this number of grid points, and
that the mass of water and osmolarity were conserved to the tolerances of the
computation.

We integrate the system until either of the following happens: (i) a final
time for 60 s is reached, or (ii) TBU occurs. We define TBU as the point at
which the AL thickness reaches a minimum value of $0.5\mu$m, and use event
detection to stop the simulation if this happens. When breakup occurs, the
TBUT becomes the final time.

##  3 Results

###  3.1 Evaporation-driven thinning

Figure 2: Solutions for spatially-uniform initial conditions except for a
Gaussian disturbance in the LL thickness. The director angle is
$\theta_{B}=0$. This figure illustrates typical dynamics of the dependent
variables in space and time for the case of TBU driven by increased
evaporation from the valley in the LL.

As an example of evaporation-driven thinning, we consider the initial
conditions described in Section 2.2.1 with a director angle of $\theta_{B}=0$.
This choice maximizes the evaporation rate for our chosen evaporation
function. Evaporative resistance in the LL is maximized for
$\mathcal{R}(\pi/2)$; $\mathcal{R}(0)$ reduces that resistance by a factor of
10.

Solutions are shown in Figure 2. The LL changes very little over time,
maintaining its original shape with a valley in the center of the domain. The
initial velocity in the LL, which ranges between $\pm 50\,\mu$m, is moderate
compared to the dynamics of other scenarios considered later. The thinning in
the LL leads to increased evaporation of the AL, which quickly thins, reaching
breakup conditions in the center of the domain by 15.9 s. AL pressure begins
to drop in the center of the domain; in general, this would help to restore
fluid to the area; however TBU occurs because that capillarity-driven
contribution is too weak.

Osmolarity $c$ increases over time, peaking in the center of the domain as a
result of increased evaporation in that area. By the final time, osmolarity in
the center of the domain is 1651 mOsM, more than five times the initial value;
and near the ends it is 789 mOsM, over two and a half times the the initial
value.

The surface surfactant concentration $\Gamma$, shown relative to its initial
value, builds up slightly in the center of the domain. We explore Marangoni-
driven flow in detail in Section 3.2; however in this example, the increase in
surfactant concentration is not large enough to contribute significantly to
the flow.

These solutions are representative of evaporation-driven thinning that leads
to breakup before the final time of 60 s. With a larger director angle,
evaporation resistance of the LL is increased, and the AL thins more slowly.
When $\theta_{B}\lesssim\pi/17$, TBU is not reached within the simulation time
of 60 s. Because evaporation is slower, the valley that develops in the AL is
less pronounced, and the resulting peak in osmolarity is smaller; see
Taranchuk and Braun [136] for more details.

####  3.1.1 Varying the Marangoni number

We now consider the effect of the Marangoni number on tear film dynamics as
the director angle $\theta_{B}$ varies from 0 to $\pi/2$. Figure 3 shows the
extreme values of AL thickness, osmolarity, and LL thickness for three values
of the Marangoni number: $\mathcal{M}$ (the default value of 187.7),
$\mathcal{M}/10$, and $10\mathcal{M}$. Solutions are also shown for the
initial conditions of Section 2.2.2, the local excess of surfactant; these are
described in further detail below in Section 3.2.1.

We see that when $\mathcal{M}\rightarrow\mathcal{M}/10$, the Marangoni effect
is weakened, and there are no cases of TBU within 60 s. The LL is the most
mobile in this case, as healing flow (driven by capillarity) allows the valley
in the LL to thicken from its initial minimum of 5 nm, to local minima of 14.7
nm to 33.3 nm, depending on the angle of the director. Osmolarity decreases
monotonically from 1323 mOsM to 467 mOsM.

For the default Marangoni number, TBU occurs before 60 s for
$\theta_{B}\lesssim\pi/17$. Around this critical director angle the LL has
some mobility, with the lipid valley thickening to a maximum of 17.2 nm;
however, for $\theta_{B}\rightarrow\pi/2$, the LL moves very little from its
initial profile.

When $\mathcal{M}\rightarrow 10\mathcal{M}$, TBU occurs before 60 s for
$\theta_{B}\lesssim\pi/9$. The LL remains tangentially immobile for all angles
of the director. Osmolarity decreases over the range of the $\theta_{B}$, but
its decline is slower than for the default Marangoni number.

Figure 3: Extreme values of AL thickness, osmolarity, and LL thickness at the
final time as $\theta_{B}$ is varied from $0$ to $\pi/2$, for three values of
the Marangoni number $\mathcal{M}$. To the left of the black dots, TBUT occurs
before 60 s; to the right, the simulations end at 60 s. The first row
corresponds to the initial condition with a valley in the LL described in
Section 2.2.1. The second row corresponds to the initial condition with a
local excess of surfactant described in Section 2.2.2.

####  3.1.2 Varying the evaporative resistance parameter

Figure 4 shows the effect of varying $\mathcal{R}_{0}$, the evaporative
resistance parameter of the LL. The default value $\mathcal{R}_{0}=17.68$ is
based on a nonlinear least squares fit to lipid thickness and permeability
measurements performed by Stapf et al. [135]; this is shown with the vertical
dashed line. With this level of resistance, only a director angle of
$\theta_{B}=0$ results in a tear breakup time (TBUT) before 60 s. However, if
$\mathcal{R}_{0}$ is decreased by a factor of 10, then TBU occurs in under 15
s for all values of the director angle. As the evaporative resistance
parameter is increased, the minimum AL thickness tends toward the initial
value of 3.5 $\mu$m.

Figure 4: Tear film breakup time or the minimum thickness of the AL after 60 s
as a function of the evaporation resistance parameter $\mathcal{R}_{0}$ for
the initial condition with a valley in the LL. The black dashed vertical line
denotes the default value of $\mathcal{R}_{0}=17.68$.

####  3.1.3 The effect of liquid crystal viscosities on thinning

Figure 5: Minimum thickness of the AL after 60 s or TBUT as a function of the
scale factor applied to the 5CB liquid crystal parameters. To the left of the
black dots, TBUT occurs before 60 s; to the right, the simulations end at 60
s.

We have chosen to set the Leslie viscosities to three times the values of 5CB,
a relatively well studied liquid crystal. By doing this, we match the dynamic
viscosity of the LL, $\mu_{2}$, to previously used value of 0.1 Pa s [100,
135]. However, there is much uncertainty around the values for the Leslie
viscosities as they have not been measured for the LL, and using a scale
factor of 3 is only one possible choice. Thus, in Figure 5, we consider the
effect of a wide range of scalings on the minimum values for both the aqueous
and LL. Note that because we are scaling the dimensional quantities, in the
nondimensional form we are effectively scaling $\mu_{2}$. Here we reduce the
Marangoni number by a factor of 10 ($\mathcal{M}\rightarrow\mathcal{M}/10$) in
order to isolate the effect of the liquid crystal viscosities on tear film
thinning and breakup.

For values of $\mu_{2}$ that are increased by less than a factor of 1000,
there is little change in either the minimum AL thickness or LL thickness at
the final time. For values of $\mu_{2}$ more than a factor of 1000 larger than
default, we begin to see breakup within 60 s for small values of $\theta_{B}$.
However, for larger director angles, the AL is fairly unaffected. On the other
hand, for small values of $\mu_{2}$ we observe a much more mobile LL for all
orientations of the director. As $\mu_{2}$ increases, the LL becomes immobile,
and the minimum thickness remains close to the initial value of 5 nm.

Figure 7 shows solutions at the final time when $\mu_{2}$ is increased by a
factor of 2000. At this value, we see the most variation in LL thickness with
respect to the director angle. For clarity, we omit $\theta_{B}=\pi/18$ and
$\pi/4$. When $\theta_{B}=\pi/2$, the minimum LL thickness is just over 12 nm,
while when $\theta_{B}=0$, the effects of capillarity result in healing flow,
which increases the minimum thickness to 31.7 nm. For the intermediate value
of $\theta_{B}=\pi/12$, the LL experiences some healing flow with a minimum of
24.1 nm, but mobility is reduced.

Figure 6: Solutions for the LL at 60 s when the initial condition has a valley
in the LL and $\mu_{2}\rightarrow 2000\,\mu_{2}$, and
$\mathcal{M}\rightarrow\mathcal{M}/10$.

Figure 7: The coefficient $A/B$ as $\theta_{B}$ varies for scalings of 3, 10,
and 16 on the Leslie viscosities $\alpha_{i},\,i=1,\cdots,6,\,i\neq 4$ for the
liquid crystal 5CB.

It is natural to explore an alternative scenario: keeping a scaling of three
on $\alpha^{\prime}_{4}$, which maintains $\mu_{2}=0.1$ Pa s, while scaling
the other Leslie viscosities, $\alpha^{\prime}_{i},\,i=1,\cdots,6,\,i\neq 4$.
By doing this, we are changing the coefficient $A/B$ on the extensional term
in Equation (4). Note that as $\alpha^{\prime}_{i}\rightarrow
0,\,i=1,\cdots,6,\,i\neq 4$, we return to the Newtonian limit, where $A/B=4$.
However, we find that if we increase the scaling too much (more than a factor
of 16), the coefficient $A/B$ becomes positive for some values of
$\theta_{B}$, making the problem ill-posed. In addition, the effect of the
scalings in this feasible range is negligible for all variables. If we look at
the coefficient $A/B$ for scalings, shown in Figure 7, we see that the change
on the coefficient is very little, even with the maximum scale factor of 16.

###  3.2 Marangoni-driven thinning

We now turn to the initial conditions of Section 2.2.2, in which concentration
gradients drive the tear film dynamics. The Marangoni effect results in rapid
tangential flow; thus we consider solutions on a short timescale, shown in
Figure 8, and a longer timescale, shown in Figure 9. Figure 8 shows the
evolution of the variables over the first 0.06 s, when $\theta_{B}=\pi/2$.
These solutions are representative of the dynamics of tear film thinning due
to surfactant concentration gradients regardless of the orientation of the
liquid crystal molecules; the director angle $\theta_{B}$ has little effect on
the short timescale. The large surfactant gradients immediately cause rapid
thinning in both the lipid and aqueous layers. As the local excess of
surfactant flattens, the Marangoni effect concludes. The LL thins
dramatically, with the valley thickness roughly half of the initial value of
50 nm. Fluid is pulled away from the excess surfactant, and builds up on the
edges of the domain. The pattern of thinning is very similar in the AL,
although to a lesser degree. During this initial redistribution of fluid, the
LL fluid velocity is very large in magnitude (around $3\times
10^{5}\,\mu$m/min), but quickly slows to the range of $\pm 5\,\mu$m/min
throughout the domain (not shown). Similarly, AL pressure (not shown) drops in
magnitude from a range of $\pm 3$ Pa during the first 1 s, to $\pm 1$ Pa, and
thinning continues over a longer timescale.

Figure 8: Solutions on a short time scale for spatially-uniform initial
conditions except for a local excess of surfactant. The director angle is
$\theta_{B}=\pi/2$. This figure illustrates typical dynamics of the dependent
variables in space and time. The surface gradients of surfactant drive all the
dynamics.

Figure 9: Solutions at the final time for spatially-uniform initial conditions
except for a local excess of surfactant. When $\theta_{B}=0$, breakup occurs
at 16.6 s; for all other values of $\theta_{B}$ the final time is 60 s. The LL
profiles (upper right) are nearly identical for all $\theta_{B}$.

Figure 9 shows solutions at the final time for fives values of $\theta_{B}$
ranging from 0 to $\pi/2$. Breakup conditions are reached when $\theta_{B}=0$
by 16.6 s; for all other values of $\theta_{B}$, the final time of 60 s is
reached. The surfactant concentration remains effectively uniform around 1.58
times the original background concentration, and over the longer timescale it
is the evaporative resistance in the LL that influence thinning dynamics. When
$\theta_{B}=\pi/2$ or $\pi/4$ thinning becomes more localized in the center of
the domain, while smaller values of $\theta_{B}$ maintain a wider valley of
thinning through time. For larger values of $\theta_{B}$, we observe that
negative pressure in the center of the AL domain pulls fluid in, and helps to
reduce thinning. This results in an overall thicker AL at the final time.
Arrows in the velocity and pressure plots of Figure 9 point to the difference
in behavior when $\theta_{B}=\pi/4$ or $\pi/2$ as compared to smaller director
angles.

On the other hand, the LL profiles are indistinguishable from each other at
the final time and look very similar to the solution at 0.06 s shown in Figure
8. Virtually all of the movement in the LL happens within the first second;
this immobilization of the fluid surface due to the Marangoni effect has been
studied in other contexts such as vertical drainage in Braun et al. [149] and
remobilization of droplet surfaces in Stebe and Maldarelli [150].

Osmolarity in the AL increases for all values of $\theta_{B}$ over the
original value of 300 mOsM; however, it is most pronounced when
$\theta_{B}=0$. While the solutions in the AL are not dissimilar for
$\theta_{B}=0$ and $\theta_{B}=\pi/18$, the osmolarity is very different. When
$\theta_{B}=0$, peak osmolarity reaches nearly 1400 mOsM by 16.6 s, more than
one and a half times the peak osmolarity at 60 s when $\theta_{B}=\pi/18$.
When $\theta_{B}=0$, thinning and TBU happen so rapidly that the solutes do
not have enough time to diffuse through the AL.

####  3.2.1 Varying the Marangoni number

The second row of Figure 3 shows extreme values for AL thickness, osmolarity,
and LL thickness as the director angle varies from 0 to $\pi/2$ for three
values of the Marangoni number: the default $\mathcal{M}$ (187.7),
$\mathcal{M}/10$, and $10\mathcal{M}$. In general, compared to the initial
condition with a valley in the LL which are discussed in Section 3.1.1,
varying the Marangoni number has a much smaller effect on the variables. This
is because with localized excess surfactant in the initial condition, the
Marangoni effect is already driving dynamics, and our default Marangoni number
is strong enough that increasing or decreasing by a factor of 10 does not have
a dramatic effect.

For all three values of the Marangoni number, TBU occurs for
$\theta_{B}<\pi/50$ when the initial condition has a local excess of
surfactant. For larger director angles, there is no breakup within 60 s, but
the AL thins much more than in the case with the lipid valley, with minimum AL
thickness ranging from 1.62 to 1.84 $\mu$m when $\theta_{B}=\pi/2$. After the
initial rapid redistribution of fluid due to the Marangoni effect, the valley
formed in the LL has a minimum value of 28 nm when $\theta_{B}=0$ for all
three cases. For $\mathcal{M}$ and $10\mathcal{M}$ the minimum remains the
same for all values of the director. For $\mathcal{M}/10$, the minimum
increases slightly as the director angle increases, with a minimum value of
31.6 nm when $\theta_{B}=\pi/2$. Osmolarity between the three cases is very
similar, with a maximum difference of 28 mOsM; in fact, the osmolarity curve
for $10\mathcal{M}$ is indistinguishable from that of $\mathcal{M}$ in Figure
3.

####  3.2.2 Varying the evaporative resistance parameter

We now explore the effect of the evaporative resistance parameter
$\mathcal{R}_{0}$. Figure 10 shows the TBU conditions as $\mathcal{R}_{0}$
varies for various orientations of the director. A small evaporative
resistance parameter results in rapid thinning of the AL with breakup just
under 5 s for all values of $\theta_{B}$. As $\mathcal{R}_{0}$ increases,
breakup time increases, until eventually the final time of 60 s is reached
before breakup occurs. The value of $\mathcal{R}_{0}$ at which this occurs
varies from 2.5 for $\theta_{B}=\pi/2$ to 25.1 for $\theta_{B}=0$. The default
value, denoted by the dashed vertical line, falls into this range where
thinning may result in breakup before 60 s, depending on the director angle.
See Figure 9 for solutions at the final time. As $\mathcal{R}_{0}$ becomes
very large, evaporation is diminished, and the minimum value tends toward 3.1
$\mu$m for all values of $\theta_{B}$.

Figure 10: Tear film breakup time or the minimum thickness of the AL after 60
s as a function of the evaporation resistance parameter $\mathcal{R}_{0}$ for
the initial condition that has a local excess of surfactant. The black dashed
vertical line denotes the default value of $\mathcal{R}_{0}=17.68$.

Figure 11 shows the relationship between the evaporative resistance parameter
and the angle of the director for the two sets of initial conditions. The
value of $\mathcal{R}_{0}$ which results in breakup conditions at 60 s is
shown in red for $\theta_{B}=0,\,\pi/18,\,\pi/12,\,\pi/4,$ and $\pi/2$. The
circles represent the solution from the evaporation-driven flow caused by a
valley in the LL; the crosses represents solutions from flow driven by the
Marangoni effect caused by localized excess surfactant. The black lines are
rational functions fit to the respective data delineating whether or not
breakup conditions are reached by the final time of 60 s. If the evaporative
resistance parameter is decreased by a factor of 10 from the default value of
17.68, then TBU will occur for all angles of the director for both sets of
initial conditions, although TBUT will vary.

If we compare the two curves, we see that the $\mathcal{R}_{0}$ which results
in TBU at 60 s for the case of evaporation-driven thinning is roughly double
that of the case of Marangoni-driven flow. In the evaporative case, when
$\theta_{B}=0$, the evaporative resistance parameter $\mathcal{R}_{0}=50$,
whereas for the Marangoni-driven thinning, $\mathcal{R}_{0}=25$. When
$\theta_{B}=\pi/2$, $\mathcal{R}_{0}$ is 5 and 2.5, respectively. This trend
roughly holds for all values of the director angle. In general, this is
expected behavior; the LL acts as a barrier to evaporation in part due to its
thickness. In the first set of initial conditions, a deficiency in the LL
leads to increased evaporation. Thus, it cannot tolerate an evaporation
resistance parameter value of any less than 50 without TBU occurring before 60
s when $\theta_{B}=0$. On the other hand, in the Marangoni-driven flow, the LL
begins with a uniform thickness of 50 nm which helps to prevent TBU even
though other aspects of the barrier may be compromised.

Figure 11: A comparison of the $\mathcal{R}_{0}$ value which results in TBU at
60 s for two sets of initial conditions: one describing evaporation-driven
flow (red circles), and one describing flow driven by the Marangoni effect
(red crosses), as a function of the director angle of the liquid crystal
molecules in the LL. The black lines are a fit to the respective data. To the
left of each black line, TBU occurs before 60 s; to the right of each black
line, TBU may occur after 60 s.

####  3.2.3 Reducing initial lipid layer thickness

The thickness of the LL ranges between 20 - 120 nm, with more recent
measurements favoring a thinner LL [140, 141, 142]. We halve the initial LL
thickness from 50 nm (a “normal” tear film thickness), to 25 nm (a “thin”
thickness). Figure 12 shows solutions for the AL and osmolarity at the final
time for a few values of $\theta_{B}$. The solutions can be compared with
those of Figure 9 which show solutions for a 50 nm thick initial LL.
Qualitatively, solutions for large values of $\theta_{B}$ are similar,
although the AL is thinner, especially at its minimum. However, for
$\theta_{B}=\pi/12$ or $\pi/18$, there are two minima, located symmetrically
about the center. For $\theta_{B}=0$, the region of thinning is steeper and
more narrow than when the LL is thicker. In all cases, the solutions for
osmolarity show higher concentrations with a thinner LL.

Table 4 shows that a 50% decrease in LL thickness results in a 40.1% decrease
in the minimum AL thickness when $\theta_{B}=\pi/4$; for other values of the
director, the percent decrease is even less. The factor increase in peak
osmolarity over the ideal osmolarity is also shown. These numbers are highest
when TBU occurs; however, with a thinner LL, peak osmolarity is well over
twice the ideal value for all angles of the director.

Table 4: Comparison of minimum AL thickness and maximum osmolarity when the initial LL thickness is reduced from 50 nm to 25 nm. The factor increase in peak osmolarity is compared to the ideal osmolarity of 300 mOsM.  | $h_{1min}$ | $c_{max}$
---|---|---
$\theta_{B}$ | 50 nm | 25 nm | % decrease | 50 nm | factor $\uparrow$ | 25 nm | factor $\uparrow$
| (TBUT) | (TBUT) |  |  |  |  |
0 | 0.5 (16.6 s) | 0.5 (8.2 s) | 0 | 1398 | 4.66 | 1457 | 4.86
$\pi/18$ | 0.6626 | 0.5 (54.8 s) | 24.5 | 922 | 3.07 | 1273 | 4.24
$\pi/12$ | 0.7538 | 0.5528 | 26.7 | 809 | 2.70 | 1134 | 3.78
$\pi/4$ | 1.3361 | 0.8010 | 40.1 | 546 | 1.82 | 768 | 2.56
$\pi/2$ | 1.6472 | 0.9926 | 39.7 | 478 | 1.59 | 659 | 2.20
|  |  |  |  |  |  |

Figure 12: Solutions at the final time for spatially-uniform initial
conditions except for a local excess of surfactant. The initial thickness of
the LL is 25 nm (half the default value). When $\theta_{B}=0$, breakup occurs
at 8.2 s; when $\theta_{B}=\pi/18$, breakup occurs at 54.8 s. For all other
values of $\theta_{B}$, the final time is 60 s.

Figure 13 shows the effect of varying the director angle by plotting
parametric changes to the minimum value of the AL thickness and the maximum
value of osmolarity at the final time as $\theta_{B}$ varies from 0 to $\pi/2$
for both a normal and thin LL. When the LL thickness is 50 nm, TBU occurs
before 60 s only for values of $\theta_{B}<\pi/50$; when the LL thickness is
25 nm, TBU occurs for $\theta_{B}<\pi/14$. The black dots on the plot mark
this change. Minimum LL thickness as $\theta_{B}$ varies is not shown, as
changes are minimal. When the initial LL thickness is 50 nm, minimum values
hovers near 28 nm; when the initial thickness is 25 nm, the minimum values
remain around 14 nm for all values of the director angle.

Figure 13: Extreme values of AL thickness, osmolarity, and LL thickness at the
end of each simulation as $\theta_{B}$ is varied from $0$ to $\pi/2$. To the
left of the black dot, TBUT occurs before 60 s; to the right of this dot, the
simulations end at 60 s.

###  3.3 Combining lipid layer and surfactant perturbations

We now compare the two scenarios which combine perturbations in both the LL
thickness and the surfactant concentration. The first combination is a bump in
$h_{2}$ with a local excess of $\Gamma$, and the second is a valley in $h_{2}$
with a local excess of $\Gamma$. Figure 14 compares solutions for AL and LL
thickness, as well as osmolarity for these two cases with $\theta_{B}=\pi/2$.
While neither case reaches breakup for this angle of the director, thinning in
the AL is more pronounced when there is a valley in the LL, with a minimum
value of 0.94 $\mu$m, as compared to 1.66 $\mu$m when the LL contains a bump.
The LL behavior is similar in both cases; it undergoes rapid redistribution of
fluid initially, and then remains tangentially immobile for the remainder of
the time. When the LL contains a bump, this extra fluid is quickly
redistributed throughout the domain; an area of thinning results, with the
minimum thickness of 28.46 nm. When the LL contains a valley, the Marangoni
effect magnifies this profile resulting in a solution similar to those of
Figure 9, but with a deeper valley: the depth is now 14.15 nm versus 28.1 nm
(when $\theta_{B}=\pi/2$).

Osmolarity in both cases reflects the respective changes in the AL. When the
LL contains a valley, the AL thinning is more pronounced, and thus osmolarity
reaches a peak of 650 mOsM, with little increase in concentration at the edges
of the domain. When the LL contains a bump, the AL thins to a lesser degree,
and the osmolarity steadily rises throughout the domain, with a smaller peak
in the center. Peak osmolarity at the final time reaches 485 mOsM. In all
scenarios that contain a local excess of surfactant initially, the Marangoni
effect results in a rapid flattening of the surfactant concentration to a more
or less uniform value of 1.58.

We can consider the influence of the director angle, by observing the final
time solutions for various angles as shown in Figure 15. When there is a bump
in the LL, we see more variation resulting from changes in the director angle.
When $\theta_{B}=0$, TBU occurs at 46.8 s, but no other angles of the director
result in breakup. When there is a valley in the LL, the director angle has
less of an impact because rapid thinning in the AL leads to TBU for
$\theta_{B}=0,\,\pi/18,$ and $\pi/12$. Breakup times are detailed further in
Table 5 which is discussed more below. Osmolarity at the final time is
strongly influenced by the director angle. When there is a bump in the LL,
peak osmolarity is almost three times higher for $\theta_{B}=0$ (1409 mOsM)
than for $\theta_{B}=\pi/2$ (485 mOsM). The difference is just over double
when there is a valley in the LL (1442 versus 650 mOsM).

Thus far we have considered a valley in the LL with depth 25 nm in a 50 nm LL.
We now vary this profile keeping a valley of half thickness. Table 5 compares
the TBUT for three different LL background thicknesses: 30, 40, and 50 nm. We
see that the dynamics are not sensitive to LL thickness when $\theta_{B}=0$.
This is because evaporative resistance is at its lowest for this angle of the
director; for all of the scenarios we have considered, $\theta_{B}=0$ results
in TBU before 60 s. LL thickness has a stronger effect for intermediate angles
of the director. For example, when $\theta_{B}=\pi/18$, decreasing the LL
thickness by 20%, from 50 nm to 40 nm, results in TBU in half the time. When
$\theta_{B}$ is large enough, TBU does not occur; for all three LL thicknesses
TBU does not occur for $\theta_{B}=\pi/4$ or $\pi/2$.

Figure 14: Left column: Solutions when $\theta_{B}=\pi/2$ for a bump $h_{2}$
combined with a local excess of $\Gamma$. Right column: Solutions when
$\theta_{B}=\pi/2$ for a valley in $h_{2}$ combined with a local excess of
$\Gamma$. In both cases $\Gamma(x,0)$ is given by Equation (16). The different
$h_{2}(x,0)$ appears in the middle row.

Figure 15: Left column: Solutions at the final time for a bump in $h_{2}$ combined with a local excess of $\Gamma$. Right column: Solutions at the final time for a valley in $h_{2}$ combined with a local excess of $\Gamma$. Table 5: TBU time comparison in seconds for various initial LL thicknesses when a valley in the LL is combined with a local excess of surfactant.  | $\theta_{B}$
---|---
LL max/min (nm) | 0 | $\pi/18$ | $\pi/12$ | $\pi/4$ | $\pi/2$
50/25 | 8.2 | 29.6 | 48.4 | - | -
40/20 | 7.4 | 15.0 | 32.2 | - | -
30/15 | 6.7 | 10.7 | 14.3 | - | -
|  |  |  |  |

##  4 Discussion and conclusion

In this paper we explore a new model for the tear film of the eye that was
derived in Taranchuk and Braun [136]. This two-layer model combines a
Newtonian AL with a nematic liquid crystal LL, and links molecular orientation
of the LL with evaporative resistance that depends on the director angle. The
LL acts as a barrier to evaporation of the AL, and thickness alone does not
explain its role in preventing evaporation [127, 128]. By representing the LL
as nematic liquid crystal, we are able to consider possible connections
between the structure and composition of the LL and its affect on tear film
thinning and breakup.

We examined two key drivers of TBU: evaporation and the Marangoni effect. We
consider several scenarios which may lead to TBU. The first is a narrow valley
in the LL, which leads to thinning of the AL due to increased evaporation. The
second is a local excess of surfactant on the AL/LL interface, for which we
see thinning on two timescales. Initially, surfactant gradients lead to rapid
outward flow in both the AL and LL. The surfactant concentration levels out
within the first second, and then evaporation takes over, continuing tear film
thinning over a longer timescale of 60 s. Finally, we considered combinations
of perturbations in both the LL and surface surfactant; one in which there is
a bump in the LL combined with a local excess of surfactant, and one in which
there is a valley in the LL combined with a local excess of surfactant. This
final case depicts a worst case scenario, where both excess surfactant and a
deficiency in the LL drive thinning, and breakup conditions are reached for
$\theta_{B}\lesssim\pi/12$.

The evaporative resistance parameter $\mathcal{R}_{0}$ directly influences how
quickly the AL thins. However, as shown in Figure 11, the minimum value of
$\mathcal{R}_{0}$ that the system can tolerate while withstanding breakup also
depends on the angle of the director, $\theta_{B}$, and the initial condition
for the LL. The evaporation function, given in Equations (9) and (10) and
derived in A.3, shows one way these quantities may be related. While we have
chosen to introduce director angle dependence in the evaporation function [29,
38, 151], we also see that the director angle has an effect on flow dynamics.
Small angles of $\theta_{B}$ cause the LL to become more mobile, which results
from liquid crystal properties in the LL.

In addition, we find that the default Marangoni number has a strong effect on
dynamics, regardless of the initial conditions or other parameter values used.
Thus, in order to observe the effect of the liquid crystal parameters,
especially the Leslie viscosities, we reduce the Marangoni number by a factor
of ten. In the absence of any measurements in the eye, we use the properties
of 5CB as a starting point for the Leslie viscosities of the liquid crystal
lipid layer. One unexpected finding is that varying the Leslie viscosities
from the 5CB values has little effect on the dynamics. Not only that, but we
are constrained in how much we can vary them, at least when using constant
scalings of 5CB. On the other hand, the other main liquid crystal parameter
$\theta_{B}$, which describes the orientation of the molecules, affects both
evaporation (by design), and flow. This is especially true for small angles.

Rather than varying all of the Leslie viscosities, we focused on varying
$\alpha^{\prime}_{4}$, which controls the dynamic viscosity of the LL (because
$\mu_{2}=\alpha^{\prime}_{4}/2$). These results are very similar to the
Newtonian case, explored in Stapf et al. [135], in that it takes a very large
increase in $\mu_{2}$ to effect the dynamics of the AL. The director angle
does have an effect on the mobility of the LL, until $\mu_{2}$ become very
large (on the order of $10^{5}$). After this point, the LL becomes virtually
tangentially immobile for all angles of the director.

Finally, we consider two initial conditions which combine perturbations in
both the LL and the surface surfactant. The first is an undesirable case where
a valley of lipid is combined with a local excess of surfactant. As expected,
we see more breakup cases; for all angles of $\theta_{B}\lesssim 10$, TBU
occurs within 60 s as a result of both the Marangoni effect and increased
evaporation due to the deficiency in the LL.

Alternatively, we consider an initial condition which combines a bump in the
LL with a local excess of surfactant. This scenario corresponds to Section 5.6
in Bruna and Breward [100], where in a limit $h_{2}$ and $\Gamma$ can both be
governed by the same equation. While we do not take that limit here, we model
an initial condition in which the LL follows the surfactant. This area of
increased thickness in the LL seem to be protective by creating a buffer
against the Marangoni effect, and thus helping to prevent excessive thinning
in the AL. Imaging of the LL has shown areas of varying thickness on occasion
[22]. This type of initial condition could resemble this scenario, but more
work is needed to link them.

Choudhury et al. [114] modeled TBU driven by dewetting due loss of mucins in
the glycocalyx using a single layer model. They consider 119 s to be a normal
time to breakup in a healthy eye, and their simulation resulted in TBU as
short as 3.3 s when there is mucin loss over a third of the corneal surface.
Dey et al. [112, 113] studied the same mechanism for TBU, but used a
continuous-viscosity model which combined the mucus and aqueous layers. Their
model predicted breakup in 2.4 to 135 s for scenarios representing
pathological to healthy conditions, respectively. in vivo measurements for
subjects with dry eye or other eye conditions have shown TBU times well under
10 s. In healthy subjects, TBU times may be longer, but would generally occur
in less than two minutes [37, 33, 126].

While our model focuses on different mechanisms for thinning and breakup, TBU
times from our simulations fall within a similar range as previous models and
experiments. We consider that a healthy eye has high evaporative resistance,
corresponding to a large director angle. For the scenarios we considered, a
director angle of $\theta_{B}=\pi/4$, or $\pi/2$ did not result in breakup
within 60 s. We did not compute time to TBU; however, at the end of 60 s,
minimum AL thicknesses were well above the threshold value of 0.5 $\mu$m for
these larger director angles. For small angles of the director in which
evaporative resistance was low and there was a deficiency in the LL or an
excess of surfactant, our model predicted breakup in as little as 6 to 10 s.
These predicted breakup times are not as rapid as the shortest TBU times
observed in eyes [126, 33].

Thinning and breakup due to a glob surfactant was studied in Zhong et al.
[111, 30]. The glob consisted of an interval of fixed concentration as well a
part that was transported outside this interval. This area of fixed
concentration, which acted as a source, was varied in size, and found to cause
TBU either under or at the edge of the glob depending on size. TBU occurred in
$<1$ s to 15 s, and it was noted that evaporation aided TBU. Our results have
central breakup and evaporation aligns with the Marangoni effect to promote
TBU. Although our TBU time is on the longer end of the range they predicted,
our treatment of surfactant is simpler. We conserve surfactant, and the local
excess of surfactant does not act as a source; nor did we experiment with
larger or smaller concentrations.

Thinning due to evaporation has been modeled previously both with a fixed
evaporation distribution [61], and with a static LL [46]. Alternatively, Luke
et al. [54] developed an evaporation-driven model that was fit to experimental
data. We have taken a mechanistic approach, and incorporated LL thickness,
permeability, and structure in our evaporation function.

We have studied a variety of initial conditions and explored the effects of
several key parameters on tear film dynamics when modeling the LL as nematic
liquid crystal. Now that we have a sense of the key features and limitations
of this model, future work could include a moving end to imitate the eyelid
during a blink. Tear breakup is often related to the LL pattern [126, 10];
more fidelity would require a moving domain [97]. In addition, it would be
valuable to determine more rheological properties for this kind of model [72,
70, 71].

## Acknowledgements

This material is based upon work supported by the University of Delaware
Graduate College through the Doctoral Fellowship for Excellence (MJT). Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the authors. RJB is grateful for SHD’s mentoring and
many examples of the right way to do things. Remembering how he got his
unofficial middle name, we have joined words wherever a possibility of using a
hyphen has occurred in this paper.

## Declarations

### Competing interests

The authors have no competing interests to declare that are relevant to the
content of this article.

##  Appendix A Dimensional equations

We give the dimensional equations (after [136]) and a very brief outline of
the derivation.

###  A.1 Governing equations

Inside the AL, $0<y^{\prime}<h^{\prime}_{1}$, we have conservation of mass and
momentum using the Navier-Stokes equations. The advection-diffusion equation
governs osmolarity transport within the AL. The equations in dimensional form
with a prime denoting a dimensional quantity, respectively, are:

$\displaystyle\nabla^{\prime}\cdot{\bf u^{\prime}}_{1}=0,$ (29)
$\displaystyle\rho_{1}({\bf u^{\prime}}_{1t}+{\bf
u^{\prime}}_{1}\cdot\nabla^{\prime}{\bf
u^{\prime}}_{1})=-\nabla^{\prime}p^{\prime}_{1}+\mu_{1}\Delta^{\prime}{\bf
u^{\prime}}_{1},$ (30) $\displaystyle c^{\prime}_{t^{\prime}}+{\bf
u^{\prime}}_{1}\cdot\nabla^{\prime}c^{\prime}-D_{1}\Delta^{\prime}c^{\prime}=0,$
(31)

where the divergence of ${\bf u^{\prime}_{1}}$ is $\nabla^{\prime}\cdot{\bf
u^{\prime}_{1}}=u^{\prime}_{1x^{\prime}}+v^{\prime}_{1y^{\prime}}$, the
gradient is $\nabla^{\prime}{\bf
u^{\prime}_{1}}=\begin{bmatrix}u^{\prime}_{1x^{\prime}}&v^{\prime}_{1x^{\prime}}\\\
u^{\prime}_{1y^{\prime}}&v^{\prime}_{1y^{\prime}}\end{bmatrix}$, and the
Laplacian of the first component of ${\bf u^{\prime}_{1}}$ is
$\Delta^{\prime}u^{\prime}_{1}=u^{\prime}_{1x^{\prime}x^{\prime}}+u^{\prime}_{1y^{\prime}y^{\prime}}$.

Inside the LL, $h^{\prime}_{1}<y^{\prime}<h^{\prime}$, we have conservation of
energy, mass, and momentum using the Ericksen-Leslie equations for liquid
crystals.

$\displaystyle\frac{\partial}{\partial x^{\prime}}\left(\frac{\partial
W^{\prime}}{\partial\theta_{x^{\prime}}}\right)+\frac{\partial}{\partial
y^{\prime}}\left(\frac{\partial
W^{\prime}}{\partial\theta_{y^{\prime}}}\right)-\frac{\partial
W^{\prime}}{\partial\theta}+\tilde{g}^{\prime}_{x}\frac{\partial
n_{x}}{\partial\theta}+\tilde{g}^{\prime}_{y}\frac{\partial
n_{y}}{\partial\theta}$ $\displaystyle=0,$ (32)
$\displaystyle-\frac{\partial\pi^{\prime}}{\partial
x^{\prime}}+\tilde{g}^{\prime}_{x}\frac{\partial n_{x}}{\partial
x^{\prime}}+\tilde{g}^{\prime}_{y}\frac{\partial n_{y}}{\partial
x^{\prime}}+\frac{\partial\tilde{t}^{\prime}_{xy}}{\partial
y^{\prime}}+\frac{\partial\tilde{t}^{\prime}_{xx}}{\partial x^{\prime}}$
$\displaystyle=0,$ (33) $\displaystyle-\frac{\partial\pi^{\prime}}{\partial
y^{\prime}}+\tilde{g}^{\prime}_{x}\frac{\partial n_{x}}{\partial
y^{\prime}}+\tilde{g}^{\prime}_{y}\frac{\partial n_{y}}{\partial
y^{\prime}}+\frac{\partial\tilde{t}^{\prime}_{yx}}{\partial
x^{\prime}}+\frac{\partial\tilde{t}^{\prime}_{yy}}{\partial y^{\prime}}$
$\displaystyle=0,$ (34) $\displaystyle\frac{\partial u_{2}^{\prime}}{\partial
x^{\prime}}+\frac{\partial v_{2}^{\prime}}{\partial y^{\prime}}$
$\displaystyle=0,$ (35)

where, working in two dimensions,

$\displaystyle\tilde{g}^{\prime}_{i}=$
$\displaystyle-\gamma_{r}N^{\prime}_{i}-\gamma_{t}e^{\prime}_{ik}n_{k},\hskip
25.0pte^{\prime}_{ij}=\frac{1}{2}\left(\frac{\partial u^{\prime}_{i}}{\partial
x^{\prime}_{j}}+\frac{\partial u^{\prime}_{j}}{\partial
x^{\prime}_{i}}\right),$ (36) $\displaystyle N^{\prime}_{i}=$
$\displaystyle\;\dot{n}^{\prime}_{i}-\omega^{\prime}_{ik}n_{k},\hskip
55.0pt\omega^{\prime}_{ij}=\frac{1}{2}\left(\frac{\partial
u^{\prime}_{i}}{\partial x^{\prime}_{j}}-\frac{\partial
u^{\prime}_{j}}{\partial x^{\prime}_{i}}\right),\hskip
25.0pt\pi^{\prime}=\;p_{2}^{\prime}+W^{\prime},$ (37) $\displaystyle
W^{\prime}=$ $\displaystyle\;\frac{1}{2}\bigg{[}K_{1}(\nabla^{\prime}\cdot{\bf
n})^{2}+K_{3}(({\bf n}\cdot\nabla^{\prime}){\bf n})\cdot(({\bf
n}\cdot\nabla^{\prime}){\bf n})\bigg{]},$ (38)
$\displaystyle\tilde{t}^{\prime}_{ij}=$
$\displaystyle\;\alpha_{1}^{\prime}n_{k}n_{p}e^{\prime}_{kp}n_{i}n_{j}+\alpha^{\prime}_{2}N^{\prime}_{i}n_{j}+\alpha^{\prime}_{3}N^{\prime}_{j}n_{i}+\alpha^{\prime}_{4}e^{\prime}_{ij}+\alpha^{\prime}_{5}e^{\prime}_{ik}n_{k}n_{j}+\alpha^{\prime}_{6}e^{\prime}_{jk}n_{k}n_{i}.$
(39)

We use summation notation here as it is standard when working with the
Ericksen-Leslie equations, with $i,j,k=1,2$, and $\dot{n}_{i}$ denotes the
convective derivative of the $i$th component of ${\bf n}$. We further assume
that the elastic constants are equal; that is $K=K_{1}=K_{3}$. These liquid
crystal quantities are defined further in Table 6.

Table 6: Liquid crystal variables and parameters [144].  Quantity | Description
---|---
${\bf n}=(\sin\theta,\cos\theta)$ | director field
$\theta(x^{\prime},y^{\prime},t^{\prime})$ | angle the director angle makes with the $y$-axis
$p^{\prime}_{2}$ | pressure in LL
$W^{\prime}$ | bulk energy density
$\pi^{\prime}=p^{\prime}_{2}+W^{\prime}$ | modified pressure
$\tilde{g}^{\prime}_{i}$ | viscous dissipation
$e^{\prime}_{ij}$ | rate of strain tensor
$N^{\prime}_{i}$ | co-rotational time flux of the director n
$\omega^{\prime}_{ij}$ | vorticity tensor
$\tau^{V^{\prime}}_{2ij}$ | viscous stress tensor
$\tau^{E^{\prime}}_{2ij}$ | elastic stress tensor
$\tau_{2}^{\prime}$ | stress tensor
$\alpha^{\prime}_{i},\,i=1,...,6$ | Leslie viscosities (Newtonian: $\mu^{\prime}=\alpha^{\prime}_{4}/2$, all other $\alpha^{\prime}_{i}=0$)
$\gamma^{\prime}_{r}=\alpha^{\prime}_{3}-\alpha^{\prime}_{2}$ | rotational/twist viscosity
$\gamma^{\prime}_{t}=\alpha^{\prime}_{6}-\alpha^{\prime}_{5}$ | torsion coefficient
$K_{1},\;K_{3}$ | elastic constants for splay and bend respectively
|

###  A.2 Boundary conditions

We now list the boundary conditions, proceeding from the cornea out towards
the air. At $y^{\prime}=0$, we have velocity continuity and a no flux boundary
condition for osmolarity transport. These are, respectively:

$\displaystyle{\bf u^{\prime}_{1}}={\bf J_{0}},$ (40) $\displaystyle
c^{\prime}\,({\bf u^{\prime}}_{1}-{\bf
s^{\prime}}_{0t^{\prime}})-D_{1}\nabla^{\prime}c^{\prime}=0.$ (41)

Note that as we assume the cornea is flat, ${\bf
s^{\prime}}_{0t^{\prime}}={\bf 0}$.

At $y^{\prime}=h^{\prime}_{1}$, we have velocity continuity, aqueous and lipid
mass conservation, and the stress balance. In addition, we have a no flux
boundary condition for osmolarity transport, and the anchoring boundary
condition for the LL. The equations, respectively, are:

$\displaystyle({\bf u^{\prime}}_{1}-{\bf
u^{\prime}}_{2})\cdot\hat{t}^{\prime}_{1}=0,$ (42) $\displaystyle\rho_{1}({\bf
u^{\prime}}_{1}-{\bf
s^{\prime}}_{1t^{\prime}})\cdot\hat{n}^{\prime}_{1}=J^{\prime}_{e},$ (43)
$\displaystyle\rho_{2}({\bf s^{\prime}}_{1t^{\prime}}-{\bf
u^{\prime}}_{2})\cdot\hat{n}^{\prime}_{1}=0,$ (44)
$\displaystyle({\bf\tau^{\prime}}_{1}-{\bf\tau^{\prime}}_{2})\cdot\hat{n}^{\prime}_{1}=-\gamma^{\prime}_{s1}\hat{n}^{\prime}_{1}\nabla^{\prime}\cdot\hat{n}^{\prime}_{1}+\nabla^{\prime}_{s_{1}}\gamma^{\prime}_{s1},$
(45) $\displaystyle c^{\prime}\,({\bf u^{\prime}}_{1}-{\bf
s^{\prime}}_{1t^{\prime}})-D_{1}\nabla^{\prime}c^{\prime}=0,$ (46)
$\displaystyle\theta=\theta_{B}.$ (47)

To resolve the stress balance into tangential and normal stress components, we
take the dot product of the stress balance and $\hat{t}^{\prime}_{1}$ or
$\hat{n}^{\prime}_{1}$ respectively. Here $\nabla^{\prime}_{s_{1}}=({\bf
I}-\hat{n}^{\prime}_{1}\hat{n}^{\prime}_{1})\cdot\nabla^{\prime}$ is the
surface gradient, where ${\bf I}$ is the identity matrix [152].

At $y^{\prime}=h^{\prime}$, we have lipid mass conservation, stress balance,
and anchoring respectively:

$\displaystyle\rho_{2}({\bf u^{\prime}}_{2}-{\bf
s^{\prime}}_{2t^{\prime}})\cdot\hat{n}^{\prime}_{2}=0,$ (48)
$\displaystyle{\bf\tau^{\prime}}_{2}\cdot\hat{n}^{\prime}_{1}=-\gamma_{2}\hat{n}^{\prime}_{2}\nabla^{\prime}\cdot\hat{n}^{\prime}_{2},$
(49) $\displaystyle\theta=\theta_{B}.$ (50)

The stress balance can be resolved into components in the same way as above,
using the unit tangent and normal vectors at this interface,
$\hat{t}^{\prime}_{2}$ and $\hat{n}^{\prime}_{2}$.

At the aqueous-lipid interface $y^{\prime}=h^{\prime}_{1}$, we have surfactant
transport and the linear equation of state for the surface tension,
respectively:

$\displaystyle\Gamma^{\prime}_{t^{\prime}}+\nabla^{\prime}_{s_{1}}\cdot(\Gamma^{\prime}u^{\prime}_{1})=D_{s}\nabla^{{}^{\prime}2}_{s_{1}}\Gamma^{\prime},$
(51) $\displaystyle\gamma^{\prime}_{s1}=\gamma_{1}-RT_{0}\Gamma^{\prime}.$
(52)

###  A.3 Osmosis and evaporation

The effect of osmosis is determined by the difference in concentration on
either side of the aqueous-corneal interface $y^{\prime}=0$,

$\displaystyle{\bf J^{\prime}}_{0}\cdot\hat{n}^{\prime}_{0}$
$\displaystyle=P_{c}\,(c^{\prime}-c_{0}).$ (53)

Here we assume an ideal osmolarity ($c_{0}=300$ mOsM) on the corneal side of
the interface.

To account for the effect of evaporation at the aqueous-lipid interface
$y^{\prime}=h^{\prime}_{1}$, we modify a boundary layer model derived
thoroughly in Stapf et al. [135]

$\displaystyle J^{\prime}_{e}$
$\displaystyle=\frac{E_{0}}{1+\frac{k_{m}}{Dk}h^{\prime}_{2}}.$ (54)

Here, $Dk$ represents resistance to evaporation provided by the LL based on
its thickness and permeability. Its value comes from a nonlinear least squares
fit to clinical data [34] performed by Stapf et al. [135], which is very
similar to that of Cerretani et al. [123]. We further modify $Dk$ to include a
dependence on the director angle of liquid crystals in the LL. We suggest the
orientation of the molecules affects evaporation through

$\displaystyle J^{\prime}_{e}$
$\displaystyle=\frac{E_{0}}{1+\frac{k_{m}}{Dk}(0.1+0.9\sin\theta_{B})h^{\prime}_{2}}.$
(55)

The maximum evaporation rate occurs when $\theta_{B}=0$, which is ten times
that of the rate for $\theta_{B}=\pi/2$. The minimum rate matches that used in
[135].

###  A.4 Deriving the leading order system

We nondimensionalize using the scalings in (1), along with

$\displaystyle W^{\prime}=\frac{K}{\delta^{2}L^{2}}W,\hskip
25.0pt\alpha_{i}^{\prime}=\mu_{2}\alpha_{i},\hskip
25.0pt\gamma^{\prime}_{t}=\mu_{2}\gamma_{t},\hskip
25.0pt\gamma^{\prime}_{r}=\mu_{2}\gamma_{r}.$

Then, we asymptotically expand all the dependent variables in powers of the
small parameter $\epsilon\ll 1$. For example,

$\Gamma(x,t)\sim\Gamma_{0}(x,t)+\epsilon\,\Gamma_{1}(x,t)+\epsilon^{2}\Gamma_{2}(x,t)+\cdots$

We substitute these expansions into the nondimensionalized equations and
collect like powers of $\epsilon$. At leading order, we are able to determine
the director angle, conservation of mass in both the AL and LL, and
conservation of surfactant concentration. At $O(\epsilon)$, we are only able
to determine intermediate variables, so we proceed to $O(\epsilon^{2})$. We
are then able to determine the AL depth averaged velocity, AL pressure, the
force balance equation for the axial velocity in the LL, and conservation of
osmolarity. This gives us the closed system of equations from Section 2.1.
This progression is typical of problems like this [135, 100].

## References

  *   * Bron et al. [2004] Bron, A.J., Tiffany, J.M., Gouveia, S.M., Yokoi, N., Voon, L.W.: Functional aspects of the tear film lipid layer. Exp. Eye Res. 78, 347–360 (2004)
  * Smith et al. [2007] Smith, J.A., et al.: The epidemiology of dry eye disease: Report of the Epidemiology Subcommittee of the International Dry Eye WorkShop. Ocul. Surf. 5, 93–107 (2007)
  * Stapleton et al. [2017] Stapleton, F., Alves, M., Bunya, V.Y., Jalbert, I., Lekhanont, K., Malet, F., Na, K.-S., Schaumberg, D., Uchino, M., Vehof, J., Viso, E., Vitale, S., Jones, L.: TFOS DEWS-II epidemiology report. Ocul. Surf. 15(3), 334–365 (2017)
  * Lemp [1995] Lemp, M.A.: Report of the National Eye Institute/Industry workshop on Clinical Trials in Dry Eyes. CLAO J 21, 221–232 (1995)
  * Lemp et al. [2007] Lemp, M.A., et al.: The definition and classification of dry eye disease: Report of the Definition and Classification Subcommittee of the International Dry Eye WorkShop. Ocul. Surf. 5, 75–92 (2007)
  * Craig et al. [2017] Craig, J.P., Nichols, K.K., Akpek, E.K., Caffery, B., Dua, H.S., Joo, C.-K., Liu, Z., Nelson, J.D., Nichols, J.J., Tsubota, K., et al.: TFOS DEWS-II definition and classification report. Ocul. Surf. 15(3), 276–283 (2017)
  * Willcox et al. [2017] Willcox, M.D.P., Arguëso, P., Georgiev, G.A., Holopainen, J.M., Laurie, G.W., Millar, T.J., Papas, E.B., Rolland, J.P., Schmidt, T.A., Stahl, U., Suarez, T., Subbaraman, L.N., Ucakhan, O.O., Jones, L.W.: TFOS DEWS-II Tear Film Report. Ocul. Surf. 15(3), 369–406 (2017)
  * Baudouin et al. [2013] Baudouin, C., Aragona, P., Messmer, E.M., Tomlinson, A., Calonge, M., Boboridis, K.G., Akova, Y.A., Geerling, G., Labetoulle, M., Rolando, M.: Role of hyperosmolarity in the pathogenesis and management of dry eye disease: Proceedings of the OCEAN group meeting. Ocul. Surf. 11, 246–258 (2013)
  * Tsubota [1998] Tsubota, K.: Tear dynamics and dry eye. Prog. Retin. Eye Res. 17(4), 565–596 (1998)
  * Braun et al. [2015] Braun, R.J., King-Smith, P.E., Begley, C.G., Li, L., Gewecke, N.R.: Dynamics and function of the tear film in relation to the blink cycle. Prog. Retin. Eye Res. 45, 132–164 (2015)
  * Mishima and Maurice [1961] Mishima, S., Maurice, D.M.: The oily layer of the tear film and evaporation from the corneal surface. Exp. Eye Res. 1, 39–45 (1961)
  * Svitova and Lin [2021] Svitova, T.F., Lin, M.C.: Evaporation retardation by model tear-lipid films: The roles of film aging, compositions and interfacial rheological properties. Coll. Surf. B: Biointerfaces 197, 111392 (2021)
  * Norn [1979] Norn, M.S.: Semiquantitative interference study of the fatty layer of precorneal film. Acta Ophthalmol. 57, 766–774 (1979)
  * Berger and Corrsin [1974] Berger, R.E., Corrsin, S.: A surface tension gradient mechanism for driving the pre-corneal tear film after a blink. J. Biomech. 7, 225–38 (1974)
  * Craster and Matar [2009] Craster, R.V., Matar, O.K.: Dynamics and stability of thin liquid films. Rev. Mod. Phys. 81(3), 1131 (2009)
  * Holly [1973] Holly, F.J.: Formation and rupture of the tear film. Exp. Eye Res. 15, 515–525 (1973)
  * Tiffany [1990a] Tiffany, J.M.: Measurement of wettability of the corneal epithelium: I. particle attachment method. Acta Ophthalmol. 68(2), 175–181 (1990)
  * Tiffany [1990b] Tiffany, J.M.: Measurement of wettability of the corneal epithelium: II. contact angle method. Acta Ophthalmol. 68(2), 182–187 (1990)
  * Gipson [2004] Gipson, I.K.: Distribution of mucins at the ocular surface. Exp. Eye Res. 78, 379–388 (2004)
  * King-Smith et al. [2004] King-Smith, P.E., Fink, B.A., Hill, R.M., Koelling, K.W., Tiffany, J.M.: The thickness of the tear film. Curr. Eye. Res. 29, 357–368 (2004)
  * Yokoi et al. [1996] Yokoi, N., Takehisa, Y., Kinoshita, S.: Correlation of tear lipid layer interference patterns with the diagnosis and severity of dry eye. Am. J. Ophthalmol. 122, 818–824 (1996)
  * Goto and Tseng [2003] Goto, E., Tseng, S.C.G.: Kinetic analysis of tear interference images in aqueous tear deficiency dry eye before and after punctal occlusion. Invest. Ophthalmol. Vis. Sci. 44, 1897–1905 (2003)
  * King-Smith et al. [2011] King-Smith, P.E., Nichols, J.J., Nichols, K.K., Braun, R.J.: A high resolution microscope for imaging the lipid layer of the tear film. Ocul. Surf. 9(4), 197–211 (2011)
  * Govindarajan and Gipson [2010] Govindarajan, B., Gipson, I.K.: Membrane-tethered mucins have multiple functions on the ocular surface. Exp. Eye Res. 90, 655–693 (2010)
  * Owens and Phillips [2001] Owens, H., Phillips, J.: Spreading of the tears after a blink: Velocity and stabilization time in healthy eyes. Cornea 20, 484–487 (2001)
  * King-Smith et al. [2008] King-Smith, P.E., Nichols, J.J., Nichols, K.K., Fink, B.A., Braun, R.J.: Contributions of evaporation and other mechanisms to tear film thinning and breakup. Optom. Vis. Sci. 85, 623–630 (2008)
  * Norn [1969] Norn, M.S.: Dessication of the precorneal film. I Corneal wetting-time. Acta Ophthalmol. 47, 865–80 (1969)
  * Lemp et al. [2007] Lemp, M.A., et al.: The definition and classification of dry eye disease: Report of the Definition and Classification Subcommittee of the International Dry Eye WorkShop. Ocul. Surf. 5, 75–92 (2007)
  * King-Smith et al. [2018] King-Smith, P.E., Begley, C.G., Braun, R.J.: Mechanisms, imaging and structure of tear film breakup. Ocul. Surf. 16, 4–30 (2018)
  * Zhong et al. [2019] Zhong, L., Braun, R.J., Begley, C.G., King-Smith, P.E.: Dynamics of fluorescent imaging for rapid tear thinning. Bull. Math. Biol. 81, 39–80 (2019)
  * Zhang et al. [2003] Zhang, Y.L., Matar, O.K., , Craster, R.V.: Analysis of tear film rupture: Effect of non-newtonian rheology. J. Colloid Interface Sci. 262, 130–148 (2003)
  * Zhang et al. [2004] Zhang, Y.L., Matar, O.K., , Craster, R.V.: Rupture analysis of the corneal mucus layer of the tear film. Mol. Simul. 30, 167–172 (2004)
  * Yokoi and Georgiev [2019] Yokoi, N., Georgiev, G.A.: Tear-film-oriented diagnosis for dry eye. Japan. J. Ophthalmol. 63, 127–136 (2019)
  * King-Smith et al. [2010] King-Smith, P.E., Hinel, E.A., Nichols, J.J.: Application of a novel interferometric method to investigate the relation between lipid layer thickness and tear film thinning. Invest. Ophthalmol. Vis. Sci. 51, 2418–23 (2010)
  * Sharma and Ruckenstein [1985] Sharma, A., Ruckenstein, E.: Mechanism of tear film rupture and formation of dry spots on cornea. J. Coll. Interface Sci. 106, 12–27 (1985)
  * Sharma et al. [1999] Sharma, A., Khanna, R., Reiter, G.: A thin film analog of the corneal mucus layer of the tear film: an enigmatic long range non-classical DLVO interaction in the breakup of thin polymer films. Coll. Surf. B 14, 223–35 (1999)
  * Begley et al. [2013] Begley, C.G., Simpson, T., Liu, H., Salvo, E., Wu, Z., Bradley, A., Situ, P.: Quantitative Analysis of Tear Film Fluorescence and Discomfort During Tear Film Instability and Thinning. Invest. Ophthalmol. Vis. Sci. 54(4), 2645–2653 (2013)
  * King-Smith et al. [2013] King-Smith, P.E., Bailey, M.D., Braun, R.J.: Four characteristics and a model of an effective tear film lipid layer. Ocul. Surf. 11, 1–10 (2013)
  * Luke et al. [2021] Luke, R.A., Braun, R.J., Driscoll, T.A., Awisi-Gyau, D., Begley, C.G.: Parameter Estimation for Mixed-Mechanism Tear Film Thinning. Bull. Math. Biol. 83(5), 56 (2021)
  * Doane [1981] Doane, M.G.: Blinking and the mechanics of the lacrimal drainage system. Ophthalmol. 88, 844–851 (1981)
  * Tomlinson et al. [2009] Tomlinson, A., Doane, M.G., McFayden, A.: Inputs and outputs of the lacrimal system: Review of production and evaporative loss. Ocul. Surf. 7, 17–29 (2009)
  * Kimball et al. [2010] Kimball, S.H., King-Smith, P.E., Nichols, J.J.: Evidence for the major contribution of evaporation to tear film thinning between blinks. Invest. Opthalmol. Vis. Sci. 51, 6294–97 (2010)
  * Braun [2012] Braun, R.J.: Dynamics of the tear film. Annu. Rev. Fluid Mech. 44, 267–297 (2012)
  * Cerretani and Radke [2014] Cerretani, C.F., Radke, C.J.: Tear dynamics in healthy and dry eyes. Curr. Eye Res. 39, 580–595 (2014)
  * Stahl et al. [2012] Stahl, U., Willcox, M., Stapleton, F.: Osmolality and tear film dynamics. Clin. Exp. Optom. 95, 3–11 (2012)
  * Peng et al. [2014] Peng, C.-C., Cerretani, C., Braun, R.J., Radke, C.J.: Evaporation-driven instability of the precorneal tear film. Adv. Colloid Interface Sci. 206, 250–264 (2014)
  * Lemp et al. [2011] Lemp, M.A., Bron, A.J., Baudouin, C., Castillo, J.M.B., Geffen, D., Tauber, J., Foulks, G.N., Pepose, J.S., Sullivan, B.D.: Tear osmolarity in the diagnosis and management of dry eye disease. Am. J. Ophthalmol. 151(5), 792–798 (2011)
  * Tomlinson et al. [2006] Tomlinson, A., Khanal, S., Ramaesh, K., Diaper, C., McFadyen, A.: Tear film osmolarity: determination of a referent for dry eye diagnosis. Invest. Ophthalmol. Vis. Sci. 47(10), 4309–4315 (2006)
  * Versura et al. [2010] Versura, P., Profazio, V., Campos, E.: Performance of tear osmolarity compared to previous diagnostic tests for dry eye diseases. Curr. Eye Res. 35(7), 553–564 (2010)
  * Gilbard et al. [1978] Gilbard, J.P., Farris, R.L., Santamaria, J.: Osmolarity of tear microvolumes in keratoconjunctivitis sicca. Arch. Ophthalmol. 96, 677–681 (1978)
  * Sullivan et al. [2010] Sullivan, B.D., Whitmer, D., Nichols, K.K., Tomlinson, A., Foulks, G.N., Geerling, G., Pepose, J.S., Kosheleff, V., Porreco, A., Lemp, M.A.: An objective approach to dry eye disease severity. Invest. Ophthalmol. Vis. Sci. 51(12), 6125–6130 (2010)
  * Dartt and Willcox [2013] Dartt, D., Willcox, M.: Complexity of the tear film: importance in homeostasis and dysfunction during disease. Exp. Eye Res. 117, 1–3 (2013)
  * Liu et al. [2009] Liu, H., Begley, C., Chen, M., Bradley, A., Bonanno, J., McNamara, N.A., Nelson, J.D., Simpson, T.: A link between tear instability and hyperosmolarity in dry eye. Invest. Ophthalmol. Vis. Sci. 50(8), 3671–3679 (2009)
  * Luke et al. [2020] Luke, R.A., Braun, R.J., Driscoll, T.A., Begley, C.G., Awisi-Gyau, D.: Parameter Estimation for Evaporation-Driven Tear Film Thinning. Bull. Math. Biol. 82(6), 71 (2020)
  * Pflugfelder [2011] Pflugfelder, S.C.: Tear dysfunction and the cornea: LXVIII Edward Jackson memorial lecture. Am. J. Ophthalmol. 152(6), 900–909 (2011)
  * Belmonte et al. [2015] Belmonte, C., Acosta, M.C., Merayo-Lloves, J., Gallar, J.: What causes eye pain? Curr. Ophthalmol. Rep. 3(2), 111–121 (2015)
  * Nichols et al. [2005] Nichols, J.J., Mitchell, G.L., King-Smith, P.E.: Thinning rate of the precorneal and prelens tear films. Invest. Ophthalmol. Vis. Sci. 46(7), 2353–2361 (2005)
  * Hamano et al. [1981] Hamano, H., Hori, M., Mitsunaga, S.: Measurement of evaporation rate of water from the precorneal tear film and contact lenses. Contacto 25(2), 7–15 (1981)
  * Peng et al. [2014] Peng, C.-C., Cerretani, C., Li, Y., Bowers, S., Shahsavarani, S., Lin, M., Radke, C.: Flow evaporimeter to assess evaporative resistance of human tear-film lipid layer. Ind. Eng. Chem. Res. 53(47), 18130–18139 (2014)
  * Dursch et al. [2018] Dursch, T.J., Li, W., Taraz, B., Lin, M.C., Radke, C.J.: Tear-film evaporation rate from simultaneous ocular-surface temperature and tear-breakup area. Optom. Vis. Sci. 95, 5–12 (2018)
  * Braun et al. [2018] Braun, R.J., Driscoll, T.A., Begley, C.G., King-Smith, P.E., Siddique, J.I.: On tear film breakup (TBU): Dynamics and imaging. Math. Med. Biol. 45, 145–180 (2018)
  * Tiffany and Dart [1981] Tiffany, J.M., Dart, J.K.G.: Normal and abnormal functions of meibomian secretion. R. Soc. Med. Int. Congr. Symp. 40, 1061–4 (1981)
  * Purslow and Wolffsohn [2005] Purslow, C., Wolffsohn, J.S.: Ocular surface temperature: A review. Eye Contact Lens 31, 117–123 (2005)
  * Tiffany [1994] Tiffany, J.M.: Viscoelastic properties of human tears and polymer solutions. In: Sullivan, D.A. (ed.) Lacrimal Gland, Tear Film, and Dry Eye Syndromes, pp. 267–70. Plenum Press, New York (1994)
  * Pandit et al. [1999] Pandit, J.C., Nagyová, B., Bron, A.J., Tiffany, J.M.: Physical properties of stimulated and unstimulated tears. Exp. Eye Res. 68, 247–53 (1999)
  * Tiffany [1991] Tiffany, J.M.: The viscosity of human tears. Intl. Ophthalmol. 15, 371–376 (1991)
  * Jossic et al. [2009] Jossic, L., Lefevre, P., Loubens, C., Magnin, A., Corre, C.: The fluid mechanics of shear-thinning tear substitutes. J. Nonnewton. Fluid Mech. 161, 1–9 (2009)
  * Tiffany and Nagyová [1998] Tiffany, J.M., Nagyová, B.: Components of tears responsible for surface tension. In: Sullivan, D.A., Dartt, D.A., Meneray, M.A. (eds.) Lacrimal Gland, Tear Film, and Dry Eye Syndromes 2, pp. 581–5. Plenum Press, New York (1998)
  * Tiffany [1987] Tiffany, J.M.: The lipid secretion of the meibomian glands. Adv. Lipid Res. 22, 1–62 (1987)
  * Leiske et al. [2011] Leiske, D.L., Leiske, C.I., Leiske, D.R., Toney, M.F., Senchyna, M., Ketelson, H.A., Meadows, D.L., Fuller, G.G.: Temperature-induced transitions in the structure and interfacial rheology of human meibum. Biophys. J. 102, 369–376 (2011)
  * Leiske et al. [2012] Leiske, D.L., Miller, C.E., Rosenfeld, L., Cerretani, C., Ayzner, A., Lin, B., Meron, M., Senchyna, M., Ketelson, H.A., Meadows, D., Srinivasan, S., Jones, L., Radke, C.J., Toney, M.F., Fuller, G.G.: Molecular structure of interfacial human meibum films. Langmuir 28, 11858–11865 (2012)
  * Rosenfeld et al. [2013] Rosenfeld, L., Cerretani, C., Leiske, D.L., Toney, M.F., Radke, C.J., Fuller, G.G.: Structural and rheological properties of meibomian lipid. Invest. Ophthalmol. Vis. Sci. 54(4), 2720–2732 (2013)
  * Miller [1969] Miller, D.: Measurement of the surface tension of tears. Arch. Ophthalmol. 82, 368–71 (1969)
  * Nagyová and Tiffany [1999] Nagyová, B., Tiffany, J.M.: Components of tears responsible for surface tension. Curr. Eye Res. 19, 4–11 (1999)
  * Mudgil and Millar [2006] Mudgil, P., Millar, T.J.: Adsorption of lysozyme to phospholipid and Meibomian lipid monolayer films. Coll. Surf. B: Biointerfaces 48, 128–37 (2006)
  * Mudgil and Millar [2008] Mudgil, P., Millar, T.J.: The effect of Meibomian lipids and tear proteins on evaporation rate under controlled in vitro conditions. Exp. Eye Res. 86, 622–8 (2008)
  * Butovich [2013] Butovich, I.A.: Tear film lipids. Experimental Eye Research 117, 4–27 (2013)
  * Rantamäki and Holopainen [2017] Rantamäki, A.H., Holopainen, J.M.: The Effect of Phospholipids on Tear Film Lipid Layer Surface Activity. Invest. Ophthalmol. Vis. Sci. 58, 149–154 (2017)
  * Galor et al. [2022] Galor, A., Sanchez, V., Jensen, A., Burton, M., Maus, K., Stephenson, D., Chalfant, C., Mandal, N.: Meibum sphingolipid composition is altered in individuals with meibomian gland dysfunction-a side by side comparison of meibum and tear sphingolipids. The Ocular Surface 23, 87–95 (2022)
  * Georgiev et al. [2017] Georgiev, G.A., Eftimov, P., Yokoi, N.: Structure-function relationship of tear film lipid layer: A contemporary perspective. Exp. Eye Res. 163, 17–28 (2017)
  * Braun et al. [2019] Braun, R.J., Driscoll, T.A., Begley, C.G.: Mathematical models of the tear film. In: Ocular Fluid Dynamics, pp. 387–432. Springer, Berlin (2019)
  * Braun et al. [2012] Braun, R.J., Usha, R., McFadden, G.B., Driscoll, T.A., Cook, L.P., King-Smith, P.E.: Thin film dynamics on a prolate spheroid with application to the cornea. J. Eng. Math. 73, 121–138 (2012)
  * Maki et al. [2010a] Maki, K.L., Braun, R.J., Henshaw, W.D., King-Smith, P.E.: Tear film dynamics on an eye-shaped domain I: pressure boundary conditions. Math. Med. Biol. 27, 227–254 (2010)
  * Maki et al. [2010b] Maki, K.L., Braun, R.J., Ucciferro, P., Henshaw, W.D., King-Smith, P.E.: Tear film dynamics on an eye-shaped domain. Part 2. Flux boundary conditions. J. Fluid Mech. 647, 361–390 (2010)
  * Li et al. [2014] Li, L., Braun, R.J., Maki, K.L., Henshaw, W.D., King-Smith, P.E.: Tear film dynamics with evaporation, wetting and time-dependent flux boundary condition on an eye-shaped domain. Phys. Fluids 26, 052101 (2014)
  * Li et al. [2016] Li, L., Braun, R.J., Henshaw, W.D., King-Smith, P.E.: Computed tear film and osmolarity dynamics on an eye-shaped domain. Math. Med. Biol. 33, 123–157 (2016)
  * Li et al. [2018] Li, L., Braun, R.J., Henshaw, W.D., King-Smith, P.E.: Computed flow and fluorescence over the ocular surface. Math. Med. Biol. 35, 51–85 (2018)
  * Driscoll et al. [2018] Driscoll, T.A., Braun, R.J., Brosch, J.K.: Simulation of parabolic flow on an eye-shaped domain with moving boundary. J. Eng. Math. 111(1), 111–126 (2018)
  * Maki et al. [2019] Maki, K.L., Henshaw, W., McManus, A., Braun, R., Chapp, D., Driscoll, T.: A model for tear film dynamics during a realistic blink. Modeling and Artificial Intelligence in Ophthalmology 2(3), 21–27 (2019)
  * Gaffney et al. [2010] Gaffney, E.A., Tiffany, J.M., Yokoi, N., Bron, A.J.: A mass and solute balance model for tear volume and osmolarity in the normal and the dry eye. Prog. Retinal Eye Res. 29, 59–78 (2010)
  * Sharma et al. [1998] Sharma, A., Tiwari, S., Khanna, R., Tiffany, J.M.: Hydrodynamics of Meniscus-induced Thinning of the Tear Film. In: Sullivan, D.A., Dartt, D.A., Meneray, M.A. (eds.) Advances in Experimental Medicine and Biology. Lacrimal Gland, Tear Film, and Dry Eye Syndromes 2, vol. 438, pp. 425–431. Springer, Berlin (1998)
  * Wong et al. [1996] Wong, H., Fatt, I., Radke, C.: Deposition and thinning of the human tear film. J. Colloid Interface Sci. 184(1), 44–51 (1996)
  * Miller et al. [2002] Miller, K.L., Polse, K.A., Radke, C.J.: Black line formation and the “perched” human tear film. Curr. Eye Res. 25, 155–62 (2002)
  * Braun and Fitt [2003] Braun, R.J., Fitt, A.D.: Modeling the drainage of the precorneal tear film after a blink. Math. Med. Biol. 20, 1–28 (2003)
  * Li and Braun [2012] Li, L., Braun, R.J.: A model for the human tear film with heating from within the eye. Phys. Fluids 24, 062103 (2012)
  * Jones et al. [2005] Jones, M.B., Please, C.P., McElwain, D.L.S., Fulford, G.R., Roberts, A.P., Collins, M.J.: Dynamics of tear film deposition and draining. Math. Med. Biol. 22, 265–88 (2005)
  * Jones et al. [2006] Jones, M.B., McElwain, D.L.S., Fulford, G.R., Collins, M.J., Roberts, A.P.: The effect of the lipid layer on tear film behavior. Bull. Math. Biol. 68, 1355–81 (2006)
  * Maki et al. [2008] Maki, K.L., Braun, R.J., Driscoll, T.A., King-Smith, P.E.: An overset grid method for the study of reflex tearing. Math. Med. Biol. 25, 187–214 (2008)
  * Aydemir et al. [2010] Aydemir, E., Breward, C.J.W., Witelski, T.P.: The effect of polar lipids on tear film dynamics. Bull. Math. Biol. 73, 1171–1201 (2010)
  * Bruna and Breward [2014] Bruna, M., Breward, C.J.W.: The influence of non-polar lipids on tear film dynamics. J. Fluid Mech. 746, 565–605 (2014)
  * Allouche et al. [2017] Allouche, M., Abderrahmane, H.A., Djouadi, S.M., Mansouri, K.: Influence of curvature on tear film dynamics. Eur. J. Mech. B Fluids 66, 81–91 (2017)
  * Mehdaoui et al. [2021] Mehdaoui, H., Abderrahmane, H.A., Bouda, F.N., Koulali, A.: 2D numerical simulation of tear film dynamics: effects of shear-thinning properties. Eur. J. Mech. B Fluids 90, 128–136 (2021)
  * Braun and King-Smith [2007] Braun, R.J., King-Smith, P.E.: Model problems for the tear film in a blink cycle: Single equation models. J. Fluid Mech. 586, 465–90 (2007)
  * Heryudono et al. [2007] Heryudono, A., Braun, R.J., Driscoll, T.A., Maki, K.L., Cook, L.P., King-Smith, P.E.: Single-equation models for the tear film in a blink cycle: realistic lid motion. Math. Med. Biol. J. IMA 24(4), 347–377 (2007)
  * Zubkov et al. [2012] Zubkov, V.S., Breward, C.J.W., Gaffney, E.A.: Coupling fluid and solute dynamics within the ocular surface tear film: a modelling study of black line osmolarity. Bull. Math. Biol. 74(9), 2062–2093 (2012)
  * Deng et al. [2014] Deng, Q., Braun, R.J., Driscoll, T.A.: Heat transfer and tear film dynamics over multiple blink cycles. Phys. Fluids 26(7), 071901 (2014)
  * Deng et al. [2013] Deng, Q., Braun, R.J., Driscoll, T.A., King-Smith, P.E.: A model for the tear film and ocular surface temperature for partial blinks. Interfacial Phenom. Heat Transf. 1(4) (2013)
  * Luke et al. [2021] Luke, R.A., Braun, R.J., Begley, C.G.: Mechanistic determination of tear film thinning via fitting simplified models to tear breakup. Model. Artif. Int. Ophthalmol. 3, 71–100 (2021)
  * Driscoll et al. [2023] Driscoll, T.A., Braun, R.J., Luke, R.A., Sinopoli, D., Phatak, A., Dorsch, J., Begley, C.G., Awisi-Gyau, D.: Fitting ODE models of tear film breakup. Model. Art. Int. Ophthal. 5(1), 1–36 (2023)
  * Zhong et al. [2018a] Zhong, L., Ketelaar, C.F., Braun, R.J., Begley, C.G., King-Smith, P.E.: Mathematical modeling of glob-driven tear film breakup. Math. Med. Biol., 1–37 (2018)
  * Zhong et al. [2018b] Zhong, L., Braun, R.J., King-Smith, P.E., Begley, C.G.: Mathematical modeling of glob-driven tear film breakup. J. Model. Ophthalmol. 2(1), 24–28 (2018)
  * Dey et al. [2019] Dey, M., Vivek, A.S., Dixit, H.N., Richhariya, A., Feng, J.J.: A model of tear-film breakup with continuous mucin concentration and viscosity profiles. J. Fluid Mech. 858, 352–376 (2019)
  * Dey et al. [2020] Dey, M., Vivek, A.S., Dixit, H.N., Richhariya, A., Feng, J.J.: A model of tear-film breakup with continuous mucin concentration and viscosity profiles–Corrigendum. J. Fluid Mech. 889, 1 (2020)
  * Choudhury et al. [2021] Choudhury, A., Dey, M., Dixit, H.N., Feng, J.J.: Tear-film breakup: The role of membrane-associated mucin polymers. Phys. Rev. E 103, 013018 (2021)
  * Ji and Witelski [2017] Ji, H., Witelski, T.P.: Finite-time thin film rupture driven by modified evaporative loss. Physica D: Nonlinear Phenomena 342, 1–15 (2017)
  * Ji and Witelski [2020] Ji, H., Witelski, T.P.: Steady states and dynamics of a thin-film-type equation with non-conserved mass. Eur. J. Appl. Math. 31, 968–1001 (2020)
  * Shi et al. [2021] Shi, X., Rodríguez-Hakim, M., Shaqfeh, E.S.G., Fuller, G.G.: Instability and symmetry breaking in binary evaporating thin films over a solid spherical dome. J. Fluid Mech. 915, 45 (2021)
  * Shi et al. [2022] Shi, X., Fuller, G.G., Shaqfeh, E.S.G.: Instability and symmetry breaking of surfactant films over an air bubble. J. Fluid Mech. 953, 26 (2022)
  * Pucker and Haworth [2015] Pucker, A., Haworth, K.: The presence and significance of polar meibum and tear lipids. Ocul. Surf. 13(1), 26–42 (2015)
  * Butovich [2017] Butovich, I.A.: Meibomian glands, meibum, and meibogenesis. Exp. Eye Res. 163, 2–16 (2017)
  * Georgiev et al. [2019] Georgiev, G.A., Borchman, D., Eftimov, P., Yokoi, N.: Lipid saturation and the rheology of human tear lipids. Int. J. Molec. Sci. 20, 3431 (2019)
  * McCulley and Shine [1997] McCulley, J.P., Shine, W.: A compositional based model for the tear film lipid layer. Tr. Am. Ophthalmol. Soc. XCV, 79–93 (1997)
  * Cerretani et al. [2013] Cerretani, C.F., Ho, N.H., Radke, C.: Water-evaporation reduction by duplex films: Application to the human tear film. Adv. Colloid Interface Sci. 197, 33–57 (2013)
  * Sledge et al. [2016] Sledge, S.M., Khimji, H., Borchman, D., Oliver, A.L., Michael, H., Dennis, E., Gerlach, D., Bhola, R., Stephen, E.: Evaporation and hydrocarbon chain conformation of surface lipid films. Ocul. Surf. 14(4), 447–459 (2016)
  * Xu et al. [2023] Xu, X., Li, G., Zuo, Y.Y.: Effect of Model Tear Film Lipid Layer on Water Evaporation. Invest. Ophthalmol. Vis. Sci. 64(1), 13 (2023)
  * King-Smith et al. [2013] King-Smith, P.E., Reuter, K.S., Braun, R.J., Nichols, J.J., Nichols, K.K.: Tear film breakup and structure studied by simultaneous video recording of fluorescence and tear film lipid layer, TFLL, images. Invest. Ophthalmol. Vis. Sci. 54, 4900–4909 (2013)
  * Fenner and Tong [2015] Fenner, B.J., Tong, L.: More to Stable Tears Than Thickness of the Tear Film Lipid Layer. Invest. Ophthalmol. Vis. Sci 56(3), 1601–1601 (2015)
  * King-Smith and Braun [2015] King-Smith, P.E., Braun, R.J.: Author Response: More to Stable Tears Than Thickness of the Lipid Layer. Invest. Ophthalmol. Vis. Sci 56(3), 1602–1602 (2015)
  * Wizert et al. [2014] Wizert, A., Iskander, D., Cwiklik, L.: Organization of lipids in the tear film: A molecular-level view. PLoS ONE 9, 92461 (2014)
  * Cwiklik [2016] Cwiklik, L.: Tear film lipid layer: A molecular level view. Biochimica et Biophysica Acta 1858, 2421–2430 (2016)
  * Bland et al. [2019] Bland, H.C., Moilanen, J.A., Ekholm, F.S., Paananen, R.O.: Investigating the role of specific tear film lipids connected to dry eye syndrome: A study on o-acyl-$\omega$-hydroxy fatty acids and diesters. Langmuir 35(9), 3545–3552 (2019)
  * Viitaja et al. [2021] Viitaja, T., Moilanen, J., Svedström, K.J., Ekholm, F.S., Paananen, R.O.: Tear film lipid layer structure: Self-assembly of o-acyl-$\omega$-hydroxy fatty acids and wax esters into evaporation-resistant monolayers. Nano Letters 21(18), 7676–7683 (2021)
  * Paananen et al. [2019] Paananen, R.O., Javanainen, M., Holopainen, J.M., Vattulainen, I.: Crystalline wax esters regulate the evaporation resistance of tear film lipid layers associated with dry eye syndrome. J. Phys. Chem. Lett. 10(14), 3893–3898 (2019)
  * Butovich et al. [2014] Butovich, I.A., Lu, H., McMahon, A., Ketelson, H., Senchyna, M., Meadows, D., Campbell, E., Molai, M., Linsenbardt, E.: Biophysical and morphological evaluation of human normal and dry eye meibum using hot stage polarized light microscopy. Invest. Ophthalmol. Vis. Sci. 55(1), 87–101 (2014)
  * Stapf et al. [2017] Stapf, M.R., Braun, R.J., King-Smith, P.E.: Duplex tear film evaporation analysis. Bull. Math. Biol. 79(12), 2814–2846 (2017)
  * Taranchuk and Braun [2024] Taranchuk, M.J., Braun, R.J.: On the effect of liquid crystal orientation in the lipid layer on tear film thinning and breakup. ArXiV (2404.13225 [physics.flu-dyn]) (2024)
  * Matar et al. [2002] Matar, O.K., Craster, R.V., Warner, M.R.E.: Surfactant transport on highly viscous surface films. J. Fluid Mech. 466, 85–111 (2002)
  * Cummings et al. [2014] Cummings, L.J., Low, J., Myers, T.G.: Extensional flow of nematic liquid crystal with an applied electric field. Eur. J. Appl. Math. 25(4), 397–423 (2014)
  * Taranchuk et al. [2023] Taranchuk, M.J., Cummings, L.J., Driscoll, T.A., Braun, R.J.: Extensional flow of a free film of nematic liquid crystal with moderate elasticity. Phys. Fluids 35(6) (2023)
  * King-Smith et al. [2010] King-Smith, P.E., Hinel, E.A., Nichols, J.J.: Application of a novel interferometric method to investigate the relation between lipid layer thickness and tear film thinning. Invest. Ophthalmol. Vis. Sci. 51, 2418–2423 (2010)
  * King-Smith et al. [2011] King-Smith, P.E., Nichols, J.J., Braun, R.J., Nichols, K.K.: High resolution microscopy of the lipid layer of the tear film. Ocul. Surf. 9(4), 197–211 (2011)
  * Goto and Tseng [2003] Goto, E., Tseng, S.C.G.: Kinetic analysis of tear interference images in aqueous tear deficiency dry eye before and after punctal occlusion. Invest. Ophthalmol. Vis. Sci. 44, 1897–1905 (2003)
  * King-Smith et al. [2004] King-Smith, E., Fink, B., Hill, R., Koelling, K., Tiffany, J.: The thickness of the tear film. Curr. Eye Res. 29(4-5), 357–368 (2004)
  * Stewart [2019] Stewart, I.W.: The Static and Dynamic Continuum Theory of Liquid Crystals: a Mathematical Introduction. CRC Press, Boca Raton, FL (2019)
  * Weast and Astle [1978] Weast, R.C., Astle, M.J.: CRC Handbook of Chemistry and Physics: A Ready-reference Book of Chemical and Physical Data. CRC Press, Boca Raton, FL (1978)
  * Tiffany [1987] Tiffany, J.M.: The lipid secretion of the meibomian glands. Adv. Lipid Res. 22, 1–62 (1987)
  * Nagyova and Tiffany [1999] Nagyova, B., Tiffany, J.M.: Components responsible for the surface tension of human tears. Curr. Eye Res. 19(1), 4–11 (1999)
  * Riquelme et al. [2007] Riquelme, R., Lira, I., Pérez-López, C., Rayas, J.A., Rodríguez-Vera, R.: Interferometric measurement of a diffusion coefficient: comparison of two methods and uncertainty analysis. J. Phys. D: Appl. Phys. 40(9), 2769 (2007)
  * Braun et al. [2000] Braun, R.J., Naire, S., Snow, S.A.: Limiting cases of gravitational drainage of a vertical free film for evaluating surfactants. SIAM J. Appl. Math. 61, 889–913 (2000)
  * Stebe and Maldarelli [1994] Stebe, K.J., Maldarelli, C.: Remobilizing surfactant retarded fluid particle interfaces: II. controlling the surface mobility at interfaces of solutions containing surface active components. J. Coll. Interface Sci. 163(1), 177–189 (1994)
  * Paananen et al. [2020] Paananen, R.O., Viitaja, T., Olzynska, A., Ekholm, F.S., Moilanen, J., Cwiklik, L.: Interactions of polar lipids with cholesteryl ester multilayers elucidate tear film lipid layer structure. Ocul. Surf. 18, 545–553 (2020)
  * Stone [1990] Stone, H.: A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface. Phys. Fluids A 2, 111–112 (1990)