Parameter Estimation for Evaporation-Driven Tear Film Thinning

Abstract

Many parameters affect tear film thickness and fluorescent intensity distributions over time; exact values or ranges for some are not well known. We conduct parameter estimation by fitting to fluorescent intensity data recorded from normal subjects’ tear films. The fitting is done with thin film fluid dynamics models that are nonlinear partial differential equation models for the thickness, osmolarity and fluorescein concentration of the tear film for circular (spot) or linear (streak) tear film breakup. The corresponding fluorescent intensity is computed from the tear film thickness and fluorescein concentration. The least squares error between computed and experimental fluorescent intensity determines the parameters. The results vary across subjects and trials. The optimal values for variables that cannot be measured in vivo within tear film breakup often fall within accepted experimental ranges for related tear film dynamics; however, some instances suggest that a wider range of parameter values may be acceptable.

Appendix

1.1 Governing Dimensional Equations

For the circular case, we use the dimensional axisymmetric coordinates $(r',z')$ to denote the position and $\varvec{u}' = (u',w')$ to denote the fluid velocity. The tear film is modeled as an incompressible Newtonian fluid on $0< r' < R_0$ and $0< z' < h'(r',t')$, where $h'(r',t')$ denotes the thickness of the film. Conservation of mass and momentum of the TF fluid and transport of solutes within the fluid are given, respectively, by

$$\begin{aligned}&\nabla ' \cdot \varvec{u}' = 0, \end{aligned}$$

(19)

$$\begin{aligned}&\varvec{u}'_{t'} + (\varvec{u}' \cdot \nabla ') \varvec{u}' = - \frac{1}{\rho } \nabla ' p' + \nu \nabla '^2 \varvec{u}', \end{aligned}$$

(20)

$$\begin{aligned}&\partial _{t'} s' + \nabla ' \cdot ( \varvec{u}' s') = D_s \nabla '^2 s', \end{aligned}$$

(21)

where $p'$ is the fluid pressure and $s'$ represents either $c'$, the osmolarity, or $f'$, the fluorescein concentration, with diffusivities $D_o$ and $D_f$, respectively. The fluid density is $\rho$ and the kinematic viscosity is $\nu$.

At the film/cornea interface $z' = 0$, we require no slip and osmosis across a perfect semipermeable membrane:

$$\begin{aligned} u' = 0, \quad w' = P_o V_w (c'-c_0). \end{aligned}$$

(22)

The membrane permeability is given by $P_o$, the molar volume of water is $V_w$, and $c_0$ is the isotonic osmolarity.

We enforce no flux of solutes across both the film/cornea and film/air interfaces:

$$\begin{aligned} D_o \varvec{n}' \cdot \nabla ' c' - \varvec{n}' \cdot \varvec{u}'c' = 0, \quad D_f \varvec{n}' \cdot \nabla ' f' - \varvec{n}' \cdot \varvec{u}' f' = 0. \end{aligned}$$

(23)

At $z' = 0'$, the outward normal is given by $\varvec{n} = (0,-1)$, and at $z = h'(r',t')$,

$$\begin{aligned} \varvec{n} = \displaystyle \frac{(-\nabla '_{II} h',1)}{(1 + |\nabla '_{II}h'|^2)^{1/2}}, \end{aligned}$$

(24)

where $\nabla '_{II}$ is the gradient in the plane of the substrate parallel to $z' = 0$.

The kinematic condition implies that the balance of the material derivative of the TF thickness and the fluid velocity in the $z'-$ direction is controlled by the evaporative mass flux, $J'$:

$$\begin{aligned} \rho \left[ \frac{ \partial _{t'} h' + u' \nabla '_{II} h' - w'}{(1 + |\nabla '_{II} h'|^2)^{1/2}} \right] = -J', \end{aligned}$$

(25)

where

$$\begin{aligned} J' = \rho \left[ v_{\min } + (v_{\max } - v_{\min }) e^{-(r'/r_w)^2/2} \right] + \alpha _0 (p'-p_v'). \end{aligned}$$

