Quantitative evaluation of biomechanical properties of optic nerve head by using acoustic radiation force optical coherence elastography

Abstract. Significance Previous studies have demonstrated that the biomechanical properties of the optic nerve head (ONH) are associated with a variety of ophthalmic diseases; however, they have not been adequately studied. Aim We aimed to obtain a two-dimensional (2D) velocity distribution image based on the one-to-one correspondence between velocity values and position using the acoustic radiation force optical coherence elastography (ARF-OCE) technique combined with a 2D phase velocity algorithm. Approach An ARF-OCE system has the advantages of non-invasive detection, high resolution, high sensitivity, and high-speed imaging for quantifying the biomechanical properties of the ONH at different intraocular pressures (IOPs) and detection directions. The 2D phase velocity algorithm is used to calculate the phase velocity values at each position within the imaging region, and then the 2D velocity distribution image is realized by mapping the velocity values to the corresponding structure based on the one-to-one relationship between velocity and position. The elasticity changes can be read directly according to the quantitative relationship between Lamb wave velocity and Young’s modulus. Results Our quantitative results show that the phase velocity and Young’s modulus of the ONH increase by 32.50% and 129.44%, respectively, with increasing IOP, which is in general agreement with the results of previous studies, but they did not produce large fluctuations with the constant change of the ONH direction. These results are consistent with the changes of elastic information in the 2D velocity distribution image. Conclusions The results suggest that the ARF-OCE technology has great potential in detecting the biomechanical properties of the ONH at different IOPs and directions, and thus may offer the possibility of clinical applications.


Introduction
Glaucoma is an optic neuropathy, which is one of the leading causes of blindness. 1,2It is estimated that by 2040, ∼112 million people worldwide will be affected by glaucoma, 11 million of whom will be blind due to glaucoma. 3Studies have shown that a range of ophthalmic disorders encompassed by glaucoma are associated with abnormal changes in intraocular pressure (IOP), but the mechanism of interaction of elevated IOP in the development of optic neuropathy is unclear. 4,5ustained elevation of IOP will cause structural changes in the posterior segment of the eye, specifically leading to outward expansion of the sclera around the optic nerve head (ONH), forcing the sieve plate and adjacent structures to undergo mechanical stresses associated with IOP, resulting in compression and remodeling. 6Therefore, a better understanding of the biomechanical properties of the ONH is of great clinical value to better understand how to effectively prevent and control the development of glaucomatous disease. 7urrently, many advanced methods are used to directly observe the ONH deformations caused by IOP changes, such as three-dimensional (3D) histomorphometry, 8 confocal microscopy, 9 microcomputed tomography, 10 ultrasound images, 11,12 magnetic resonance images, 13 and optical coherence tomography (OCT). 14,15However, these methods have some limitations because the ONH and the lamina cribrosa (LC) tissues can not only displace posteriorly but also anteriorly under the influence of IOP.Therefore, it has also been observed that the LC tissues do not undergo significant displacement and deformation in a subset of patients. 16Therefore, there is an urgent need for an elasticity method to quantify the biomechanical properties of the ONH tissue.
Elastography is used as a non-destructive imaging tool to quantify the biomechanical properties of ocular tissues.Among the commonly used elastography techniques, magnetic resonance elastography 17 and ultrasound elastography 18 are the most widely used ones, which is mainly because of their high imaging depth and wide field of view, but their relatively low spatial resolution limits the detection of small changes of the velocity.Although Brillouin microscopy has higher spatial resolution and has been used to study the biomechanical properties of ocular tissues, the correlation between Brillouin frequency shift and Young's modulus remains unclear due to the uncertainty of Poisson's ratio, limiting its clinical application. 19Optical coherence elastography is an emerging elastography technique that quantifies the biomechanical properties of ocular tissues by detecting the propagation of elastic waves induced by external excitation using the OCT technique. 20,21][24] In this paper, acoustic radiation force optical coherence elastography (ARF-OCE) technology is used to detect the ONH elasticity, which has been well-validated in our previous study for ocular tissue elasticity imaging. 25,26In this study, we first investigated the biomechanical properties of the ONH at different IOPs, meanwhile the corresponding information of the ONH at different directions have also been study since it may have an influence on the ONH elasticity, which has not been researched before.In the data processing, the phase velocity algorithm is used to calculate the phase velocity at each position of the imaging region, meanwhile the 2D velocity distribution image is realized by mapping the velocity values to the structure based on the one-toone relationship between the velocity value and the position.Therefore, changes of the velocity and elasticity information can be directly read from the structure map based on the Lamb wave velocity-Young's modulus relationship.Furthermore, the average Young's modulus values of the samples are acquired.

