Understanding tear film dynamics is a prerequisite for advancing the management of Dry Eye Disease (DED). In this paper, we discuss the use of optical coherence tomography (OCT) and statistical decision theory to analyze the tear film dynamics of a digital phantom. We implement a maximum-likelihood (ML) estimator to interpret OCT data based on mathematical models of Fourier-Domain OCT and the tear film. With the methodology of task-based assessment, we quantify the tradeoffs among key imaging system parameters. We find, on the assumption that the broadband light source is characterized by circular Gaussian statistics, ML estimates of 40 nm +/− 4 nm for an axial resolution of 1 μm and an integration time of 5 μs. Finally, the estimator is validated with a digital phantom of tear film dynamics, which reveals estimates of nanometer precision.
© 2013 OSA
Dry Eye Disease (DED) is a multifactorial disease of the tears and the ocular surface. Symptoms include discomfort, visual disturbance, and tear film instability, and in some cases damage to the ocular surface . DED, a serious public health problem that affects more than 100 million people worldwide, degrades the quality of life seriously [2,3]. However, therapies for DED remain elusive for two reasons: there is no quantitative, highly repeatable, objective diagnosis method, and there is little or no correlation between signs and symptoms . In this work, we hypothesize that tear film instability is a key indicator of DED and can be quantified by tear film dynamics, specifically the thinning rate of the tear film across the cornea.
Optical coherence tomography (OCT) is a biomedical imaging method that has emerged in the last two decades . In OCT, one collects the back-reflected, or scattered, light from tissue. The method is similar to ultrasound, but it operates in the near infrared region of the spectrum. OCT has been used extensively in ophthalmology, especially for imaging the retina. In this work, we propose to use OCT to monitor tear film dynamics. Considering the speed of changes in tear film thicknesses, Fourier-Domain OCT (FD-OCT) is tailored to this challenge.
Statistical decision theory has been widely used in various biomedical fields in which random processes associated with light emission and detection must be accounted for. Statistical decision theory was introduced into time domain OCT for classification tasks [6,7]. We recently reported on key results of using statistical decision theory to estimate tear film thickness . In this paper, we present the detailed mathematical framework that underlies these results, and we further report on the limit of the estimation. To our knowledge, this work represents the first time a maximum-likelihood (ML) estimator has been used to interpret OCT data.
In section 2 of this paper, we present a mathematical model of an FD-OCT system in the context of the ML estimation. We present a model of the tear film in section 3. Next, in section 4, we detail the generation of simulated spectra that serve as inputs to the ML estimator. The principle of the ML estimator is then detailed in section 5, while its performance is investigated in section 6. In the last section of the paper, the ML estimator is validated using a digital phantom.
2. Mathematical model of FD-OCT in the context of the ML estimation
In this work, we consider a free-space setup of an FD-OCT system that is shown schematically in Fig. 1. It is based on a Michelson-like interferometer geometry, in which the detector is a spectrometer.
The broadband light source emits an electric field that can be regarded as a superposition of plane waves and expressed as
The electric field on the xth pixel of the detector can be written asEq. (3) may be rewritten as
The number of electrons, Ng(x,Δt), accumulated at the xth pixel during one detector integration interval Δt can be written asEq. (7) may be rewritten as6] thatEq. (11) may be rewritten asEq. (10), we can express the ensemble average of Np(x,Δt) asEq. (15) can be expressed as Eq. (11), applying the complex Gaussian moment theorem  to the first term on the right-hand side of Eq. (18) yieldsEq. (18), we obtainEq. (5) and Eq. (12), thatEq. (20) we can rewrite it asEq. (17) and Eq. (23), we can express the covariance KNp(x, x’;Δt) as
3. Tear film model
The normal tear film consists of three parts: the lipid layer, the aqueous layer, and the mucin layer. Because the aqueous layer has the largest volume in the tear film, here we only consider the aqueous layer ; we take this as a first approximation, to be refined in future investigations. We model the tear film as a one layer structure with two interfaces, the interface of the air and the tear film and the interface of the tear film and the cornea; these are shown schematically in Fig. 2, where n1, n2, and n3 are the refractive indexes of the air, tear film, and cornea, respectively. In the figure, l is the thickness of the tear film, and r1 and r2 are the reflection coefficients of the two interfaces.
The interface of the tear film and the cornea is considered rough (and is illustrated as such), as it is known that the corneal surface has microplicae and glycocalyx . In our model, we use the root mean squared deviation in the height of the surface, σ, to characterize the surface roughness. We assume that deviations of the corneal surface from a smooth plane can be described by a stochastic process with a Gaussian distribution. We assume that the incident light strikes the corneal surface at a right angle. Using the perturbation theory developed in , the following equation has been established 
4. Generate simulated OCT imaging data sets
To generate one simulated image for a tear film with thickness l, we assume that one instance of the output at the xth pixel follows a normal distribution:
5. Principle of the ML estimator
On the assumption that all pixels are independent of each other, the probability that a particular spectrum Nlg is generated by a tear film of thickness d can be expressed as
The ML estimator estimates parameters by maximizing P(Nlg|d) with respect to d, which is equivalent to minimizing the following expression
Figure 3 illustrates the principle of our ML estimation. We consider spectra that are Gaussian-shaped, with 1 μm, 2 μm, and 4 μm axial resolutions, respectively. Axial resolution, in this paper, is conventionally defined as the half width half maximum of the axial point spread function [7,14]. The source power is set at 1 mW to be in compliance with the ANSI standard for eye safety in this spectral region . The integration time is 5 μs. Ground truth of the tear film thickness is set to be 0.5 μm, and we assume a root mean squared height of the corneal surface of 150 nm  to characterize the corneal roughness.
The plots in red frames of Fig. 3 show the negative log of the conditional probability, which characterizes the probability that one observed spectrum Nlg is generated by a tear film thickness d. Simulations for three different optical axial resolutions are shown. Results show that the conditional probability oscillates when the possible thickness changes. The distance between two adjacent peaks is half of the center wavelength in the sample, or 300 nm (in the tear film). For a tear film whose thickness is 0.5 μm the three resolutions of 1, 2, and 4 μm yield estimates of 505 nm, 485 nm, and 545 nm, respectively, where each estimate is determined by the position of the minimum point in the corresponding conditional probability plot shown in Fig. 3. We will characterize precision of this measurement in section 6.1. Note that with this method we can obtain estimates with precision greater than the optical axial resolution of the system. It is clear from Fig. 3 that for low resolution (i.e., 4 μm) the estimator misses the minimum value because nearby thicknesses yield values comparable to the minimum.
6. Investigation of the performance of the ML estimator
6.1 Impact of detector integration time
Imaging speed is a key factor in the performance of the system. The two parameters that may affect the imaging speed, the repetition rate of the pulsed source and the detector integration time, must be considered. In OCT, a first requirement is that the source speed is at least as great as the detector’s response rate. Therefore we will study the impact of the detector integration time. To quantify how the system performance changes with integration time, we set up a simulation in which we consider a fixed tear film thickness of 0.5 μm and we repeat the estimation 5000 times. We then compute the root mean squared error (RMSE) to evaluate the performance for different optical axial resolutions as a function of integration time. We plot the resulting RMSE as a function of integration time in Fig. 4.
These results show that for a fixed axial resolution, the RMSE decreases exponentially with increasing integration time, eventually reaching a plateau value in the nanometer range. For the same integration time, higher axial resolution results in a higher precision. There are two costs to consider in pushing the system to highest axial resolution: one is the cost of the light source, and the other is the possible reduction in imaging depth. Indeed, as the axial resolution increases, the spectral bandwidth increases as well, but may still be collected on an array of the same size. Therefore, the spectral resolution typically decreases, unless a more costly detector is used; if one uses the same detector, the depth of imaging is typically reduced.
6.2 Limit of the estimation
In estimating thickness, it is important to identify the limit of precision of the ML estimation. To study this limit, we vary the thickness of the ground truth in the simulation and calculate the relative error, which is defined as the ratio of the RMSE to the tear film thickness. In this simulation, the integration time is fixed to be 5 μs and the RMSE is calculated based on 5000 simulated data sets. Results are shown in Fig. 5. For axial resolutions of 1 μm and 4 μm, the relative error diminishes as the tear film thickness increases. For a maximum 10% relative error (shown as a red dashed line in Fig. 5), the ML estimator can estimate thicknesses as low as 0.04 μm and 1.3 μm for axial resolutions of 1 μm and 4 μm, respectively, where each estimate corresponds to the intersection of the red dashed line and the curve of relative error.
7. Validation of the ML estimator with a digital phantom
Maki et al. were the first to develop a comprehensive mathematical model for the tear film dynamics on an eye-shaped domain. They studied the relaxation of the human tear film after a blink on a stationary eye-shaped domain (corresponding to a fully open eye) using a mathematical model based on lubrication theory . The model provides the ground truth of our simulations. We evaluate the performance of the ML estimator assuming an optical axial resolution of 1 μm and integration times of 5 μs and 20 μs, respectively, which is plotted in Fig. 6. Results show that the ML estimator is unbiased, as expected. We further learn that as the integration time increases by a factor of 4, the precision improves by only a factor of 2. Results show uncertainty in the estimates that are less than 10 nm and 5 nm for integration times of 5 μs and 20 μs, respectively.
We have developed and presented a mathematical model of FD-OCT with a single layer (aqueous) tear film model laid on a rough corneal surface. To benchmark the framework of ML estimation in FD-OCT, we implemented an ML estimator to interpret simulated OCT data sets that were generated from the mathematical models. We quantified the tradeoff between the detector integration time and the precision of the estimation. Results show that a combined axial resolution of 1 μm (defined conventionally as the half width at half maximum of the axial point spread function) and an integration time of 5 μs allow estimates in the nanometer range (e.g. 40 nm). Finally, the estimator was validated with a digital phantom of tear film dynamics. Results reveal estimates of nanometer precision.
In future studies we will first further develop the theoretical framework presented in this paper to account for the two-layer structure of the tear film for both its lipid and the aqueous/mucin layers. We will also next consider the statistics of real sources that will be used in our future experiments and compare theoretical results against the performance benchmarked in this paper. In validating the theoretical framework in experiments, physical phantoms will also need to be developed with high accuracy. The theoretical foundation presented in this paper using a digital phantom is guiding step by step the development of an optimal OCT system to quantify tear film dynamics. Key theoretical results of an optimized OCT system developed from this guiding theoretical framework will be fully validated in experiments.
This research was supported by the NYSTAR Foundation C050070, NIH grants RC1-EB010974 and R37-EB000803, and by Research to Prevent Blindness.
References and links
1. M. A. Lemp, C. Baudouin, and J. Baum et al.., “The definition and classification of dry eye disease: report of the definition and classification subcommittee of the International Dry Eye WorkShop (2007),” Ocul. Surf. 5(2), 75–92 (2007). [CrossRef] [PubMed]
4. D. A. Sullivan, K. M. Hammitt, D. A. Schaumberg, B. D. Sullivan, C. G. Begley, P. Gjorstrup, J. S. Garrigue, M. Nakamura, Y. Quentric, S. Barabino, M. Dalton, and G. D. Novack, “Report of the TFOS/ARVO symposium on global treatments for dry eye disease: an unmet need,” Ocul. Surf. 10(2), 108–116 (2012). [CrossRef] [PubMed]
5. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef] [PubMed]
6. J. Rolland, J. O’Daniel, C. Akcay, T. DeLemos, K. S. Lee, K. I. Cheong, E. Clarkson, R. Chakrabarti, and R. Ferris, “Task-based optimization and performance assessment in optical coherence imaging,” J. Opt. Soc. Am. A 22(6), 1132–1142 (2005). [CrossRef] [PubMed]
7. A. C. Akcay, E. Clarkson, and J. P. Rolland, “Effect of source spectral shape on task-based assessment of detection and resolution in optical coherence tomography,” Appl. Opt. 44(35), 7573–7580 (2005). [CrossRef] [PubMed]
8. J. Huang, K. S. Lee, E. Clarkson, M. Kupinski, K. L. Maki, D. S. Ross, J. V. Aquavella, and J. P. Rolland, “Phantom study of tear film dynamics with optical coherence tomography and maximum-likelihood estimation,” Opt. Lett. 38(10), 1721–1723 (2013). [CrossRef] [PubMed]
9. H. H. Barrett and K. J. Myers, Foundations of Image Science (Wiley, Hoboken, 2004), Chap. 14.
10. H. H. Barrett and K. J. Myers, Foundations of Image Science (Wiley, Hoboken, 2004), Chap. 8.
12. T. A. Leskova, A. A. Maradudin, and I. V. Novikov, “Scattering of light from the random interface between two dielectric media with low contrast,” J. Opt. Soc. Am. A 17(7), 1288–1300 (2000). [CrossRef] [PubMed]
13. Y. Ashtamker, V. Freilikher, and J. C. Dainty, “Ambiguity of optical coherence tomography measurements due to rough surface scattering,” Opt. Express 19(22), 21658–21664 (2011). [CrossRef] [PubMed]
15. American National Standards Institute, “American national standard for safe use of lasers,” ANSI Z136.1–2007 (2007).
16. H. B. Chen, S. Yamabayashi, B. Ou, Y. Tanaka, S. Ohno, and S. Tsukahara, “Structure and composition of rat precorneal tear film. A study by an in vivo cryofixation,” Invest. Ophthalmol. Vis. Sci. 38(2), 381–387 (1997). [PubMed]
17. K. L. Maki, R. J. Braun, P. Ucciferro, W. D. Henshaw, and P. E. King-Smith, “Tear film dynamics on an eye-shaped domain II: flux boundary conditions,” J. Fluid Mech. 647, 361–389 (2010). [CrossRef]