Dynamics and mechanisms for tear breakup (TBU) on the ocular surface

1.   Introduction

The tear film is critical for nourishment and protection for the ocular surface ^[1,2] as well as providing a smooth optical surface essential for good vision. ^[3,4]. A sketch of the lacrimal supply and drainage system is shown in Figure 1. The tear film is painted onto the ocular surface during the upstroke of the upper eyelid during the blink; it is a thin liquid film with multiple layers that, in health, establishes itself rapidly after each blink.

At the anterior interface with air is an oily lipid layer that primarily retards evaporation ^[7,8], which helps to retain a smooth, well-functioning tear film ^[9]. Posterior to the lipid layer is the aqueous layer, which consists mostly of water ^[10]. At the corneal surface, 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 hydrophilic glycocalyx on the ocular surface prevents the tear film from dewetting ^[11,12,13]. The overall thickness of the tear film is a few microns ^[14], while the average thickness of the lipid layer is on the order of tens to 100 nanometers ^{[9,15,16,17,18]} and the thickness of the glycocalyx is a few tenths of a micron ^[2]. This overall structure is reformed on the order of a second after each blink in a healthy tear film.

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. Doane ^[5] proposed the mechanism of tear drainage in vivo whereby tear fluid is drained into the canaliculi through the puncta during the opening interblink phase. Water lost from the tear film due to evaporation into air is an important process also ^[7,19,20]. 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 ^[18]. Some water is supplied from the ocular epithelia via osmosis when the osmolarity increases above the isotonic level (the so-called hyperosmolar state) ^[18,21,22]. 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 ^[23], which affects millions ^[24,25]. An accepted view in these DEWS reports is that repeated exposure to sufficient hyperosmolarity in tears, as may occur in tear breakup, can result in inflammation and ocular surface damage in a vicious cycle that may can lead to dry eye disease. There is still considerable ongoing research to more fully understand the causes and etiology of dry eye disease; some mathematical modeling results are discussed in this paper help to quantify the relationship between tear breakup and hyperosmolarity.

The characteristic ranges of osmolarity on the ocular surface may be summarized as in ^[26] and ^[27]; the following closely follows those sources. The healthy tear film maintains homeostasis with the blood in the range 296-302 mOsM (or, equivalently, mOsm/L or Osm/m$^3$) ^[28,29,30], while healthy blood is in the range 285-295 mOsM ^[31]. In dry eye disease, the lacrimal system may be unable to maintain this homeostasis and osmolarity values in the meniscus rise to 316-360 mOsM ^{[29,32,33,34]}, and even higher values may be attained over the cornea. Using in vivo experiment and sensory feedback, Liu et al. ^[35] estimated peak values of 800-900 mOsM in tear break up (TBU). Similar or higher values were computed from mathematical models of TBU in ^[18,36], or for models of the whole ocular surface ^[27]. The estimates from these models take on added significance since the osmolarity in TBU or in most locations cannot be measured directly; the only available clinical method used in the clinic measures the osmolarity in the inferior meniscus ^[28].

The ongoing supply and drainage of tear fluid affects the distribution and flow of the tear film. A number of methods have been used to visualize and/or measure tear film thickness and flow, including interferometry ^[14,37,38], optical coherence tomography ^[39], fluorescence imaging ^{[40,41,42,43]} and many others. There are different ways that fluorescence may be used to visualize the tear film ^[44]. In the dilute limit, the fluorescein concentration is below the critical concentration, and the intensity of the fluorescence from the tear film is proportional to its thickness. In the concentrated (or self-quenching) limit, the intensity drops as the tear film thins in response to evaporation, and the thickness is roughly proportional to the square root of the intensity for a spatially uniform (flat) tear film ^[45,46]. We use the shorthand FL to indicate fluorescence or fluorescein below.

Recently, the DEWS reports ^[25] and others (summarized in ^[18]) have argued that TBU is driven, at least in many cases, by evaporation. Evaporation causes increased osmolarity in the tear film, leading to irritation, inflammation and damage to the ocular epithelium. Currently, the authors are unaware of any directly measurements of the osmolarity in regions of TBU over the cornea. Theoretical efforts have responded to this situation by creating models that incorporate evaporation and osmosis, as well as other effects, into the mathematical models to better understand the dynamics of the process at a small scale. These models emphasize understanding the overall dynamics of the thickness, flow, and transport of solutes inside the film. These predictions can produce insights beyond what direct measurements can currently provide, and include significant new results about the distribution of osmolarity in TBU. We believe that these results impact the understanding of osmolarity dynamics in the tear film as well as the understanding of FL imaging of TBU dynamics.

A variety of mathematical models have incorporated various important effects of tear film dynamics as has been recently reviewed ^[18,21,47]. The most common assumptions for these models are a Newtonian tear fluid and a flat cornea ^[48,49]. This paper emphasizes local partial differential equation (PDE) models for tear film dynamics in TBU. Beyond previous reviews, we include some recent results using non-Newtonian fluids and non-uniform properties of the ocular surface. We also review new results for fitting local thin film models to in vivo imaging data; these are reviewed in the Results section. Opportunities for improvements are suggested in the Discussion section. Conclusions follow thereafter.