Materials
To conduct ONH elasticity experiments, fresh porcine eyes (n ¼ 4) from a standard slaughterhouse are obtained, which are kept in cold saline for transfer to the laboratory.Moreover, the experiments are completed within 24 h after removal.Before the experiment, the excess tissue around the eye is removed by a scalpel to keep the ONH flush with the outer surface of the surrounding sclera.Subsequently, the entire eyes are fixed on a custom-built stent so that the cornea are facing down and the ONH are facing up.In addition, the eyeball is placed in sterile phosphate-buffered saline (PBS) to maintain the normal physiological state of the eye tissue, which also serves as a transport medium for ultrasonic wave.

ARF-OCE System Design
A customized ARF-OCE system is designed to assess the biomechanical properties of the ONH.The ARF-OCE system consists of a phase-sensitive OCT system and an acoustic radiation force excitation system, as shown in Fig. 1.The OCT system uses a swept-source laser (Axsun Technologies Inc., Billerica, Massachusetts) with a central wavelength of 1300 nm and a scan repetition rate of 50 kHz.The axial and lateral resolution are measured to be 6.7 and 20 μm, respectively.In order to stabilize the phase in OCT data, a fiber Bragg grating was added to produce wavelength swept trigger for the data acquisition card to record the interference signal.The optical phase stability of ARF-OCE system is calculated as ∼167.3mrad.The light from the laser is split by the optical fiber coupler (90:10), with 90% of the output light into the sample arm and 10% output light into the reference arm.The back-reflected and back-scattered light from both arms is collected to generate an interference signal by the 50:50 coupler, which is then detected by a balanced photodetector.The signal of the detector is digitized by using a waveform digitizer acquisition card (Alazar Technologies Inc., Quebec, Canada) and processed by the computer.
The ultrasound excitation system is composed of an ultrasound transducer, a function generator, and a power amplifier.A custom-built ultrasound transducer with a center frequency of 4.5 MHz, −6 dB lateral acoustic beam width of 390 μm, and a focal length of 35 mm is used to excite tissues, which is placed at a position orthogonal to the OCT beam and can cover the OCT imaging depth in the ONH (∼0.4 mm).The trigger signal of the swept-source laser is utilized to synchronize the function generator to produce a 4.5 MHz sine wave with a duration of 500 μs, which is amplified by ∼46 dB to drive the ultrasound excitation to produce an elastic wave in the tissue.
The IOP control system mainly consists of a micro syringe pump, a catheter, a 23G needle, and a digital pressure sensor, which connects the syringe filled with electrolyte complex intraocular flush solution, the pressure sensor, and the eye through a medical T-connector.The different IOPs are controlled by adjusting the composite electrolyte intraocular rinse solution fed into the eye by the micro-infusion pump and monitored in real-time using the digital pressure sensor.

Data Acquisition
The M-B protocol is used to visualize the propagation of elastic waves in the ONH tissue, as shown in Fig. 2, which is performed in the transverse direction from P 1 to P n .At each lateral position, 1000 A-lines are acquired to record the intensity and phase change over time (M-mode), where the ARF excitation is generated between the 101st and 125th A-lines.After completing Fig. 1 The schematic of the ARF-OCE system.
one M-mode data acquisition, the galvanometer scanner moves the detection beam to an adjacent position and repeats the same operation.

