Glaucoma refers to a group of eye diseases in which the optic nerve degenerates. Due to the degeneration of optic nerve fibers, glaucoma causes progressive, irreversible death of retinal ganglion cells, eventually leading to blindness. As the second leading cause of blindness, a number of parameters have been studied as possible causal factors in glaucoma [1]. Parameters include elevated intraocular pressure (IOP), low ocular perfusion pressure, increased scleral elasticity, age, ethnicity, myopia, local vascular abnormalities, and alterations in biomechanical properties of the optic nerve head (ONH) among other alternatives [2]. However, the mechanism by which retina ganglion cells are damaged in glaucoma remains controversial presumably because individual patients demonstrate a wide range of sensitivities to these risk factors. For example, although elevated IOP has been considered as the primary risk factor for the development of glaucoma, it still cannot be used as a reliable indicator either of the glaucomatous status or likelihood of progressive ONH changes. Both ocular hypertensives (people with high IOP in the absence of glaucoma) and normal or low tension glaucoma (glaucoma with normal or low IOP) are very common.

Evidence increasingly suggests that abnormal biomechanical properties of the ONH may play an important role in the development of glaucoma [3]. Ganglion cell axons form into bundles and pass through pores in the lamina cribrosa of the ONH before exiting the eye. Because the ONH is a discontinuous region (weak point) in the corneoscleral shell [3], we postulated that pulse-induced changes in IOP might lead to pulsatile deformation. In the presence of abnormal ONH biomechanical properties, stress and strain resulting from pulse-induced motion could directly damage the retinal ganglion cells [4], disturb the capillary circulation perfusing the ONH, or obstruct nutrient transport to [5] or cause chronic progressive deformation of ONH structures.

Characterization of pulse-induced axial ONH movement in vivo would provide a valuable tool to evaluate ONH biomechanical properties because such properties determine ONH responses to IOP forces impinging on it. The extent to which ONH mechanical properties determine susceptibility to damage from IOP is unknown because functional measurement tools have been lacking. Pulse-induced movement of the ONH offers one such tool, but no currently available technology is capable of measuring ONH movement, probably because it is too small (typically a few micrometers).

Since reported in 1991 [6], OCT has become a routine diagnostic tool in ophthalmology [7]. One recent report indicated that structural OCT may be used to monitor pulse-induced movement of the ONH by measuring the change in distance between the cornea and the ONH using OCT images acquired by a newly developed Fourier domain dual-OCT system [8]. However, the axial resolution of the reported system was relatively low (7 μm), while the measured maximal ONH movement was only $\sim 10\mathrm{\mu m}$ for normal subjects and $\sim 14\mathrm{\mu m}$ for glaucoma subjects. It is difficult for the technique to distinguish axial motion of the entire globe from axial motion specific to the ONH. In addition, the use of the pixel-based structural OCT images also imposes limits on measurement accuracy.

For accurate evaluation of pulse-induced ONH movement, we postulate that phase-sensitive OCT (PhS-OCT) using time-lapse B-scans could measure micrometer-scale movement of fundus tissue and isolate the ONH component. PhS-OCT was first reported to measure tissue Doppler [9] using the phase information between adjacent A-scans. It was quickly realized that the evaluation of phases between adjacent A-scans is not amenable to measurement of the small tissue movement at a level of less than one wavelength (i.e., micrometer scale). To solve this problem, PhS-OCT was further developed by evaluating phases between adjacent B-scans. The dramatic increase of the time interval between the B-scans results in sensitivity enhancement to the nanometer scale [10]. However, measurement of the small pulsatile ONH movements using PhS-OCT presents imposing problems; inevitable gross patient movement causing a bulk tissue-motion imaging artifact is the primary challenge. In this Letter, we propose a novel algorithm to meet this challenge by using phase information from retina tissue near the ONH as a reference to compensate for the bulk-tissue motion artifact; then phase changes due to ONH pulsatile movement are extracted.

The spectral domain (SD)-OCT system used here is similar to the one previously reported [11]. Briefly, the system is equipped with a dual high-speed spectrometer detection unit, capable of providing 500 kHz A-scan line rate [11] with an axial resolution of $\sim 7\mathrm{\mu m}$ in air and 95 dB system sensitivity at 0.5 mm imaging depth. The system sensitivity to the tissue movement was evaluated at 0.3 nm, sufficient to measure the small movement of the ONH. Before imaging ONH pulsatile movements, we first produced a fundus image using a traditional $4\mathrm{mm}\times 4\mathrm{mm}$ 3D scan of the ONH region. The result is shown in Fig. 1(a). Scanning positions are marked, indicating locations used to measure pulsatile ONH movements and blood flow in the central retinal artery (CRA) (see below).

 