2.   Materials and methods

This paper reviews mathematical models and theories for human tear film dynamics. We typically refer the reader to the original works for derivations and detailed discussions of solution methods. Generally, multiple scale perturbation methods are applied to the equations of fluid dynamics and transport of mass and energy to derive the so-called thin film or lubrication eqautions for the quantities of interest. Reviews of the sizable literature on thin films include ^[50,51,52]. Recent reviews of tear film models and dynamics include ^{[18,21,47,53]}. The resulting systems are typically dissipative, in one- or two-dimensional domains, and may be solved via a variety of numerical methods; a good summary is given in Maki et al. (^[54], p. 364).

Quite a few parameters appear in the equations. Table 1 gives typical dimensional values that are needed in these models (after ^[55]). The thickness of the tear film ($d$) is about 100 times less the length scale along the tear film ($\ell$); this small ratio enables the use of the thin film approximation to derive governing equations for the dependent variables.

3.   Results

Overall flows over a two-dimensional (2D) eye-shaped domain with stationary boundaries have been able to capture thinning near the lid margins due to the presence of a meniscus there, and flow within the meniscus along each lid margin. Tear film models are often formulated on a one-dimensional (1D) domain oriented vertically through the center of the cornea with stationary ends corresponding to the eyelid margins. We refer to models on this kind of domain as 1D models. Surface tension, viscosity, gravity and evaporation are often incorporated into 1D models; wetting forces have been included as well. Due to space considerations we refer the reader to a recent review ^[21].

3.1.   Overall tear flow over the open eye

In this section we briefly discuss two-dimensional (2D) results over the open eye shape. For a stationary boundary (no blinking), thin film fluid (lubrication) models have been developed for the dynamics of the tear film. A series of papers developing these models have used a piecewise quartic polynomial fit to the lid margins in an image of an eye as the boundary. The first papers used time-independent pressure ^[65] and flux ^[54] boundary conditions. The thickness of the tear film and the flows inside it were computed. Thinning occurred near the lid margin, and flow around the menisci was observed. Subsequent work added time-dependent to the boundary conditions ^[66], evaporation and osmolarity transport ^[27], and fluorescein transport and fluorescence ^[67]. Evaporation and osmolarity transport revealed that hyperosmolarity occurred near the eyelid margins (boundaries) in the so-called "black line" where the tear film is thinnest. Adding fluorescein and fluorescence captured elements of a commonly used procedure for imaging the tear film in the clinic, and set the stage for quantitative comparison with these imaging techniques.

3.2.   Blinking

A sample of papers that review blinking includes ^[68,69]. For a recent review of tear film dynamics and blinking observations, together with some mathematical modeling, the reader is referred to ^[18]. Of particular interest with regard to visualizing tear film dynamics are high-resolution imaging of the lipid layer in both space and time during blinks, simultaneous imaging with fluorescence, and a survey of methods for imaging the tear film. The latter work also includes first efforts to closely compare computed and observed FL intensities in TBU.

3.2.1.   2D domains

Recently, an approach that incorporates blinking of a 2D eye-shaped domain has been developed. A model problem treats the eyelid margins as sectors of circles and uses conformal mapping methods to map the eye shape to more convenient domains for computing and moving the upper lid ^[70]. That work solved a model second-order nonlinear diffusion problem on the blinking domain to high accuracy. Combining the moving domain with thin film fluid flow equations is in progress (TA Driscoll and RJ Braun, in preparation). Preliminary results for one model, corresponding to a model for pure water without evaporation, appeared in ^[47].

3.2.2.   1D domains

The situation for tear film dynamics on moving 1D domains representing the opening of the eyelid and the complete blink cycle has been treated in a number of papers. Some models that combine the eyelid upstroke and interblink include ^[71,72,73], and a two layer model with a dynamic lipid layer ^[74]. Several models have included full blink cycles, including fluid dynamics only ^[75]; flow and heat transfer ^[76,77]; and flow, evaporation and osmolarity ^[78]. For a recent summary of these models, the reader is referred to ^[47].

We briefly discuss two papers that included the lid upstroke and interblink but simplified the lipid layer to its polar lipid component by treating it as an insoluble surfactant. In Aydemir et al. ^[71], the boundary conditions on the insoluble surfactant were of no-flux type; the aqueous flux at the ends was specified. This choice fixed the total amount of surfactant on the tear film surface. They were able to analyze some limits, e.g., the early time dynamics, via a perturbation approach. They obtained good agreement between the analysis and the numerical solutions of the full model. Maki et al. ^[79] used Dirchlet boundary conditions on the insoluble surfactant in an effort to study a hypothesized lipid reservoir at the lid margin ^[1]. They also studied different conditions on the aqueous layer. The model could be analyzed in some limits, yielding good agreement between theory and numerical solutions for parts of the dynamics. The agreement with expectations of tear film dynamics appears to be better for the no-flux conditions in most aspects of the results; however, taken together with the results for the surfactant distributions in two-layer models, there are still aspects of lipid layer dynamics that warrant further study.

