Measuring mechanical wave speed, dispersion, and viscoelastic modulus of the cornea using optical coherence elastography.

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.


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][2][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][6][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][9][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][12][13][14][15][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][21][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][24][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][27][28][29][30]. There is a direct theoretical relationship between complex shear wave velocity c s 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, c s has frequency dependence (dispersion), and measurement of c s 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][30][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.

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][35][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

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 u z (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, n m and n 0 are the refractive indices of the sample medium and the air above the sample, respectively, and φ(r surf , 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 u z (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., EcoFlex TM 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).

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 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)). 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: u z (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 , U z , and S are the Fourier transformations of γ u,s , u z , 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 ∼ f r. 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 u z ∼ e ikr would be given by ϕ = kr = (2π/c) f r. 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.
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 c = 1/Re(c −1 s ) [23], which can be expressed as: where | µ| = µ 2 + µ 2 and tan δ = µ /µ . The Rayleigh surface wave in incompressible materials has a velocity c R ≈ 0.955 c s . 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. 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).  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.

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]: 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/cm 3 , µ 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 (A n ) or quasi-symmetric (S n ), 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 A 0 represents the flexural plate wave. At high frequencies, it asymptotically transforms to become the Scholte wave with a velocity c Sch ≈ 0.846 c s , which is a surface wave traveling at the cornea-fluid interface [23]. The S 0 branch becomes the Rayleigh surface wave traveling the cornea-air boundary at a velocity c Ray ≈ 0.955 c s at high frequencies, whereas the other branches tend towards the bulk shear wave velocity. We fitted the A 0 -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 √ r 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. 9(d). We found rather strong negative correlation with the weight (Pearson's r 2 = 0.56) and mild to weak correlation with the curvature (r 2 = 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 (r 2 = 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 (r 2 = 0.29) with a slope close to unity (dotted black line, Fig. 9d).

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 A 0 branch of the theoretical model, as previously observed [30]. It appears unlikely that this is the signature of the S 0 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 A 0 wave is only an approximation. In principle, the contact stimulus on the surface could excite a combination of multiple modes including the S 0 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 A 0 and S 0 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.