Fig. 1. (a) OCT fundus image of ONH, (b) typical cross-sectional image, and (c) corresponding phase difference map. R, retina; C, choroid; PL, pre-lamina; LC, laminar cribrosa. $\text{Bar}=500\mathrm{\mu m}$.

Download Full Size | PDF

To measure the ONH movements, a protocol using repeated B-scans (i.e., time-lapsed mode) was designed. 600 A-lines were captured to form one B-scan (covering $\sim 3\mathrm{mm}$ in length). The scanning position is illustrated in Fig. 1(a) (red line with the arrow indicating the scanning direction). Figure 1(b) is a typical cross-sectional image. To assure success of ONH motion extraction, the primary requirement of the OCT scanning region should include both the ONH and peripapillary retina. In the experiments, the imaging rate was 500 frames per second. At this rate, the maximum detectable velocity without phase wrapping is $\sim 105\mathrm{\mu m}/\mathrm{s}$. Considering that the human heart rate fundamental frequency is $\sim 1\mathrm{Hz}$, 2600 B-frames were acquired for each dataset within $\sim 5.2\mathrm{s}$ to provide coverage of $\sim 5$ heart cycles.

The phase difference map between adjacent B-frames was calculated, as shown in Fig. 1(c). Within the tissue region, the map shows obvious tissue motion ($\mathrm{\Delta }\phi$), generated by both the bulk tissue motion ($\mathrm{\Delta }{\phi }_m$) and the ONH movement ($\mathrm{\Delta }{\phi }_0$). Here, $\mathrm{\Delta }\phi =\mathrm{\Delta }{\phi }_m+\mathrm{\Delta }{\phi }_0$. To obtain $\mathrm{\Delta }\phi$ from the phase-difference map, a histogram method reported in [12] was performed on each A-line scan. According to the relationship between $\mathrm{\Delta }\phi$ and the tissue velocity, $\mathrm{\Delta }\phi$ can be presented as $\mathrm{\Delta }\phi =4n\pi v\mathrm{\Delta }t/\lambda +\mathrm{\Delta }{\phi }_0$, where $n$ is the refractive index of the sample, $v$ is the velocity of the tissue motion, $\lambda$ is the central wavelength of the OCT system (842 nm) and $\mathrm{\Delta }t$ is the time interval between adjacent B-frames (2 ms). For in vivo imaging of a tissue that is in constant motion (including bulk and localized motion), the tissue velocity $v$ is a function of time. Because the subject’s head was secured on a slit-lamp headrest during imaging, the bulk tissue motion should be small. Under this condition, the bulk tissue velocity should be a slowly varying function of time. Within a short time period of 2 ms, the acceleration of the bulk tissue movement can be considered constant, meaning that it can be described by a first-order polynomial function, $v=at+v_i$, where $a$ is the acceleration and $v_i$ is the initial velocity at the beginning of the B-scan. Accordingly, $\mathrm{\Delta }\phi$ due to the tissue motion is also a function of time and can be represented as $\mathrm{\Delta }\phi =4n\pi at\mathrm{\Delta }t/\lambda +\mathrm{\Delta }{\phi }_0+{\phi }_i$, where ${\phi }_i$ is due to the initial tissue velocity at the beginning of the B-scan, which will not affect the calculation of ONH movement. A typical $\mathrm{\Delta }\phi$ evaluated from adjacent B-scans is shown in Fig. 2(a), marked as the blue curve. Because of the inevitable bulk motion of the eye, the evaluated $\mathrm{\Delta }\phi$ is continuously decreasing within one B-scan, demonstrating that the acceleration of the bulk movements is almost constant.

 

Fig. 2. (a) Typical curve of the raw $\mathrm{\Delta }\phi$ between adjacent B-scans and (b) corresponding tissue motion after compensating the bulk tissue-motion. (c) Typical curve of the raw $\mathrm{\Delta }\phi$ indicating the problem of phase wrapping (the inset shows $\mathrm{\Delta }\phi$ map across the B-scan) and (d) corresponding $\mathrm{\Delta }\phi$ curve after phase unwrapping. The spikes were due to the blood flow.

Download Full Size | PDF