We note that Zubkov et al. ^[80] studied flow near the lid margin. They compared flows in the meniscus using 2D Navier-Stokes equations and lubrication approximations near the domain boundary corresponding to the lids. The two approaches showed qualitative agreement there, but the lubrication theory was unable to capture some fine details that the Navier-Stokes model exhibited. Thus, models that rely on lubrication theory for fine detail of meniscus flows should be viewed as qualitative or model problems. We also note that no models to our knowledge have taken the fine details of the geometry of the cross section of the lid margin into account.

3.3.   Tear breakup (TBU)

We now turn to solving local problems for the failure of the tear film via TBU. One solves for the thickness, osmolarity and fluorescein concentration that vary with both time and one space dimension. There are no menisci in these local problems; they are applied to a small region of interest within the open eye domain. The shape of TBU regions observed using FL imaging may be complex, and they are dynamic. PDE models of TBU have been solved only for simple shapes so far. Axisymmetric circular shapes, corresponding to "spots, " have been solved in an axisymmetric polar coordinate system. Linear shapes, corresponding to "streaks, " have been solved in Cartesian coordinates. We limit our discussion to these simplified geometries here. Such models can approximate and produce comparable data for simple TBU regions seen in vivo ^[18].

According to driving mechanism, there are at least two types of TBU: evaporative and rapid glob-driven. Figure 2 shows sketches of the different effects that are present and their relative importance in TBU when these respective mechanisms are dominant. The evaporative mechanism is thought to be from defects in the lipid layer; a local thin region of the lipid layer is sufficient to drive faster localized thinning and subsequent TBU ^[82]. The Marangoni effect is caused by differences in lipid layer composition, in particular for the surface active components at the lipid/aqueous layer interface. The Marangoni effect is known to be active just after a blink ^[48,83], where the formation of new area during eye opening causes a lower surfactant concentration near the upper lid. It has also been been hypothesized to be the driving force in some instances of rapid TBU ^[84,85], where uneven lipid spreading after the blink can cause divergent tangential flow, and associated local thinning, of the tear film.

A third mechanism is thought to be important as well: dewetting. The ocular surface may be defective, causing it to be unwettable ^[86,87]; the idea has been discussed in the context of TBU types ^[88] and for dry eye disease diagnosis ^[89]. In this case, a pressure gradient is developed inside the tear film due to van der Waals forces pushing fluids from between the air/tear and tear/cornea interfaces ^[90]. This mechanism typically causes accelerated decrease of thickness as the film becomes thinner (e.g., ^[52]); however, the level of detail seen in carefully controlled physical systems may not be readily achieved in the relatively messy tear film system. This could make observation of accelerated thinning difficult in humans, though animal models may be more amenable to modification for detailed observation of dewetting-driven TBU.

We first consider evaporative TBU, then proceed to rapid TBU.

3.4.   Thinning driven by evaporation

3.4.1.   One-layer models with fixed distributions

The thin film approximation, or lubrication theory results in a simplified system of equations to solve for the tear film dynamics. The following equations apply to axisymmetric spot TBU; the derivation and the Cartesian versions of the equations may be found in ^[91].

One must solve

These equations are solved for the tear film thickness $h(r, t)$, the pressure $p(r, t)$, the osmolarity $c(r, t)$, and the fluorescein concentration $f(r, t)$. The domain is $0 < r < R_O$, where $R_O$ is the outer radius of the domain. The nondimensional parameters are summarized in Table 2, with some typical values for the scalings used in ^[91].

These equations apply to axisymmetric spot TBU; those for streak TBU are similar. The pressure is treated as a dependent variable to be found as part of the numerical treatment. Treating the pressure this way requires only first order derivatives in the boundary conditions, which is favorable for the numerical methods that we employ ^{[36,65,91,92,93]}.

The boundary conditions are symmetry at the origin, homogenous Neumann at $R_O$, and the initial conditions are uniform in space:

The initial FL concentration $f_0$ is relative to the critical FL concentration $f_{cr} = 0.2 \%$.

The dynamics are driven by the evaporative flux $J_e$, which is scaled with $v_{max}$, the maximum volume flux measured from experiment. The evaporative flux is assumed to be composed of a time-independent distribution at the initial time, which is modified by the film pressure. The latter includes capillary and wetting terms and are motivated by one-sided models of evaporation ^[60]. Empirical evidence for a fixed distribution may be found in Figure 6 of ^[82], where a very thin region of lipid forms, then remains roughly stationary, leading to TBU by 20 seconds post blink. For in vivo experiments with controlled humidity and little airflow, measurements yield constant thinning rates over a significant time interval ^[56]. Thus, we expect that our assumed form of evaporation rate applies best for comparison with data from controlled laboratory conditions. Nondimensionally, we have

We show results for $J(r)$, the evaporation rate with a Gaussian distribution inside a disk, namely

The standard deviation $r_w$ will be used to indicate the radius of the evaporation distribution. Here $v_b = v_{min}/v_{max}$ is the nondimensional version of the background evaporation rate $v_{min}$, and the scaling is based the peak evaporation rate $v_{max}$.

The FL intensity $I$ is computed from the thickness $h$ and the fluorescein concentration $f$ via the nondimensional equation ^[44,45]