(26)

Here, $\alpha _0$ is effectively $\alpha /K$ from Ajaev and Homsy (2001), $v_{\min }$ and $v_{\max }$ are background and peak thinning rates, respectively, $r_w$ is the standard deviation that corresponds to the width of the evaporation distribution, and $p_v'$ is atmospheric pressure. The contribution of pressure to evaporation may be ignored with little consequence.

The normal stress condition at $z' = h'(r',t')$ is given by

$$\begin{aligned} p'-p'_v = - \left( \sigma _0 \nabla '_s \cdot \varvec{n}' + \frac{A^*}{h'^3}\right) , \end{aligned}$$

(27)

where $\sigma _0$ is the surface tension, $\nabla '_s = (I - \varvec{n}' \varvec{n}') \cdot \nabla$ (Stone 1990), and $A^*$ is the Hamaker constant.

1.2 Derivation of Tear Film Equations, Spot Case

Using the scalings (2), (3), we nondimensionalize the governing equations as in Braun et al. (2018). At leading order, conservation of mass and momentum of the fluid on $0< z < h(r,t)$ are given by

$$\begin{aligned}&\frac{1}{r} \partial _r(ru) + \partial _z w = 0, \end{aligned}$$

(28)

$$\begin{aligned}&\partial _z^2 u = \partial _r p, \quad \partial _z p = 0. \end{aligned}$$

(29)

The leading order boundary conditions at $z = 0$ are

$$\begin{aligned} u = 0, \quad w = P_c(c-1). \end{aligned}$$

(30)

The leading order boundary conditions at $z = h(r,t)$ are

$$\begin{aligned}&\partial _t h + u \partial _r h - w = -J, \end{aligned}$$

(31)

$$\begin{aligned}&p = - \frac{1}{r} \partial _r (r \partial _r h) - \frac{A}{h^3}. \end{aligned}$$

(32)

Integrating (28) over the vertical domain, applying the Leibniz rule, and using (31) and (30) to substitute in for the first three resulting terms gives

$$\begin{aligned} \partial _t h + J - P_c(c-1) + \frac{1}{r } \partial _r (r h {\bar{u}}) = 0, \end{aligned}$$

(33)

where

$$\begin{aligned} \displaystyle {\bar{u}}(r,t) = \frac{1}{h} \int _0^h u(r,z,t) dz \end{aligned}$$

(34)

is the depth-averaged fluid velocity. The radial velocity u in the tangentially immobile case is given by

$$\begin{aligned} u(z) = \frac{1}{2}( z^2 - zh) \partial _r p. \end{aligned}$$

(35)

For solutes, we keep all powers of $\epsilon$ before assuming an expansion in this small parameter. The nondimensional solute transport equations are

$$\begin{aligned} \partial _t s + u \partial _r s + w \partial _z s = \text {Pe}_s^{-1} \left[ \frac{1}{r} \partial _r (r \partial _r s) + \epsilon ^{-2} \partial _z^2 s \right] , \end{aligned}$$

(36)

where s denotes either osmolarity c or fluorescein concentration f, and $\hbox {Pe}_s = \displaystyle \frac{v_{\max } \ell }{\epsilon D_s}$ is the Péclet number $\hbox {Pe}_c$ or $\hbox {Pe}_f$, respectively. We continue the derivation for the osmolarity c following Jensen and Grotberg (1993).

The solute boundary condition at $z = 0$ is

$$\begin{aligned} \epsilon ^{-2} \text {Pe}_c^{-1} \partial _z c = wc, \end{aligned}$$

(37)

and the boundary condition at $z = h(r,t)$ is

$$\begin{aligned} \partial _z c = \text {Pe}_c\epsilon ( u \partial _r h - w) c + \epsilon ^2 \partial _r h \partial _r c . \end{aligned}$$

(38)

Assume that c(r, z, t) can be expanded as:

$$\begin{aligned} c = c_0(r,z,t) + \epsilon ^2 c_1(r,z,t) + O(\epsilon ^4). \end{aligned}$$

(39)

After substituting this expression for c into (36), the leading order equation is given by

$$\begin{aligned} \partial _z^2 c_0 = 0, \end{aligned}$$

(40)

and thus $c_0 = c_0(r,t)$. The next order in $\epsilon$ results in

$$\begin{aligned} \partial _z^2 c_1 = \text {Pe}_c (\partial _t c_0 + u \partial _r c_0) - \frac{1}{r} \partial _r (r \partial _r c_0). \end{aligned}$$

(41)

Integrating (41) over the vertical domain gives

$$\begin{aligned} \partial _z c_1(r,h,t) - \partial _z c_1(r,0,t) = \text {Pe}_c (h \partial _t c_0 + h {\bar{u}} \partial _r c_0) - \frac{1}{r} h \partial _r (r \partial _r c_0). \end{aligned}$$

(42)

The terms involving $c_1$ can be eliminated by identifying the boundary conditions at $O(\epsilon ^2)$; these result in an equation for $c_0$. We drop the subscript to give our leading order PDE for osmolarity:

$$\begin{aligned} h(\partial _t c + {\bar{u}} \partial _r c) = \frac{1}{r}\text {Pe}_c^{-1} \partial _r (r \partial _r c) + Jc - P_c(c-1)c. \end{aligned}$$

(43)

The evolution equation for f may be obtained similarly:

$$\begin{aligned} h(\partial _t f + {\bar{u}} \partial _r f) = \frac{1}{r}\text {Pe}_f^{-1} \partial _r (r \partial _r f) + Jf - P_c(c-1)f. \end{aligned}$$

(44)

1.3 Tear Film Equations, Streak Case

The derivation of the problem in the linear case for streaks is similar to the axisymmetric case, and more details may be found in Braun et al. (2015, 2018). The nondimensionalization is the same in both cases.

The problem is solved on $0< x < X_L$ and $0< z < h(x,t)$, where h is the TF thickness. Homogeneous Neumann boundary conditions are applied at $x = 0$ and $x = X_L$. The fluid velocity coordinates in the (x, z) directions are given by (u, w). Nondimensionally, the system is given by

$$\begin{aligned}&\partial _t h - J - P_c(c-1) + \partial _x(h {\bar{u}}) = 0, \end{aligned}$$

(45)

$$\begin{aligned}&{\bar{u}} = - \frac{h^2}{12} \partial _x p, \end{aligned}$$

(46)

$$\begin{aligned}&p = - \partial _x^2 h - A h^{-3}, \end{aligned}$$

(47)

$$\begin{aligned}&h(\partial _t c + {\bar{u}} \partial _x c) = \text {Pe}_c^{-1} \partial _x (h \partial _x c) - P_c(c-1) c + Jc, \end{aligned}$$

(48)

$$\begin{aligned}&h(\partial _t f + {\bar{u}} \partial _x f) = \text {Pe}_f^{-1} \partial _x (h \partial _x f) - P_c(c-1) f + Jf. \end{aligned}$$

(49)

Nondimensionally, the evaporation distribution is given by

$$\begin{aligned} J(x,t) = v_b + (1-v_b) e^{-(x/x_w)^2/2} + \alpha p. \end{aligned}$$

(50)

The parameters $v_b$ and $\alpha$ are identical to that in the spot case, and r and $r_w$ have simply been replaced with x and $x_w$.

1.4 TBU Images

The last image from each trial is shown in Fig. 18 with fitted TBU highlighted. These correspond to the results given in Tables 3 and 4.

Fig. 18

The last image in each trial: a S1v2t7, b S1v2t8, c S1v2t10, d S8v2t3, e S27v2t9, f S28v1t3 (contrast of f enhanced 2$\times$). The bright rectangle, called the Purkinje, is due to the reflection from the light source (Color figure online)