In order to extract $\mathrm{\Delta }{\phi }_0$ due to the ONH movement, the $\mathrm{\Delta }\phi$ curve was partitioned into two regions [in Fig. 2(a)]: peripapillary retina (the yellow square) and ONH (the green square), which is done with the help of the OCT structural image [e.g., Fig. 1(b)]. We consider the motion within the retinal region to result from the bulk tissue motion; thus, the motion occurring in the retinal region can be used as a reference to extract the motion in the ONH region. Note that the acceleration of the tissue bulk movement is assumed to be constant between adjacent B-scans. Note also that $\mathrm{\Delta }{\phi }_0\sim =0$ within the retinal region. We determined that acceleration of the bulk tissue movement could be evaluated through directly fitting a first-order polynomial function to the $\mathrm{\Delta }\phi$ values in the retinal region [the red line in Fig. 2(a)]. Therefore, with known scan timing, the bulk tissue movement within the region of ONH can be extrapolated from this fitted polynomial function, shown as the broken line in Fig. 2(a). Consequently, the $\mathrm{\Delta }{\phi }_0$ due to the ONH movement can be extracted by a simple subtraction of the extrapolated bulk tissue movement from the $\mathrm{\Delta }\phi$ within the ONH. The result is given in Fig. 2(b). The resultant values within the reference region (the yellow square) remain close to zero, whereas, within the ONH region, they deviate from the reference, demonstrating that ONH movement separates from bulk tissue motion.

However, if the tissue motion is relatively large, the evaluated $\mathrm{\Delta }\phi$ would be phase-wrapped because of the $2\pi$-modulo in the sinusoidal function. Figure 2(c) gives an example of typical $\mathrm{\Delta }\phi$ results when the phase wrapping occurs. The blue curve is the original $\mathrm{\Delta }\phi$ evaluated from each A-line. Because of the phase wrapping, there is an abrupt step jump. This step jump has to be corrected before proceeding to extract the ONH movement. To unwrap the $\mathrm{\Delta }\phi$ curve, we applied a phase-unwrapping algorithm, leading to the result shown in Fig. 2(d). Compared to Fig. 2(a), Fig. 2(d) presents a similar $\mathrm{\Delta }\phi$ curve after the unwrapping procedure, with which the algorithms discussed above can then be used to extract the desired $\mathrm{\Delta }{\phi }_0$ due to the ONH movement.

After all the $\mathrm{\Delta }\phi$ curves were evaluated from 2600 B-frames using the algorithms described above, the final results were stacked to produce a 2D motion map. The result is illustrated in Fig. 3(a), in which the right part demonstrates obvious oscillatory patterns caused by the pulse-induced ONH movements. Figure 3(b) plots the velocity curve of pulse-induced ONH movements extracted from the position marked by the black line in Fig. 3(a), corresponding to the location denoted by “*” in Fig. 1(a). For better visualization, the y axis of Fig. 3(b) is reversed, demonstrating a magnitude from $-30$ to $30\mathrm{\mu m}/\mathrm{s}$. Figure 3(c) is the map indicating displacement of the ONH tissue obtained through integrating Fig. 3(a) over the time lapse $t$. Figure 3(d) plots the displacement curve of the ONH extracted from the same position as in Fig. 3(b). Following the method reported in [13], we performed Fourier analysis of Figs. 3(b) and 3(d), giving a fundamental frequency of 1.2 Hz. The second- and third-order harmonics are also presented.

 

Fig. 3. (a) Velocity map of ONH movements and (b) a typical ONH velocity curve at the position marked by the black line in (a). (c) Displacement map of ONH and (d) displacement curve corresponding to (c). (e) Typical Doppler flow map and (f) dynamic blood flow measured from CRA. $\text{Bar}=500\mathrm{\mu m}$.

Download Full Size | PDF

Axial movements of the ONH as measured above were correlated with pulsatile blood flow in CRA using repeated B-scans ($\sim 5.2\mathrm{s}$ duration) to measure the dynamic flow velocity in the CRA. To form one B-scan covering $\sim 3\mathrm{mm}$ across the ONH that contained the CRA, 2500 A-lines were captured [the scanning location is marked by the green line in Fig. 1(a)]. Figure 3(e) shows a typical blood flow map for one B-scan obtained by calculating the phase differences between adjacent A-lines, i.e., a phase-resolved algorithm. To evaluate the pulsatile flow, the phase difference values of the CRA [blue circle in Fig. 3(e)] were integrated over time. The result is given in Fig. 3(f), where, as expected, blood flow within the CRA is clearly pulsatile. Comparison of the charts of Figs. 3(b) and 3(f) reveals that the two curves are similar; both have a steep increase at the beginning of each cycle followed by a slow decay after the cycle peak. The Fourier frequency analyses determined the fundamental frequency to be 1.2 Hz, equivalent to the heartbeat of the subject. CRA pulse frequency correlates 100% with the fundamental frequency found in the ONH tissue movement, a consistent finding in five normal subjects.

In conclusion, we demonstrate that our high-speed SD-OCT system quantitatively measures pulsatile axial movements of the ONH in human subjects. We have shown that the pulsatile ONH movement is correlated with the pulsatile blood flow in the CRA. Pulse-induced ONH movement is determined by biomechanical properties subject to characterization by this PhS-OCT approach. Measurement of these properties may be clinically valuable in both predicting ONH susceptibility and in monitoring changes in such properties that may occur as ONH damage progresses.