Here $\phi$ is the nondimensional Napierian extinction coefficient, and $I_0$ is a normalization constant that takes the optical system into account. As mentioned previously, for fixed thickness $h$ the FL intensity decreases linearly at small concentration in the so-called dilute regime. Sufficiently above a critical concentration $f_{cr} = 0.2\%$, the intensity decreases quadratically with increasing $f$, which we call the self-quenching regime ^[45]. For a tear film thinning by evaporation, the situation is slightly more complicated. For a flat tear film, the product $hf$ is constant, and so in the dilute regime, $I$ is constant during thinning, while in the self-quenching regime, the quadratic decrease for increasing $f$ still holds ^[46].

3.4.2.   Numerical methods

The numerical method used to solve this model (and the other results shown herein) is an application of the method of lines that used a Chebyshev spectral discretization in space ^[94]. The resulting differential-algebraic system being solved using $\texttt{ode15s}$ in Matlab (MathWorks, Natick, MA, USA). This solver uses adaptive BDF time stepping with variable size time steps and variable order of 1 to 5; this solver is well suited to these stiff, high-order problems. With sufficient resolution of the spatial discretization, we found that we did not need to use positivity-preserving schemes that have met with success elsewhere, e.g., ^[95]. We note that finding TBU times was enabled by using event detection capability built into the solver.

3.4.3.   Results with fixed distributions

In this case, we show results from ^[91] for the problem with PDEs given by Eqs (3.1)-(3.4) and (3.7) subject to the homogeneous Neumann boundary conditions and uniform initial conditions Eq (3.6). For illustration, the Gaussian distribution has peak thinning rate $v_{max} = 20\; \mu$m/min, background rate $v_{min} = 1\; \mu$m/min, and width $r_w = 1$; this width is the same as the length scale along the film ($0.35$ mm in this case). The time scale is $d/v_{max} = 10.5$ s. This approximates a hole in the lipid layer allowing localized fast evaporation compared to a slow background rate. In this example, wetting forces are active ($A\neq 0$) and there is a pressure correction to the evaporation rate ($\alpha\neq 0$). The thickness is scaled with typical film thickness $d$, the evaporative flux with the maximum evaporation rate $\rho v_{max}$, the osmolarity $c$ is scaled with the isotonic concentration $c_0$, and the fluorescein concentration $f$ with its critical concentration.

Figure 3 shows tear film thickness $h$, evaporative flux $J_e$, osmolarity $c$ and FL concentration $f$ as a function of space for different times. We see that the tear film thickness decreases to a low level that we interpret as TBU and that this takes slightly longer than one unit of time; this confirms the suitability of our choice of time scale for a thinning rate distribution of this width. When the thickness of the tear film gets small enough, the evaporative flux is reduced according to the model. The osmolarity and FL concentrations develop elevated values in the TBU region, with the peak in $f$ larger than the peak in $c$ due to its smaller diffusivity. However, these peak values of osmolarity are not sufficient to stop the thinning, which progresses to TBU ^[36]. The evaporative flux $J_e$ decreases at the thinnest $h$ values due to the presence of the wetting term $Ah^{-3}$. This prevents $h$ from decreasing below about $0.25\ \mu$m to mimic the thickness of the glycocalyx ^[96]. If desired, we could compute solutions beyond the initial appearance of TBU, but we do not treat the spreading of TBU in this model because it is limited by the assumed evaporation distribution. While the TBU time is often used as a measure of tear film instability, the progression of TBU as measured by increasing area of TBU during the interblink is an important clinical aspect ^[97]. We note that when different evaporation models are used, the thinning may stop at much smaller values than we use here ^[36,98]. For a more complete discussion, see ^[91].

In Figure 4, we show TBU time as a function of spot size $r_w$ for different peak thinning rates $v_{max}$. The TBU time is defined as when the tear film thickness reaches twice the wetting equilibrium thickness $h_{eq} = (A\alpha)^{1/3}$ where evaporation stops due to the balance with the wetting term. In Figure 4, we use $h_{eq} = 0.25\; \mu$m, so the TBU time was taken to be when the tear film reached $0.5\; \mu$m. This size was chosen because the mathematical equilibrium is only slowly approached, which artificially extends the TBU time. In these results, osmosis is active with $P_o = 12.1\; \mu$m/s permeability ^[18].

The results show that for relatively large spots with $r_w > 0.8$ and $w_0\ge 10\mu$m/min, the TBU time is about one unit of time. This indicates that the dimensional time to get to TBU is about $d/v_{max}$ which is the time for a flat film to evaporate away. If the peak thinning rate is small, such as $v_{max} = 5\; \mu$m/min, then the thinning is slow and the TBU time increases because there is sufficient time for capillary flow and osmosis to supply fluid to the thinning region ^[91].

In the range of moderate spot sizes $0.4 \le r_w \le 0.8$ and with $v_{max} \ge 10\; \mu$m/min, the TBU time increases as $r_w$ decreases. In this regime, capillary flow is increasing, which brings more fluid into the TBU region. However, osmolarity diffuses out rapidly enough that the peak value of $c$ does not approach the high levels that occur for the flat film case ^[36], and thinning can proceed to TBU in similar times to the large spots.

