Open Access
Add to BibSonomy
Abstract
Acoustic wave velocity measurement based on optical coherence tomography (OCT) is a promising approach to assess the mechanical properties of biological tissues and soft materials. While studies to date have demonstrated proof of concept of different ways to excite and detect mechanical waves, the quantitative performance of this modality as mechanical measurement has been underdeveloped. Here, we investigate the frequency dependent measurement of the wave propagation in viscoelastic tissues, using a piezoelectric point-contact probe driven with various waveforms. We found that a frequency range of 2–10 kHz is a good window for corneal elastography, in which the lowest-order flexural waves can be identified in post processing. We tested our system on tissue-simulating phantoms and ex vivo porcine eyes, and demonstrate reproducibility and inter-sample variability. Using the Kelvin-Voigt model of viscoelasticity, we extracted the shear-elastic modulus and viscosity of the cornea and their correlation with the corneal thickness, curvature, and eyeball mass. Our results show that our method can be a quantitative, useful tool for the mechanical analysis of the cornea.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Changes in mechanical properties of corneal tissues have been linked to structural and geometrical changes clinically observed in corneal degenerative diseases, such as keratoconus [1–3]. Measurement of corneal biomechanics offers opportunities for improving diagnosis and treatment of the diseases. Biomechanical characterization also has the potential to improve the accuracy of routine screening intraocular pressure (IOP) measurements by reducing errors due to variations in corneal mechanical properties [4]. Furthermore, quantitative measurement of corneal tissue stiffness could potentially improve refractive surgeries [5–7].
However, it remains challenging to measure the mechanical properties of corneal tissues quantitatively. Mechanical measurement techniques, such as strip extensometry and eye inflation tests, have been used extensively in laboratory settings, but are destructive and can not easily be applied in vivo [1, 8–10]. Analyzing corneal responses to air puff tonometers, such as the ocular response analyzer (Reichert) and Corvis ST (Oculus), can probe biomechanical changes in normal and pathologic eyes of live patients [11–16], but does not provide a direct, quantitative readout of elastic modulus [17]. Ultrasound elastography has been proposed to map the stiffness of the cornea, but it has relatively low spatial resolution [18, 19]. Brillouin microscopy can map the mechanical properties of tissues with high spatial resolution [20–22], but it measures longitudinal modulus rather than shear modulus that is more directly related to tissue deformation and stiffness.
Optical coherence elastography (OCE) is an emerging technique for localized measurements of shear elastic modulus [23–25]. Traveling-wave OCE is one embodiment of OCE, in which the propagation of mechanical waves in tissues is visualized by OCT. This is particularly well suited for the cornea with well-defined geometry and structure [26–30]. There is a direct theoretical relationship between complex shear wave velocity cs and shear modulus μ:
where ρ is material density. Like other soft materials, the corneal tissue is viscoelastic, described by complex-valued, frequency-dependent shear modulus. As a result, cs has frequency dependence (dispersion), and measurement of cs at a single frequency is insufficient to describe the viscoelastic properties. Furthermore, the speed of waves propagating in the cornea is influenced by the interfaces of the cornea with the air at the anterior surface and with the aqueous humor at posterior surface. The mechanical waves traveling along the cornea are described as guided (Lamb) waves [31]. OCE described in this paper measures the phase velocity of the guided waves c. To compute c and ultimately complex-valued μ from c, it is essential to consider the geometrical effect or waveguide dispersion of the waves, together with material dispersion.A few recent studies have reported measurements of wave dispersion in the cornea using OCE [29–31]. However, there are noticeable discrepancies among the results between different groups and different approaches. The focus of our work described here is to characterize the frequency dependence of wave velocity and its impact on elastic modulus estimation. The following sections describe experimental and analytical methods to generate, measure and analyze the mechanical wave propagation in the cornea, and to measure wave velocity over a frequency range from 0 to 10 kHz.
2. OCE system
Elastography measurements were performed using a home-built, swept-source OCT system previously described [32]. Briefly, this system has an A-line rate of 45 kHz, axial resolution of 15 μm and transverse resolution of 30 μm (full width at half maximum in air) using a polygon swept laser with a tuning range of 80 nm and a center wavelength of 1280 nm. The acquisition of OCT interference fringes is phase-stabilized using the reference signal of an external Mach-Zehnder interferometer [33]. The optical beam is scanned using a two-axis galvanometer scanner (Cambridge Technology, 6210H). The noise-limited sensitivity of the system to vibration is approximately 4 nm [32].
Mechanical stimulation was achieved using a home-built piezoelectric probe (Fig. 1(a)). It consists of a hemispheric ceramic tip with a diameter of 2 mm (Thorlabs, PKCESP) glued on a piezoelectric transducer (PZT) (Thorlabs, PA4CEW). The probe is mounted on a translation stage and brought to a physical contact with a sample. The broadband piezoelectric probe allows to use virtually any stimulus waveform to generate mechanical waves. Three representative waveforms are shown in Fig. 1(b). For a given peak-to-peak ratio, pure tones (first row) provide the maximum amount of power in an individual frequency. An impulse waveform (second row) allows broadband measurement in the frequency domain. A chirped signal (third row) provides nearly uniform, maximum power density over a broad frequency range for a fixed peak-to-peak amplitude [34–36]. Among the three types of stimulus, the most straightforward approach is using pure tones at discrete, selected frequencies, whereas the chirp stimulus is particularly useful to obtain a continuous dispersion curve over frequency.
To utilize the full Nyquist bandwidth (22.5 kHz), the OCT system was operated in the M-B mode: m consecutive A-lines are acquired at any specific transverse location (x coordinate). Each line is formed of 1024 pixels along the z axis. At each location along the x axis, a complete time series (M-scan) is acquired for a single or repeated mechanical stimulus. Typically, each individual stimulus was 5 ms long (225 lines) and repeated 64 times to improve measurement sensitivity through averaging. After the M-scan, the OCT beam is moved to the next location along the x axis, and the measurement is repeated. A single-line scan typically comprises 96 transverse positions over a range of 10 mm. The complete resulting measurement is thus comprised of 1024 pixels along the z axis, 64 repetitions of a 5 ms stimulus, and 96 pixels along the x axis, which yields a total measurement time of approximately 30 s. The galvanometer scanners have a finite stabilization time of approximately 200 μs, during which the mechanical jitter of the scan mirror causes spurious phase fluctuations in the OCT signal. To circumvent this problem, A-lines obtained during this transient period were simply ignored in post-processing.
Fig. 1 (a) Schematic of the contact probe and measurement setup. (b) Representative stimulus waveforms used in this study (left) and their frequency contents.
3. Data processing and validation using an elastomer
3.1. Displacement field
A Fourier transform of raw OCT inteferometric data from the optical wavenumber domain to axial position (z), repeated for each time (t) of A-line and lateral (x, y) point, produces a complex OCT tomogram A (r, t), where r = (x, y, z). The magnitude of this multidimensional array, typically in log scale, produces structural OCT images. The phase ϕ (r, t) of the complex A(r, t) is used to compute the axial component of the displacement field uz(r, t) using the following set of equations:
where ROI is a small neighborhood around r, λ0 is the average free-space wavelength of the OCT beam, nm and n0 are the refractive indices of the sample medium and the air above the sample, respectively, and ϕ(rsurf, t) is the phase at the air-sample boundary. Note that although the displacement field u(r, t) is a vector quantity, phase-sensitive OCT can only detect its component along the optical beam uz(r, t).To visualize wave propagation along the transverse coordinates, a z-projection of the displacement field was used. First, the top sample surface is identified using an automatic segmentation algorithm (or manually when the algorithm failed in low SNR portions of the images). Then, a ±5 pixels window around the top surface is defined and the signal in the region is summed (Eq. (2)). This step acts as brightness-weighted averaging of the displacement field.
Figure 2 shows the projected displacement field measured in an elastomer sample obtained with three different stimulus waveforms. The tissue-simulating silicone rubber phantom was prepared by mixing part A and part B precursors (Smooth-On Inc., EcoFlexTM 00-10) with a 1:1 volume ratio and casting the mixture onto a cylindrical container with a diameter of 10 cm and a height of 5 cm. The sample was removed from the mold after curing overnight. Mechanical waves are launched at the contact with the PZT probe and propagate radially in the surface of the homogeneous medium. The magnitude of the displacement is about 100 nm at the origin and decreases with r. It is clearly seen from the single-tone profiles that, as the frequency increases, the attenuation increases and the wavelength decreases. As a result, the number of detectable wavelengths is nearly constant (about 2–3 full cycles).
Fig. 2 Displacement field of elastic waves at the surface of an elastomer phantom for different waveforms. (a) Several time frames selected from a full data set for single 100 μs-long Gaussian impulse stimulus. (b) Corresponding time frames for a Schroeder-type chirp [34]. (c) Snapshots for pure tones at 1, 2, 4, 6, and 8 kHz, respectively.
3.2. Dispersion measurement
From the measured displacement field, the wave velocity is calculated via a couple of Fourier transform steps. Briefly, a two-dimensional (2D) Fourier transform of the displacement field (time t to frequency f and distance r to wavenumber k) produces a dispersion relation k(f) in the f-k plane. The phase velocity c is obtained from c = 2πf/k. For a chirped tone or any arbitrary stimulus waveform, the processing is identical to an impulse stimulus except that an additional step is needed to compute a cross-correlation function.
In Figure 3, we illustrate more detail of the data processing using actual experimental data acquired with a chirp stimulus (0–10 kHz bandwidth) on the isotropic, semi-infinite elastomer sample. We show intermediate steps in Figs. 3(b)–3(d) to illustrate the principle of data processing. However, in practice, one can skip these steps and move directly to the f-k domain (Fig. 3(e)) using a 2D Fourier transform of the raw data in the t-x domain (Fig. 3(a)).
Fig. 3 Processing steps to compute the mechanical wave dispersion of a sample applied to a tissue-simulating phantom. (a) Displacement field as a function of the time and transverse propagation distance (t-r plane). (b) Impulse response calculated from cross-correlation with the stimulus waveform. (c–d) A Fourier transform of the cross-correlation function over time, showing the magnitude, (c), and phase, (d), in the f-r plane. (e) Magnitude of the Fourier transform of the Fourier-domain cross-correlation function over space, revealing the dispersion curve of the lowest-order guided wave (bright curve). (f) Phase velocity obtained from the dispersion map as a function of frequency.
Figure 3(a) displays the measured displacement field at the top surface of the sample (e.g. Fig. 2(b)) in the t-r space: uz(r, t), where r denotes the propagation distance from the center of the vibrating tip along the surface. In this experiment, the data was obtained using a line scan along the x axis, so r = x. For the cornea, we note that r would be defined along the curved surface.
The chirped displacement signal can be compressed to an impulse-like waveform using a cross-correlation operation:
where τ is a time delay variable, ★ denotes the correlation operation, s(t) is the stimulus waveform, and γu,s, is the cross-correlation between u and s. This step makes the chirp approach equivalent to the impulse response measurement, but offers an advantage of enhanced signal-to-noise ratio (for a given peak mechanical displacement) because the power in the chirp is distributed over a longer period of time. This strategy is commonly used in radar imaging [37], and was recently proposed in the context of traveling-wave OCE [35]. There are several options to use as the reference waveform s(t). It can be either the displacement of the unloaded vibrating probe tip measured using the OCT, the displacement field measured at some reference point in the sample (alternatively, other reference material), or simply the voltage waveform applied to the probe if the frequency response of the device is sufficiently flat. We have found that these three options yielded identical results in practice. Figure 3(b) shows the calculated cross-correlation γu,s(r, τ), illustrating the equivalence between the chirp and impulse stimuli.The Fourier-domain function of the cross-correlation field (or the displacement field for the case of impulse stimulus) is given by
where Γu,s, Uz, and S are the Fourier transformations of γu,s, uz, and s, respectively, and * denotes complex conjugation. The magnitude and phase of the cross-spectrum Γu,s(r, f) are plotted in Figs. 3(c) and 3(d). We highlight that, in practice, we performed the correlation operation in the frequency domain using Eq. (6), which implicitly enforces the periodicity of each time series and naturally describes the repeated stimulus of our M-B scan acquisition configuration.The patterns in the f-r plane appear to follow hyperbolic curves suggesting that Γu,s is a function of ∼ fr. In other words, the magnitude and phase are approximately constant for a given f · r product. This behavior is expected for the phase considering that the equi-phase contours of a non-dispersive plane wave uz ∼ eikr would be given by ϕ = kr = (2π/c) fr. Any deviation from the hyperbolic relation is due to dispersion. The hyperbolic relation of magnitude indicates the frequency dependence in attenuation. In fact, the attenuation coefficients of compression and shear waves in soft viscoelastic materials, including tissues, have been measured to increase nearly linearly with frequency in the ultrasonic frequency range. [38,39].
Careful inspection of the impulse response plot (Fig. 3(b)) reveals that the mechanical stimulus excites not only the main Rayleigh-type wave (red), but also a weaker, faster-traveling wave (blue) that is out of phase with the main wave. A detailed investigation of the origin of the fast waves is outside the scope of this article. It is important to notice that the interference of the fast wave with the main wave shows up in the amplitude and phase of the Fourier-domain correlation function (Figs. 3(c) and 3(d)). In previous studies, the frequency-dependence of corneal wave velocity was calculated using the slope of the phase profile over distance [27,29,30]. However, we found that this method becomes unreliable, particularly at frequencies above 5 kHz because the interference causes the phase profile over distance to become nonlinear. The nonlinear slope yields position-dependent velocities even in homogeneous materials. When the amplitude of the fast wave is comparable with or greater than that of the main wave, the wave velocity can be greatly overestimated.
Our solution to this problem is to perform a second Fourier transform along the spatial axis (transforming distance r to spatial frequency or wavenumber k). The result is shown in Fig. 3(e). The magnitude of the Fourier transform in the f-k plane reveal the dispersion curves k(f) of all wave components. The main wave component can be tracked using a simple bounded max search at every frequency and separated from other waves. Zero-padding by a factor 8 was used in the Fourier transform so that the maxima could be found with greater accuracy.
Finally, the phase velocity c = 2πf/k is obtained from the dispersion curve. The result is plotted in Fig. 3(f). In principle, the group velocity can also be obtained from the local slope of the dispersion curve. Again, we note that the dispersion curve (Fig. 3(e)) can be directly computed as
where U(f, k) is the 2D Fourier transform of the displacement, and S*(f) is the complex-conjugate Fourier transform of the reference stimulus waveform.To evaluate the fidelity of the measurement, we measured the wave velocity using an impulse stimulus and pure tones at several discrete frequencies. All the results are plotted in Fig. 4(a). We can see that the data obtained with the chirp (cyan circles) are consistent with those obtained using a short 100 μs Gaussian pulse (magenta circles). The chirp stimulus tends to allow slightly less noisy readings due to the improved signal to noise ratio. We also verified that both approaches based on broadband stimuli produce the same results with pure tones (yellow squares). This equivalence is expected for linear materials, which is a reasonable assumption for low amplitude waves traveling in viscoelastic materials. The measurement data at frequencies below 2 kHz have some apparent noise or artifact mainly due to the fast waves that were not well separated in the f-k plane and therefore not completely removed. The data in 2–10 kHz show high fidelity. In this high frequency range, the measured velocities increase from ∼ 7.4 to 9.6 m/s. The frequency dependence is attributed to the material dispersion of the elastomer. The waveguide dispersion is negligible because the 5 cm thickness of the sample is much greater than the wavelength in the frequency range (3.7 mm at 2 kHz to 0.96 mm at 10 kHz). In this regime, the lowest-order Lamb waves propagate along the surface, known as the Rayleigh surface wave.
Fig. 4 (a) Measured phase velocity of silicone rubber using three different stimulus waveforms: chirp (Schroeder-type [34], with frequency components between 0 and 10 kHz by steps of 180 Hz), impulse (Gaussian-shaped with a 4σ pulse length of 100 μs), and single tones at frequencies 2, 4, 6, 8 and 10 kHz. (b) Comparison to wave speeds measured by rheometry in a low frequency range of 0.1–100 Hz). Inset, complex shear modulus measured using a commercial shear rheometer.
Standard mechanical test equipment is operated at lower frequencies below 100 Hz. Therefore, validation of our measurement is not straightforward. As indirect assessment, we measured the shear storage modulus (μ′) and loss modulus (μ″) of the elastomer material using a commercial shear rheometer (AR-G2, TA Instruments) in the frequency range of 0.1 Hz and 100 Hz. The data is shown in Fig. 4(b). Using the measured complex shear modulus, μ = μ′ + iμ″ and Eq. (1), we obtain the complex bulk shear speed. The real-part phase velocity is given by [23], which can be expressed as:
where and tan δ = μ″/μ′. The Rayleigh surface wave in incompressible materials has a velocity cR ≈ 0.955 cs. We calculated the Rayleigh wave phase velocity from the rheometer data and plot the result in Fig. 4(b), along with the OCE data. The two sets of data are connected remarkably well, indirectly supporting the accuracy of our OCE measurement. As described above, the rapid frequency dependence in the OCE data at 0.2–2 kHz is likely an artifact due to fast-traveling waves.4. Elastography of the cornea
4.1. Displacement field and dispersion
Figure 5(a) shows a schematic of the setup for corneal elastography. Porcine eyes (10 pairs) were obtained less than 8 h postmortem, and connective tissues such as fat, muscles and the optic nerve were dissected out for consistent boundary conditions and weight measurements. After placing an eye in a custom-made holder mounted on a manual 3-axis positioning stage, the IOP is set to 15 mmHg using a water column. When the IOP is stabilized, the eye is aligned with the OCT beam, and the PZT probe is brought to the eye about 4 mm away from the corneal apex. After the probe tip makes a contact with the corneal epithelium, it was gently advanced by ∼ 100 μm further to ensure the mechanical contact to be maintained during stimulus actuation. Experiments were performed on a total of 20 eyes (10 pairs).
Fig. 5 (a) Schematic of elastography experiment performed on a porcine eye. (b) Time-dependent displacement field in the cornea from a impulse stimulus. The displacement field is overlaid over the standard OCT intensity image.
Figure 5(b) shows the vibrography images of the cornea in response to an 0.1-ms-long impulse stimulus. At the onset (0 ms), the downward (positive displacement) movement initiates from the point of contact and propagates along the corneal tissue. The displacement field is quite uniform across the depth of the cornea, indicating the guided wave is predominantly the lowest-order, anti-symmetric (flexural) Lamb wave. At 1 ms, the leading edge (red region) of the wave reaches the corneal apex 4 mm away from the probe tip (i.e. c ≈ 4 m/s).
We acquired OCE dataset using the same chirp waveform and processing steps described in Section 3. Figure 6 shows the measured displacement field in the t-r domain, the magnitude of the cross-correlation function in the f-r domain, and the dispersion map in the f-k domain. Two dispersive branches are apparent in Figs. 6(a) and 6(c), corresponding to the lowest-order Lamb wave and a fast wave, respectively. The resulting interference patterns are visible in Fig. 6(b). The identity of the fast wave(s) is not well understood but is thought to come from the excitation of compression waves, which have bulk velocities of ∼ 1620 m/s in the cornea [22] and ∼ 1540 m/s in the aqueous humor, and possibly also from high-order guided waves. The dispersion curve of the dominant Lamb wave (bright curve) was isolated and used to compute the phase velocity dispersion of the cornea.
Figure 6(d) shows the velocity dispersion data obtained from a single eye, measured 3 times. Between each measurement, the eye was removed from the holder and put back in place, and the contact of the PZT probe was re-established, introducing substantial difference in the contact position and OCT scan line. The result demonstrates a good repeatability of measurement with a standard deviation of ±3.9 %.
Figure 6(e) shows the velocity dispersion curves for all of the 20 measured eyeballs. The data show an inter-sample variation of ±7.2 % in wave velocity for the porcine eyes at a constant IOP. Besides the small quantitative differences, all the curves show common features. At low frequencies below 1.5 kHz, the velocity increases with decreasing frequency. We believe this is an artifact due to the finite size of the cornea. We will discuss it more later. At higher frequencies above 2 kHz, the velocity linearly increases with frequency. In the following section, we analyze this frequency range to extract viscoelastic properties of the corneal tissue.
Fig. 6 Phase velocity measurement of porcine corneas at the same IOP level of 15 mmHg. (a) Displacement field. (b) Cross-spectrum magnitude. (c) Dispersion map obtained by a 2D Fourier transform of the displacement field. (d) Frequency-dependent phase velocity of the same porcine cornea measured 3 times. (e) Frequency-dependent phase velocity of 20 different porcine eye samples.
4.2. Viscoelasticity
One promising approach to extract the biomechanical parameters of corneal tissue is to fit the phase velocity dispersion curves with those predicted by a Lamb-wave waveguide model, such as the modified Rayleigh-Lamb model combined with a viscoelatic material model such as the Kelvin-Voigt (KV) model [29, 40]. This approach has been applied to data obtained with air puff stimuli with narrowband frequency contents below 1 kHz [29]. To good approximation, the cornea can be modeled as an infinite, flat plate with a thickness h and density ρ, bounded by the air (or vacuum) on the top surface and water at the bottom surface. The material properties are assumed to be isotropic and homogeneous. The shear modulus in the Kelvin-Voigt model is described as [41]:
where the real part μ0 corresponds to shear modulus at zero frequency (μ = μ0 at f = 0), and the imaginary part originates from shear viscosity η. The bulk shear wave phase velocity is given by [23]: where .The traveling-wave analysis in OCE described here measures the complex wavenumber k of the wave, whose real part defines the wave velocity c and imaginary part describes the wave attenuation α. Applying the surface boundary conditions to the elastodynamic wave equation yields a 5 by 5 matrix equation whose determinant must be zero for guided waves (see equation 17 of reference [29]). In Fig. 7(a), we show a computational result of this model for h = 0.8 mm, ρ = 1 g/cm3, μ0 = 25 kPa, η = 0, and λ = 2.2 GPa. By taking η = 0, the result is consistent with the well-known dispersion curves of various Lamb modes in isotropic, elastic plates. The modes are labeled as either quasi-antisymmetric (An) or quasi-symmetric (Sn), where the index n indicates the order of the mode. Due to the presence of different media at the top and bottom boundaries, the modes do not have pure odd or even symmetry. At low frequencies, the fundamental branch A0 represents the flexural plate wave. At high frequencies, it asymptotically transforms to become the Scholte wave with a velocity cSch ≈ 0.846 cs, which is a surface wave traveling at the cornea-fluid interface [23]. The S0 branch becomes the Rayleigh surface wave traveling the cornea-air boundary at a velocity cRay ≈ 0.955 cs at high frequencies, whereas the other branches tend towards the bulk shear wave velocity.
Fig. 7 Parameter estimation using the Lamb wave model. (a) Example Lamb wave spectrum for typical porcine cornea properties, neglecting viscosity. (b) Lamb wave model fitted to the measured porcine cornea phase velocities. (c) Fit error estimation using the initial guess variation method. Each marker represents the fit results for a different set of initial guess parameters. Probability density function (PDF) contours are also shown to help visualization. Each panel uses a different frequency ranges for fitting.
We fitted the A0-mode wave solution to the average phase velocity curve in the 2–10 kHz range obtained from the 20 porcine cornea samples, using μ0 and η as the only two free parameters and the trust-region-reflective least-square algorithm. The result is shown in Fig. 7(b), where the best fit was obtained with μ0 = 20.5 kPa and η = 0.28 Pa s. We found that the model and the data agree well over the fitting range. The positive slope in phase velocity is accounted for by viscous damping.
To determine the precision of μ0 and η values, we have run the fitting calculation 100 times using different initial guess values of μ0 randomly chosen between 10 and 50 kPa and of η between 0.05 and 0.45 Pa s. The fitting algorithm converges to a slightly different solution depending on the initial guess. The distribution of solutions in the μ0-η space is shown in Fig. 7(c). When the entire data set in the 2–10 kHz range was used, the analysis yields a narrow distribution with standard deviations of <10 % for both μ0 and η. However, partial data for narrower frequency ranges produce more inaccurate estimation, especially for η. Note that the distributions of μ0 and η are not independent. This result highlights that broadband excitation stimulus is necessary for measuring both the stiffness and viscosity from velocity dispersion.
The wave attenuation coefficient α can also be obtained from the KV model using the following equation [23]:
We measured the attenuation of the displacement field at each frequency by fitting a decaying exponential to the magnitude expressed as a function of r. The 1/r-type diffraction loss was taken into account by scaling the amplitude by prior to fitting. The attenuation value obtained is shown in Fig. 8 for each of the 20 porcine corneas as grey lines and with the mean and 2-standard-deviation statistics shown in green. The measured attenuation in the range of 0–4 kHz compares well with the near-quadratic frequency-dependence predicted by Eq. (11) for μ0 = 20.5 kPa and η = 0.28 Pa s (dotted red line in Fig. 8). To account for waveguide loss, we solved the Lamb wave characteristic equation using the KV-model material property and calculated attenuation from the complex propagation constant (solid red line in Fig. 8). The reasonably good correspondence in the 0–3.5 kHz range supports the validity of our analysis based on the KV and Lamb wave model. There is a small offset of ∼ 0.15 mm−1 between the data and theory, which may be attributed to the deviations of the cornea from the KV and Lamb models. Note that because the attenuation was measured in the f-r, rather than f-k, domain, the interference of the fast waves was not filtered out, resulting in unreliable attenuation measurement in 4–10 kHz.Fig. 8 Wave attenuation measurement of porcine corneas as a function of frequency. Gray lines represent the measured wave attenuation data for each of the 20 porcine eyes at an IOP of 15 mmHg, measured by fitting an exponential decay function to the wave magnitude at each frequency. Green line and shaded region represent the mean and 2 standard deviation error. The solid red curve represents the attenuation predicted by the KV and Lamb models with μ0 = 20.5 kPa and η = 0.28 Pa s previously obtained. The dotted line represents the attenuation in the absence of waveguiding as predicted by Eq. (11).
4.3. Correlation of wave velocity with anatomical parameters
Figure 9 shows the distribution of shear modulus across the 20 eyes, plotted as a function of the eyeball weight, central corneal thickness, and central anterior radius of curvature. Left and right eyes of the same animal are represented with left- and right-pointing triangles of matching color in Figs. 9(a)–9(c), and are plotted against each other in Fig. 9(d). We found rather strong negative correlation with the weight (Pearson’s r2 = 0.56) and mild to weak correlation with the curvature (r2 = 0.11). The origin of the clear weight dependence could be physiological, but other weight-dependent geometrical effect might have contributed via waveguide dispersion [42]. We found no correlation between shear modulus and corneal thickness (r2 = 0.048), which is expected given that the Rayleigh-Lamb model takes into account the effects of thickness on wave velocity. Other parameters influencing thickness, such as corneal swelling, could potentially be appreciated, but the relatively narrow range of corneal thicknesses and the absence of significant correlation suggests that this parameter was adequately controlled during our experiments. Comparison between the left and right eye showed a moderately strong correlation (r2 = 0.29) with a slope close to unity (dotted black line, Fig. 9d).
Fig. 9 Measured static shear modulus, μ0 of 10 pairs of porcine eyes (20 eyes) as a function of (a) weight, (b) central corneal thickness (CCT), and (c) anterior corneal radius of curvature. (d) Comparison between the left and right eyes of matching pairs.
5. Discussion and conclusion
We found that PZT-based OCE is well suited to generate and measure the propagation of broadband mechanical waves in samples, with advantages of the wide bandwidth, flat frequency response, and low-cost of piezoelectric actuators. Using this technique, we have measured the phase velocity dispersion, stiffness, and viscosity of porcine corneas with high fidelity, although we do not have alternative, gold-standard values to validate the data quantitatively. For clinical applications, non-contact approaches such as air puffs and acoustic radiation force [43,44] are appealing. However, using local anesthetics, contact with the cornea is a well tolerated medical procedure used in tonometry, pachymetry, and some embodiements of OCE [45]. The PZT-based contact probe may also be a useful tool for other potential applications of OCE in dermatology and other medical areas, as well as laboratory measurement.
Our result indicate that 2–10 kHz is a suitable frequency window for traveling-wave OCE for the cornea. Below 2 kHz, the data deviates strongly from the A0 branch of the theoretical model, as previously observed [30]. It appears unlikely that this is the signature of the S0 mode, because the velocity dispersion and the spatial pattern of displacement field deviate from those of the symmetric mode. One possible explanation is that low frequency excitation is susceptible to excite resonance modes of the cornea or the whole eyeball [46]. Resonance modes are standing waves, making them appear to have infinite wave velocity. These fast waves interfere in the phase velocity analysis, causing overestimation of the shear wave velocity. While the f-k domain processing approach proposed here is designed to reject these waves, it is less effective when multiple branches overlap at low frequency below 2 kHz for the cornea. Another explanation is that at frequencies below 1 kHz, the wavelength becomes comparable with the size of the cornea. The wave propagation could therefore be more influenced by any inhomogeneity in the mechanical properties of the cornea and reflection from the cornea-sclera interface.
Fitting the velocity dispersion with the A0 wave is only an approximation. In principle, the contact stimulus on the surface could excite a combination of multiple modes including the S0 and higher order A and S modes in the far field if the excitation frequency is above their cutoffs. In the near field, an infinity of complex modes are also required to fully describe the wave field in the vicinity of the source, where the bulk shear and compression waves have not yet developed into true guided waves. Our f-k domain approach can separate the different wave components if the propagation distance is long enough, but this is often not the case due to viscous damping and the finite extent of the cornea itself. Furthermore, the wave velocity is similar between the A0 and S0 modes above 5 kHz. The Lamb wave model also neglects the effects of IOP-induced stress, corneal curvature, anisotropy, and inhomogeneity of the material properties, resulting in biased estimations of shear modulus and viscosity [42].
In the KV model, we calculate a 10-dB attenuation distance of the shear wave, L, defined as the propagation distance over which the displacement decreases to 10%:
The ratio of L to the wavelength λ = c/f, or the number of wave cycles (N = L/λ) within the effective propagation distance is: If we impose N > 1 as a necessary condition for reliable measurement, we get ζ < 1.31, which is (2π f η) < 0.85μ0. For porcine corneas, using our measurement data of μ0 = 20.5 kPa and η = 0.28 Pa s, we find the high-frequency limit to be 9.9 kHz. Although the 10-dB attenuation distance is somewhat arbitrary and will depend on the sensitivity of the instrument and the magnitudes of spurious fast waves, our empirical experience suggests that this condition is reasonable. Another difficulty encountered at higher frequencies above 10 kHz is low displacement amplitude, which for a fixed power, decreases with frequency.The above criteria for optimal frequency ranges imply an important condition for samples. For highly viscous samples, the high-frequency limit would come down to a lower value. When this value is lower than the low-frequency limit imposed by the size of the sample, the traveling-wave analysis would not yield reliable measurement, and different algorithms may be needed to extract both the elastic and viscous properties.
The shear modulus values of porcine corneas we measured are comparable with previous values measured by other OCE techniques using air puff [29,31] and acoustic radiation force [30], as well as ultrasound elastography [19]. There are discrepancies of up to a factor of two among reported values, part of which could be attributed to differences in wave frequency, IOP, processing algorithm, and tissue models (e.g. linear elastic).
In summary, this study describes the generation and measurement of broadband mechanical waves in the cornea to characterize the phase velocity dispersion in the cornea and estimate its viscoelastic material parameters. Our results suggest that a wide frequency band improves the accuracy of elasticity and viscosity estimation by fitting the phase velocity with a Lamb wave model. This research shows that the combination of phase-sensitive OCT vibration sensing with a piezoelectric actuator is a promising approach to measure the mechanical properties of the cornea and other samples.
Funding
National Institutes of Health (NIH) (P41-EB015903, R01-EY025454); National Science Foundation (NSF) (CMMI-1562863)
Acknowledgments
The authors wish to thank Dr. Zeinab Hajjarian for her help with rheometry measurements and Prof. Laurent Demanet for helpful discussions.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. T. T. Andreassen, A. Hjorth Simonsen, and H. Oxlund, “Biomechanical properties of keratoconus and normal corneas,” Exp. Eye Res. 31, 435–441 (1980). [CrossRef] [PubMed]
2. I. S. Nash, P. R. Greene, and C. S. Foster, “Comparison of mechanical properties of keratoconus and normal corneas,” Exp. eye research 35, 413–424 (1982). [CrossRef]
3. C. Edmund, “Corneal elasticity and ocular rigidity in normal and keratoconic eyes,” Acta ophthalmologica 66, 134–140 (1988). [CrossRef] [PubMed]
4. J. Liu and C. J. Roberts, “Influence of corneal biomechanical properties on intraocular pressure measurement,” J. Cataract. & Refract. Surg. 31, 146–155 (2005). [CrossRef]
5. P. D. Jaycock, L. Lobo, J. Ibrahim, J. Tyrer, and J. Marshall, “Interferometric technique to measure biomechanical changes in the cornea induced by refractive surgery,” J. Cataract. & Refract. Surg. 31, 175–184 (2005). [CrossRef]
6. R. R. Krueger and W. J. Dupps, “Biomechanical effects of femtosecond and microkeratome-based flap creation: prospective contralateral examination of two patients,” J. refractive surgery 23, 800–807 (2007).
7. W. J. Dupps and C. Roberts, “Effect of acute biomechanical changes on corneal curvature after photokeratectomy,” J. refractive surgery 17, 658–669 (2001).
8. T. Seiler, M. Matallana, S. Sendler, and T. Bende, “Does Bowman’s layer determine the biomechanical properties of the cornea?” Refract. & corneal surgery 8, 139–142 (1992).
9. A. Elsheikh, D. Wang, and D. Pye, “Determination of the modulus of elasticity of the human cornea,” J. refractive surgery 23, 808–818 (2007). [CrossRef]
10. A. Elsheikh, D. Alhasso, and P. Rama, “Biomechanical properties of human and porcine corneas,” Exp. Eye Res. 86, 783–790 (2008). [CrossRef] [PubMed]
11. S. Shah, M. Laiquzzaman, R. Bhojwani, S. Mantry, and I. Cunliffe, “Assessment of the biomechanical properties of the cornea with the ocular response analyzer in normal and keratoconic eyes,” Investig. Ophthalmol. Vis. Sci. 48, 3026–3031 (2007). [CrossRef]
12. D. Ortiz, D. Piñero, M. H. Shabayek, F. Arnalich-Montiel, and J. L. Alió, “Corneal biomechanical properties in normal, post-laser in situ keratomileusis, and keratoconic eyes,” J. Cataract. & Refract. Surg. 33, 1371–1375 (2007). [CrossRef]
13. A. K. Lam, D. Chen, and J. Tse, “The usefulness of waveform score from the ocular response analyzer,” (2010).
14. A. Kotecha, “What biomechanical properties of the cornea are relevant for the clinician?” Surv. Ophthalmol. 52, S109–S114 (2007). [CrossRef] [PubMed]
15. J. Kerautret, J. Colin, D. Touboul, and C. Roberts, “Biomechanical characteristics of the ectatic cornea,” J. Cataract. Refract. Surg. 34, 510–513 (2008). [CrossRef] [PubMed]
16. D. A. Luce, “Determining in vivo biomechanical properties of the cornea with an ocular response analyzer,” J. Cataract. Refract. Surg. 31, 156–162 (2005). [CrossRef] [PubMed]
17. V. S. De Stefano and W. J. Dupps, “Biomechanical diagnostics of the cornea,” Int. Ophthalmol. Clin. 57, 75–86 (2017). [CrossRef] [PubMed]
18. K. W. Hollman, S. Y. Emelianov, J. H. Neiss, G. Jotyan, G. J. Spooner, T. Juhasz, R. M. Kurtz, and M. O’Donnell, “Strain imaging of corneal tissue with an ultrasound elasticity microscope,” Cornea 21, 68–73 (2002). [CrossRef] [PubMed]
19. M. Tanter, D. Touboul, J. L. Gennisson, J. Bercoff, and M. Fink, “High-resolution quantitative imaging of cornea elasticity using supersonic shear imaging,” IEEE Transactions on Med. Imaging 28, 1881–1893 (2009). [CrossRef]
20. G. Scarcelli and S. H. Yun, “Confocal Brillouin microscopy for three-dimensional mechanical imaging,” Nat. Photonics 2, 39–43 (2008). [CrossRef]
21. G. Scarcelli, S. Besner, R. Pineda, P. Kalout, and S. H. Yun, “In vivo biomechanical mapping of normal and keratoconus corneas,” JAMA ophthalmology 133, 480–482 (2015). [CrossRef] [PubMed]
22. P. Shao, T. G. Seiler, A. M. Eltony, A. Ramier, S. J. J. Kwok, G. Scarcelli, R. Pineda II, and S.-H. Yun, “Effects of corneal hydration on Brillouin microscopy in vivo,” Investig. Opthalmology & Vis. Sci. 59, 3020 (2018). [CrossRef]
23. M. A. Kirby, I. Pelivanov, S. Song, Ł. Ambrozinski, S. J. Yoon, L. Gao, D. Li, T. T. Shen, R. K. Wang, and M. O’Donnell, “Optical coherence elastography in ophthalmology,” J. Biomed. Opt. 22, 121720 (2017). [CrossRef]
24. S. Marcos, G. Scarcelli, A. Ramier, and S. H. Yun, “Probing ocular mechanics with light,” in Encyclopedia of Modern Optics (Second Edition), (Elsevier, Oxford, 2018), pp. 140–154. [CrossRef]
25. K. V. Larin and D. D. Sampson, “Optical coherence elastography – OCT at work in tissue biomechanics [Invited],” Biomed. Opt. Express 8, 1172 (2017). [CrossRef] [PubMed]
26. R. K. Manapuram, S. R. Aglyamov, F. M. Monediado, M. Mashiatulla, J. Li, S. Y. Emelianov, and K. V. Larin, “In vivo estimation of elastic wave parameters using phase-stabilized swept source optical coherence elastography,” JBO Lett. 17, 15–18 (2012).
27. S. Wang and K. V. Larin, “Noncontact depth-resolved micro-scale optical coherence elastography of the cornea,” Biomed. optics express 5, 3807–3821 (2014). [CrossRef]
28. Y. Qu, T. Ma, Y. He, J. Zhu, C. Dai, M. Yu, S. Huang, F. Lu, K. K. Shung, Q. Zhou, and Z. Chen, “Acoustic radiation force optical coherence elastography of corneal tissue,” IEEE J. Sel. Top. Quantum Electron. 22, 288–294 (2016). [CrossRef]
29. Z. Han, J. Li, M. Singh, C. Wu, C.-h. Liu, R. Raghunathan, S. R. Aglyamov, S. Vantipalli, M. D. Twa, and K. V. Larin, “Optical coherence elastography assessment of corneal viscoelasticity with a modified Rayleigh-Lamb wave model,” J. Mech. Behav. Biomed. Mater. 66, 87–94 (2017). [CrossRef]
30. Ł. Ambroziński, S. Song, S. J. Yoon, I. Pelivanov, D. Li, L. Gao, T. T. Shen, R. K. Wang, and M. O’Donnell, “Acoustic micro-tapping for non-contact 4D imaging of tissue elasticity,” Sci. Reports 6, 38967 (2016). [CrossRef]
31. Z. Han, S. R. Aglyamov, J. Li, M. Singh, S. Wang, S. Vantipalli, C. Wu, C.-h. Liu, M. D. Twa, and K. V. Larin, “Quantitative assessment of corneal viscoelasticity using optical coherence elastography and a modified Rayleigh–Lamb equation,” J. Biomed. Opt. 20, 020501 (2015). [CrossRef]
32. A. Ramier, J. T. Cheng, M. E. Ravicz, J. J. Rosowski, and S.-H. Yun, “Mapping the phase and amplitude of ossicular chain motion using sound-synchronous optical coherence vibrography,” Biomed. Opt. Express 9, 5489 (2018). [CrossRef] [PubMed]
33. B. Braaf, K. A. Vermeer, V. A. D. Sicam, E. van Zeeburg, J. C. van Meurs, and J. F. de Boer, “Phase-stabilized optical frequency domain imaging at 1-μm for the measurement of blood flow in the human choroid,” Opt. Express 19, 20886 (2011). [CrossRef] [PubMed]
34. M. Schroeder, “Synthesis of low-peak-factor signals and binary sequences with low autocorrelation,” IEEE Transactions on Inf. Theory 16, 85–89 (1970). [CrossRef]
35. T.-M. Nguyen, S. Song, B. Arnal, E. Y. Wong, Z. Huang, R. K. Wang, and M. O’Donnell, “Shear wave pulse compression for dynamic elastography using phase-sensitive optical coherence tomography,” J. Biomed. Opt. 19, 016013 (2014). [CrossRef]
36. T.-M. Nguyen, B. Arnal, S. Song, Z. Huang, R. K. Wang, and M. O’Donnell, “Shear wave elastography using amplitude-modulated acoustic radiation force and phase-sensitive optical coherence tomography,” J. biomedical optics 20, 16001 (2015). [CrossRef]
37. J. R. Klauder, A. C. Price, S. Darlington, and W. J. Albersheim, “The theory and design of chirp radars,” Bell Syst. Tech. J. 39, 745–808 (1960). [CrossRef]
38. P. Wells, “Absorption and dispersion of ultrasound in biological tissue,” Ultrasound Medicine & Biol. 1, 369–376 (1975). [CrossRef]
39. E. L. Madsen, H. J. Sathoff, and J. A. Zagzebski, “Ultrasonic shear wave properties of soft tissues and tissuelike materials,” The J. Acoust. Soc. Am. 74, 1346–1355 (1983). [CrossRef] [PubMed]
40. I. Z. Nenadic, M. W. Urban, S. A. Mitchell, and J. F. Greenleaf, “Lamb wave dispersion ultrasound vibrometry (LDUV) method for quantifying mechanical properties of viscoelastic solids,” Phys. Medicine Biol. 56, 2245–2264 (2011). [CrossRef]
41. M. A. Meyers and K. K. Chawla, Mechanical Behavior of Materials (Cambridge University, Cambridge, 2008). [CrossRef]
42. Z. Han, J. Li, M. Singh, S. R. Aglyamov, C. Wu, C.-h. H. Liu, and K. V. Larin, “Analysis of the effects of curvature and thickness on elastic wave velocity in cornea-like structures by finite element modeling and optical coherence elastography,” Appl. Phys. Lett. 106, 233702 (2015). [CrossRef] [PubMed]
43. S. Wang, K. V. Larin, J. Li, S. Vantipalli, R. K. Manapuram, S. Aglyamov, S. Emelianov, and M. D. Twa, “A focused air-pulse system for optical-coherence-tomography-based measurements of tissue elasticity,” Laser Phys. Lett. 10, 075605 (2013). [CrossRef] [PubMed]
44. Ł. Ambroziński, I. Pelivanov, S. Song, S. J. Yoon, D. Li, L. Gao, T. T. Shen, R. K. Wang, and M. O’Donnell, “Air-coupled acoustic radiation force for non-contact generation of broadband mechanical waves in soft media,” Appl. Phys. Lett. 109, 043701 (2016). [CrossRef]
45. V. S. De Stefano, M. R. Ford, I. Seven, and W. J. Dupps, “Live human assessment of depth-dependent corneal displacements with swept-source optical coherence elastography,” PLOS ONE 13, e0209480 (2018). [CrossRef] [PubMed]
46. B. I. Akca, E. W. Chang, S. Kling, A. Ramier, G. Scarcelli, S. Marcos, and S. H. Yun, “Observation of sound-induced corneal vibrational modes by optical coherence tomography,” Biomed. Opt. Express 6, 3313 (2015). [CrossRef] [PubMed]