Results
First, we used a Lamb wave model to quantify the biomechanical properties of the ex vivo ONH at normal IOP (10 mmHg), as shown in Fig. 3, where the areas to the left and right of the yellow dashed line in Fig. 3(a) shows of the ONH and sclera, respectively.Figures 3(b)-3(f) show the propagation process of elastic waves at five different times.Furthermore, the calculation process of the phase velocity at a specific depth is shown in Fig. 4. The spatial-temporal maps are processed with the two-dimensional Fourier transform (2D-FFT) to obtain the k-space distribution of the frequency with respect to the wave number.According to Eq. ( 4), the phase velocity can be calculated by the frequency at each point divided by the maximum intensity value at that frequency.Due to the viscous characteristic of the ONH, the phase velocity profile shows a diffuse nature, as shown in Fig. 4(d).The phase velocity values in the dispersion curve are then extracted based on the main frequency corresponding to each position. 26Accordingly, the phase velocity values of each position in the imaging region are achieved, which are then mapped to the corresponding structure to generate the 2D velocity distribution image, as shown in Fig. 5(a), where the velocity changes can be directly observed.In addition, the change of the biomechanical properties can also be read directly from the distribution image based on the relationship between velocity and Young's modulus, as shown in Eq. (6).
To investigate the effect of IOP changes on the biomechanical properties of the ONH, we synchronized and integrated a customized IOP control system on the ARF-OCE system, in which the IOP was controlled at 10, 15, 20, 25, 30, and 35 mmHg, respectively.Figures 5(a 1.In order to visualize the changing trend, the average phase velocity and Young's modulus of the ONH at different IOPs were plotted in Fig. 6, which was fitted with a second-order polynomial function with a goodness-of-fit R-squared of 99.6% due to the wide  range of tissue deformation and nonlinear effects.(Each phase velocity value, modulus value, and error bar in the figure denotes the mean and standard deviation of three repeated measurements, respectively.)It can be found from Figs. 5 and 6 that the phase velocity and Young's modulus increase by 32.50% and 129.44%, respectively, with the increase of IOP, which is in general agreement with the results of previous studies. 27These results are consistent with the changes in Fig. 5.
Subsequently, to better investigate the biomechanical properties of the ONH in different directions, the eye was rotated at 45 deg intervals along the center of the ONH and marked with sterile marker lines to facilitate imaging of the ONH by the system.Figures 7(a 2. Figure 8 plots the average phase velocity and Young's modulus obtained at different directions.The comparison of the study results shows that the phase velocity and Young's modulus did not produce large fluctuations with the constant change of the ONH direction.

Discussion
In previous research work, several groups have investigated the biomechanical properties of the ONH by tensile test, 28 AFM, 29 and ultrasound elastography. 30However, tensile test measurements require cutting the ONH or peripapillary sclera tissue into strips meanwhile AFM measurements need micron-sized specimens. 31,32Therefore, the test conditions of these different  research methods lead to a wide range of deviations of the obtained Young's modulus from the in vivo measurements.To further quantify the biomechanical properties of the ONH and to perform non-invasive imaging of the ONH tissue.Ma et al. used a 3D cross-correlation scatter tracking algorithm to measure the tissue deformation in the posterior segment of the eye, but this method only provided stress-strain results and did not allow for quantitative measurements of Young's modulus. 33Therefore, Zhou et al. developed the confocal ARF-OCE method to obtain the average elastic modulus of the ONH. 34Although it has demonstrated that the ARF-OCE technique   can obtain the average velocity of elastic waves, the wave velocity of each position cannot be achieved and thus its 2D distribution image cannot be realized.In this work, we combined the ARF-OCE system with the phase velocity algorithm to generate 2D phase velocity distribution images of the ONH at different IOPs and detection directions.Our method not only can obtain the average Young's modulus but also can offer the possibility to read the elasticity changes directly from the structure map.Therefore, this study further demonstrated the superiority of the ARF-OCE system in detecting the biomechanical properties of the ONH.
Although the technique has been used to quantitatively analyze the biomechanical properties of the ONH with elevated IOP and altered detection direction, there are still some challenges and limitations that need to be addressed for the successful translation of the technology for clinical applications.First, faster lasers are required to reduce the impact of artifacts caused by body shaking, breathing, and heartbeat on the imaging quality. 35In addition, due to the faster propagation of elastic waves in stiffer objects, the wave propagation may not be tracked if the acquisition speed of the imaging system is not fast enough. 36Second, PBS was used as the ultrasonic propagation medium; however, it will inevitably lead to a bad user experience and thus an air-coupled ultrasound transducer should be employed to better fulfil requirements of in vivo measurements.Meanwhile, developing a novel OCE system based on a broadband laser with a central wavelength of 850 nm and a zoom ring transducer will be used to detect ONH in the retina through the lens.In addition, the mechanical index (MI) of the ARF applied in our study should be well below the FDA ophthalmic MI standard of 0.23. 37Third, the effective depth resolution of the system can be improved by increasing the sampling points in the depth direction.Finally, a detailed analysis of wave propagation models in bounded media is essential to assess tissue biomechanical properties.In future studies, more sophisticated and advanced wave propagation models will be developed to quantify the viscoelasticity of tissues by evaluating the dispersion of elastic wave propagation.

Conclusion
In this paper, we have demonstrated that the ARF-OCE technology can be used to quantify the biomechanical properties of the ONH at different IOPs and detection directions.More importantly, the 2D velocity distribution images can be realized based on the one-to-one correspondence between the velocity value and the position, where changes of the velocity and Young's modulus can be directly read.Moreover, our quantitative result shows that the Young's modulus of the ONH increases significantly with increasing IOP, but it basically does not change with the change of the ONH detection direction.As a result, this technology has great potential in detecting the biomechanical properties of the ONH at different IOPs and directions, and thus may offer the possibility of clinical applications.
6 Appendix: Data Processing

Displacement Map
In our experiments, the phase-resolved Doppler algorithm was adopted to extract the phase information that varies with time.With the M-mode OCT image, Doppler phase shift profile of the adjacent A-lines can be calculated using the following equation: 38 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 4 ; 2 1 1 where F m is the complex signal at the given position of the M-mode, F mþ1 is F m at the next time point, and F Ã is the conjugate complex of F. Im () and Re () are the imaginary and real parts of the OCT complex signal, respectively.Thus, according to phase change Δφ, the axial displacement Δd can be calculated by using the following equation: 27 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 1 2 5 where λ 0 is the center wavelength of the laser source, and n is the refractive index.

Velocity and Modulus Calculation
Since ONH is a thin plate structure at the tissue structure level and is small compared to the transverse wave wavelength.Lamb waves usually propagate in a medium with upper and lower boundaries, thus the lamb wave is identified as the elastic waves propagating in the ONH.The lamb wave can only propagate in the low frequency range due to the viscosity in soft tissues.Moreover, the zeroth order mode can travel at any frequency and is most commonly detected in the low-frequency range as compared with non-zeroth order mode Lamb waves.For a thin plate immersed in an incompressible fluid, the dispersion equation based on the zero-order antisymmetric Lamb wave mode is where , k L ¼ ω∕c L is the Lamb wavenumber, ω is the angular frequency, c L is the frequency dependent phase velocity, k S ¼ ω ffiffiffiffiffiffiffiffiffi ρ∕U p is the shear wavenumber, U ¼ μ þ iωη is the shear modulus, μ and η are shear elasticity and shear viscosity, respectively, and h is the half thickness of the ONH.
The spatial-temporal Doppler displacement map of the entire ONH tissue is first obtained from the intensity and phase of the raw data.Due to the presence of the low-frequency noise, the spatial-temporal Doppler displacement map is filtered by a high-pass filtering algorithm and then converted to the wavenumber/frequency domain (k-space map) by 2D-FFT.Based on the wavenumber-frequency domain map, the phase velocity of the Lamb wave is obtained by dividing the frequency by the corresponding maximum wavenumber as follows: 40 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 7 ; 4 7 9 where ω is the angular frequency, and κ L is the Lamb wave number.Since the thickness of the ONH tissue is much smaller than the wavelength of the elastic wave generated by the excitation, the Lamb wave velocity V L in the vacuum can be expressed as 41 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 7 ; 4 0 8 where f is the frequency, h is the ONH thickness, and V s is the shear wave velocity.Since the ONH tissue is submerged by liquid, the corresponding Lamb wave velocity can be corrected by multiplying a factor of 1∕ ffiffi ffi 2 p the total longitudinal wave leakage and the total transverse wave reflection at both boundaries (V L ¼ V L_vacuum ∕ ffiffi ffi 2 p ).According to E ¼ 3ρV 2 , the corresponding Young's modulus can be calculated based on the Lamb wave velocity by the following equation: 42 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 3 0 0 where ρ is the density of the ONH, V L is the Lamb wave velocity.

Disclosure
The authors declare no competing financial interests.

Fig. 2
Fig. 2 The M-B scanning protocol for the OCE system.(a) The OCT trigger signal for each M-scan image, (b) the trigger signal of the ultrasonic transducer, (c) the signals for controlling the X-axis galvanometer scanner, and (d) the position of the galvanometer in relation to the ARF.
)-7(d) show the 2D phase velocity distribution obtained at different directions (color scale indicates the phase velocity values in m/s).For changing the detection direction, the mean Young's modulus were calculated to be 177.29 AE 7.77 kPa at 0 deg, 183.11AE 15.91 kPa at 45 deg, 182.14 AE 15.20 kPa at 90 deg, and 180.24AE 22.46 kPa at 135 deg, which are collected in Table

Fig. 6
Fig. 6 The biomechanical properties of the ONH at different IOPs.(a) Averaged phase velocity and (b) average Young's modulus.

Fig. 8
Fig. 8 The biomechanical properties of the ONH at different directions.(a) Averaged phase velocity and (b) average Young's modulus.

Table 1
The result of ex vivo ONH.
a (M ± SD) represents mean and standard deviation.

Table 2
The result of ex vivo ONH.