Finally, in the range of small spots with $r_w\le 0.4$, the TBU time is very long because capillary flow into the TBU region (healing flow due to capillarity) is fast relative to the diffusion of solutes out of the TBU region. Increased advection causes osmolarity to reach values near the flat film case, and then osmosis and capillarity cooperate to slow thinning dramatically. However, the theoretical TBU times are long compared to in vivo results that we have fit with TBU models. This suggests that a different or additional mechanism may be present in small TBU spots that is not currently in the models.

Besides affecting TBU time via flow and osmolarity distributions, these dynamics have important consequences for imaging the tear film from the distributions of $f$ in TBU.

Figure 5 shows results at different values of spot size and initial fluorescein concentration for evaporative TBU. It is clear that the best agreement between thickness and a function of intensity, for this mechanism, is for the lower left part of the array for evaporative TBU. This area is where $f_0$ is at the critical concentration or higher, and the spot sizes are near or larger than the length scale $\ell$ (about 0.35 mm in this case). For further details about variation with these parameters, the interested reader is referred to ^[91].

3.4.4.   Parameter estimation from in vivo FL imaging

In ^[55], the evaporation-driven model above was optimized to match intensity data extracted from in vivo experiments on subjects with a healthy tear film. The instances of TBU selected had simple spot (circular) and streak (linear) shapes for directly matching to the mathematical models. The TBU instances were also selected to be relatively slow to develop, as is expected for evaporation-driven TBU (e.g., ^[82]). The details of the instances used and the method for extracting the intensity data are given in their paper. The theoretical model is solved computationally, the intensity evaluated, and the process repeated until optimal parameters for matching the experimental intensity are found. Matlab's $\texttt{lsqnonlin}$ was used with either Levenburg-Marquardt or trust region reflective options (MathWorks, Natick, MA, USA). Detailed results of the fits can be found in ^[55].

Here we focus on the thinning rate found from the best fits and compare those results with experimentally determined thinning rates. Results are shown in Figure 6.

The fitting to the experimental intensity determines the peak evaporation rates $v_{max}$ and the background evaporation rate $v_{min}$. The resulting computed thinning rate (denoted $dh/dt$) corresponds to the experimentally measured thinning rate from Nichols et al. ^[56]. We note that the experimental data is a point measurement from the central cornea and did not specifically target TBU instances. The results for thinning rate fall within the experimental range, and their maximum evaporation rate can fall above the maximum experimentially-measured thinning rate. Since flow within the film can reduce the thinning rate, particularly for evaporative thinning, this is an expected result.

3.4.5.   Aqueous and mucin layers

Sharma and Ruckenstein developed mathematical models with van der Waals driven rupture using linear ^[99] and nonlinear ^[100,101] theories. They included a separate (fluid) mucus layer between the aqueous layer and the epithelium, and there was an insoluble surfactant on the interface between the two layers. These two-layer models were shear-dominated in both layers, and a surfactant transport governed the dynamics of the surfactant on the interface between them. The expressions for the fluid and surfactant fluxes are complicated; see those papers for details. The mucus layer was unstable to van der Waals forces in the model, and the models could all give reasonable TBU time ranges.

The two-layer film theory was generalized to include van der Waals forces in both mucus and aqueous layers, and surfactant transport on the aqueous-air interface ^[102,103]. The mucus layer was treated as a power-law fluid with a fit to experimental data for whole tears ^[104]. The aqueous layer was assumed to be Newtonian, and the lipid layer was simplified to the transport of an insoluble surfactant. They also found that the tear film could be unstable, with TBU driven by van der Waals forces. Thinner mucus layers led to reduced TBU times, while increased Marangoni effect led to increased TBU times. Related papers, which may apply to eyes via analogy with lung surfactants, include ^[105]. For a comprehensive review of related bilayer work, see ^[52].

Dey et al ^[106,107] posed a model with continuous viscosity variation through the depth of the tear film. It replaced separate aqeuous and mucin layers with variable mucin concentration and concentration-dependent viscosity. The concentration is governed by a transport equation. They computed linear and nonlinear stabilty results for the model, and found that for the fixed average velocity, a mucin distribution that caused higher viscosity near the ocular surface stabilized the tear film by increasing TBU time.

3.4.6.   Aqueous and lipid layers

The TBU model by Peng et al. ^[36] adds a stationary non-polar lipid layer atop a dynamic aqueous layer. The evaporation rate is determined by two factors. It is inversely proportional to the lipid layer thickness and it incorporates the resistance to water vapor transport through the boundary layer in air just outside the tear film. That model was the first to find that the osmolarity rises dramatically in TBU, up to about 10 times the isotonic value. This result supports the notion that hyperosmolarity and evaporative TBU can contribute to ocular surface inflammation and dry eye. This paper also found that the osmolarity diffuses out of the TBU region, and this prevents osmosis from stopping TBU when there is spatial variation in the tear film quantities. The paper also investigated various conditions in the air outside the tear film.

Stapf et al. ^[98] extended ^[36] to include a mobile nonpolar lipid layer and an insoluble surfactant representing polar lipid at the lipid/aqueous interface. It was based on an earlier model of Bruna and Breward ^[74] that was used for one dimensional blinking. Both the nonpolar lipid and aqueous layer were assumed to be Newtonian; the lipid layer is dominated by extension, and the aqueous layer is dominated by shear. A nonlinear least squares fit was carried out for the dependence of the evaporation on the thickness the lipid layer from measured data ^[108]. They found that evaporative TBU can occur within this model with concentrated surfactant; the surfactant spread, causing a local thin region in the lipid layer, which then enhances evaporation there, leading to localized thinning of the tear film. The results presented in ^[98] were a TBU process dominated by evaporation; the early thinning of the lipid layer was very rapid, but did not drive TBU. Beyond that motion driven by the Marangoni effect, the lipid layer dynamics were rather slow compared to the other processes such as evaporation and aqueous flux due to capillarity. This result is consonant with simultaneous imaging experiments ^[82]: those observations contained cases where the lipid layer is essentially quiescent while thinning occurs due to evaporation through holes in it. The fixed lipid layer model of ^[36] was motivated by such observations, but the model with two mobile layers could capture the rapid development of the lipid imperfections that lead to increased evaporation processes that are slower.

3.5.   One-layer models with tangential flow driving thinning

The above two-layer models generate tangential flows that help drive thinning toward TBU. In this section, we turn to single-layer models which may have either an insoluble surfactant or a boundary condition playing the role of the lipid layer.

3.5.1.   Mobile insoluble surfactant driving the flow

A model for TBU with an insoluble surfactant was developed by Siddique et al. (2015). In that model, the surfactant is assumed to be a monolayer. The evaporation rate linearly decreases with the surface concentration of a surfactant that represented the polar lipid layer of tear film. The nondimensional evaporation rate was given by

where $\delta$, $\beta$ and $\bar{K}$ are constants. The results from that model showed initial rapid thinning occurs, but then the tear film could develop a competition between capillarity-driven healing flow and an outward flow from the Marangoni effect. The upshot of the model was that for the parameters tested, TBU occurred in too long a time for comparison with experiments of interest.

Results for a simplified model are now presented in this review. The evaporation rate is assumed independent of surfactant concentration, and entire tear film surface is mobile with an insoluble surfactant. Flow is driven by an initial "pile" of excess surfactant at the center of the film, and the surfactant is governed by a transport equation on the surface of the film. The low surface tension there causes a rapid divergent tangntial flow that thins the tear film in the center of the domain. Typically such thinning is much faster than the evaporative thinning ^[82].

Here, we give the equations in Cartesian coordinates for the streak model. Nondimensionally, the system is given by

here,

In addition to the dependent variables in the evaporative case, we must also find the insoluble surfactant concentration $\Gamma(r, t)$. The surface velocity $u_s$ also appears in the surfactant transport equation. The problem is solved on $0 < x < X_0$, where $X_0$ is the outer end of the domain. Homogeneous Neumann boundary conditions are applied at $x = 0$ and $x = X_0$. The initial conditions are the same as for the spot case with $r$ and $R_0$ replaced by $x$ and $X_0$. We give results here for a uniform evaporation rate (constant in space and time). Nondimensional parameters arising from scaling the problem are given in Table 3.

The initial condition for $\Gamma$ is a high concentration of unity near the origin and a low concentration outside $X_I$, using the transition function $B$. The initial condition uses the transition function $B(x)$:

where

Here, $X_I$ defines the half-width of the region of higher concentration of lipid, and $X_W = 0.1$ is the (fixed) transition width.

Computational solution of this problem yields dynamics where the thickness decreases rapidly to a limited degree due to the outward divergent flow. That flow is driven by surfactant concentration gradients and the Marangoni effect, but the surfactant distribution rapidly relaxes toward a more uniform distribution, and this lessens outward tangential flow as time increases. This slows thinning and results in longer TBU times when compared to experiment or other TBU models to follow. We find that the thickness and intensity lag experimental results later in the trial; this model has difficulty keeping up with experimentally observed thinning. A version of these results for a representative set of dimensional parameter values are shown in Figure 7. In the upper left, we see that the thickness and intensity are quite similar in this case, with the intensity slightly lagging the thickness. The osmolarity and FL concentrations (upper right), are also quite similar. The pressure variation decreases with increasing time (middle left), and the surfactant distribution rapidly flattens (middle right). The depth-averaged and surface velocities (bottom row) are both primarily outward and their magnitude decreases steadily as the surface concentration of surfactant flattens.

Thus in both models described in this subsection, the insoluble monolayer with a linearized equation of state appears to relax too rapidly to drive the tear film to TBU. We now turn to a model which addresses this issue via a hypothesized "glob" of excess lipid that is tangentially immobile on the tear film surface.

3.5.2.   Glob-driven rapid TBU

Zhong et al. ^[84] developed a model that hypothesized that a localized region of excess lipid, a so-called lipid glob, could supply polar lipid (surfactant) for a period of time long enough to drive rapid TBU. Experimental evidence suggests that this occurs in some kinds of rapid TBU ^[82]. Their model has a region under the glob that is deformable, is tangentially immobile and has constant high concentration of surfactant. Outside the glob, the surface is mobile, and a standard surfactant transport equation conserves surfactant. No flux conditions apply at the center of the glob and at the edge of the film away from the glob.

These different conditions on the surface are combined by weighting them with a smooth blend function $B(r)$ (for a spot) that has a narrow transition at the edge of the glob. For the streak geometry, the Cartesian model applies, and the glob edge is at the fixed location $r = R_I$. The nondimensional blend function is given by

One may choose $R_w$ to be small so that the transition width of $B(R; R_I, R_W)$ is very narrow compared to the domain size. $B(r; R_I, R_W)$ then approximates a step function which is zero to the left of $R_I$ but one to the right. The equations on the two subdomains are then combined into one continuous equation over the whole domain.

Once the boundary conditions are blended using this function, then lubrication theory can be applied to the thin film throughout the domain to obtain simplified equations for the dependent variables. The reduced equations in this case govern the tear film thickness $h(r, t)$, pressure $p(r, t)$ the solute concentrations $c(r, t)$ and $f(r, t)$ and surfactant surface concentration $\Gamma(r, t)$:

and

where the osmotic flux is given by

and the depth-averaged velocity component along the film is given by

One also has a surfactant transport equation governing the surface concentration of $\Gamma$,

where the surface velocity component along the film is given by

The boundary conditions are as follow. At the glob center, $r = 0$, we enforce symmetry with homogeneous Neumann conditions to obtain the following equations for use in numerical solution:

At $r = R_O$, no-flux conditions require

Finally, we choose uniform initial conditions except for the surfactant concentration:

We also have the consistent pressure initial condition $p = -A$. For the initial $\Gamma$, $\Gamma = 1$ in the glob and 0.1 outside the glob with the smooth, narrow transition between them. All dynamics are driven by the surfactant distribution and any nonzero evaporative flux $J_e$. The variable step size and variable order solver ($\texttt{ode15 s}$) that we used for the spatially-discretized system can readily handle the computations on short time intervals that we find for this model.

Different choices for the evaporative flux $J_e(x)$ were made in Zhong et al. ^[84,85]. We only give a sample of those results here using a constant evaporation rate; the reader is referred to the original papers for more complete results.

Figure 8 portrays a spot TBU result for our default parameters from Tables 1 and with no evaporation ($J_e = 0$). In this case, TBU occurs and is induced only by the different composition of lipid inside the glob. The locally-elevated surfactant concentration there spreads out immediately over the tear/air interface, and the surfactant concentration variation generates a significant shear stress, which in turn drives a strong divergent tangential flow. A significant depression of the tear/air interface forms at 0.15 s. In 0.59 s, the aqueous layer thins from 3.5 to 0.25 $\mu$m, the thickness which we defined as TBU. This TBU time is near unity, which indicates that the time scales we chose work well. The Marangoni flow dominates the thinning process in this case, which leads to TBU in less than a second. In some in vivo observations with simultaneous fluorescence and lipid layer interferometry ^[82], the dark spot has already appeared by 0.14 s (practically instantaneously); the bright circular glob in the lipid layer spreads rapidly but its center remains relatively bright and thus relatively thick. This experimental observation corresponds well with the model.

Within an appropriate range of glob size, $R_I$, TBU can happen in about a second or less, successfully capturing the short time scale in rapid TBU. Figure 8 shows that $h$ thins to 0.25 $\mu$m at $0.57$ s, which is the threshold for TBU we defined for this model. At the same time, $\Gamma$ spreads outside the glob (to the right), so the glob acts as a source that drives tangential flow away from the glob. Due to the short time and our assumed uniform evaporation profile, $c$ and $f$ increased only slightly under the glob. The solute concentration increases were quite small compared to the evaporative TBU cases discussed above.

Comparing the first and last columns of Figure 8 shows that $I$ is a good approximation to $h$ in this case. This is because the intensity change is all from the thickness change since $f$ changes very little. Detailed parametric study showed that the best results for imaging the tear film with rapid TBU appeared to be around $f_0 = 0.7$ (or 0.14$\%$) for this fast thinning case. This value of FL concentration is below the critical concentration ^[85,109]. For related cases, a value slightly lower could also work best. Thus, for rapid thinning, using a more dilute initial fluorescein concentration, together with the intensity itself as indicator of thickness, appeared to betst approximate thickness of the tear film. See ^[85] for a more detailed explanation.

Equations (3.25) and (3.27) imply that the Marangoni contribution to flow is directly proportional to $\partial_r \Gamma$. As the surfactant spreads out from the glob, $\partial_r \Gamma$ on the tear/air interface decreases, which results in a weaker Marangoni flow. Meanwhile, the contribution to aqueous flow from capillarity increases due to the deformation of the tear film. In some cases, the aqueous layer thins around the edge of the glob at short times, but then thickens at longer times. When that happens, the divergent flow from the Marangoni effect is insufficient to lead to TBU because the capillary contribution to flow dominates the dynamics at later times. This phenomenon matches our expectations for the model and some experimental observations where TBU does not occur.

3.5.3.   Mixed-mechanism models and parameter estimation

Recently, Luke et al. ^[57] chose to explore parameters and scalings that intentionally combined both evaporation and tangential flow at comparable levels. The models are rescaled versions of the models developed by Zhong et al. ^[85] but then wrapped in optimization to determine parameter values that best match experimental observations in a least squares sense. The models were solved computationally, and the computed FL intensity was calculated from the resulting $h$ and $f$ distributions as functions of space and time. The 2-norm of the difference between computed ($I_{th}$) and experimentally measured ($I_{ex}$) FL intensities was minimized over the glob size ($R_I$), Marangoni effect ($\Delta \sigma_0$), and evaporation parameters ($J_e$). In terms of a continuous optimization, they solved

where the argmin designates the minimization over the parameters. Eight instances of TBU were fit this way, and tangential flow was found to contribute in all cases they selected. In a small number of cases, evaporation only played a small role, while there was a fairly even mix in five cases. We consider this to be strong evidence that divergent tangential flow is important in TBU instances that are significantly under 10 s for TBU time.

The value of the osmolarity in TBU is of particular interest to clinicians due to its potential for causing pain and inflammation ^{[110,111,112]}. Results are shown for both evaporative and mixed-mechanism TBU in Figure 9.

This figure shows that rapid TBU can still reach significant osmolarity values even when evaporation is moderate. Those results argue against eliminating evaporation from TBU modeling except for the extremely rapid cases (under 1 s for TBU time). The figure also shows that the osmolarity in the fits does not continue to grow with time as the trial duration grows. The evaporative cases peaked in the 650-900 mOsM range, with a typical value around 700 mOsM. These results suggest that the high osmolarities obtained from ODE models for spatially uniform tear films ^[18,47] should be viewed as upper bounds that are not very tight when compared to PDE models. We also note that these fits are from a limited sample, and so the range in vivo may be larger than we have found here.

4.   Discussion

The tear film is a complex fluid structure that is highly dynamic ^[46]. There are recent studies that solve model problems for flows all over the exposed ocular surface. The problems use thin film approximations, even in the meniscus; they have also used simplified material properties. The models still capture a number of aspects of observed tear film dynamics in vivo including flow around the eyelid margins. However, the models do not yet closely link fine scale dynamics for phenomena such as tear break up to the overall flow or blinking dynamics. There is clearly room for progress in this area.

Blinking dynamics and tear film dynamics have been studied more extensively in one-dimensional domains. The situation is similar in that a number of aspects of the overall flows have been described and explained by the models. However, linking small scale dynamics of TBU, for example, have not been linked to the overall flow. And simplified material parameters have been used in these models as well, but more complex models have been studied that included two-layer films. There is room for including more realistic parameters here also.

TBU has been studied in local thin film models. The models have differed in complexity, but they have had some successes. Evaporative models have predicted high levels of osmolarity ^[36,91] that comport well with clinical estimates based on subject feedback ^[35]. More recently, parameters of TBU have been estimated by fitting theoretical models to in vivo experiments on normal subjest using FL imaging ^[55,57]. The thinning rates found by parameter estimation have agreed well with in vivo measurements ^[56] that did not specifically target TBU. Additional insights were gained about thinning rates versus evaporation rates. The relationship between the of TBU and osmolarity are complex, however. Very fast instances of TBU may have a small increase in osmolarity, and more slowly developing instances caused by evaporation may have quite high osmolarity. But there are intermediate cases where both outward tangential flow and evaporation contribute to TBU, and the resulting hyperosmolarity takes on an intermediate value. Developing fast, automated methods to fit the TBU data may lead to inferring the hyperosmolarity in TBU and to quantitative understanding of the hyperosmolarity in regions of TBU. These kinds of models offer the possibility of inferring the mechanism in specific instances of TBU as well; we believe that this close connection between theory and experiment is crucial to better understanding of the tear film in general, and of TBU and dry eye disease in particular.

To date, it has been easier to image the aqueous and lipid layers in the tear film ^[81]. This has contributed substantially to the modeling efforts described in this work. Imaging the glycocalyx and the ocular surface remains more challenging than the other parts of the tear film ^[2,113]. It may be that advances in technology for that purpose lead to advances in modeling and pinpointing TBU thought to be driven primarily by dewetting or other effects chiefly due to the epithelial surface ^[114]. These remain areas for advancement.

5.   Conclusions

While mathematical modeling of tear film dynamics has made progress in explaining and predicting some aspects of the tear film, there remain many open problems and challenges. Better models in both 1D and 2D domains with blinking would be welcome, as well as the use of realistic material properties. We also believe that it is crucial to work closely with clinicians to better understand the tear film. Nothing is simple about the tear film, yet there is a need to simplify the mathematical descriptions of the tear film to make progress. Making appropriate modeling decisions based on insight from data and from clinicians is essential. Testing the results on data leads to better understanding for all involved.

Acknowledgments

This work was supported by National Science Foundation grant DMS 1909846. The content is solely the responsibility of the authors and does not necessarily represent the official views of the funding source. The authors gratefully acknowledge P.E. King-Smith for helpful conversations and sharing of data.

Conflict of interest

The authors declare that they have no conflicts of interest.