An approach to viscoelastic characterization of dispersive media by inversion of a general wave propagation model

In the characterization of elastic properties of tissue using dynamic optical coherence elastography, shear/surface waves are propagated and tracked in order to estimate speed and Young's modulus. However, for dispersive tissues, the displacement pulse is highly damped and distorted during propagation, diminishing the e®ectiveness of peak tracking approaches, and leading to biased estimates of wave speed. Further, plane wave propagation is sometimes assumed, which contributes to estimation errors. Therefore, we invert a wave propagation model that incorporates propagation, decay, and distortion of pulses in a dispersive media in order to accurately estimate its elastic and viscous components. The model uses a general ̄rst-order approximation of dispersion, avoiding the use of any particular rheological model of tissue. Experiments are conducted in elastic and viscoelastic tissue-mimicking phantoms by producing a Gaussian push using acoustic radiation force excitation and measuring the wave propagation using a Fourier domain optical coherence tomography system. Results con ̄rmed the e®ectiveness of the inversion method in estimating viscoelastic parameters in both the viscoelastic and elastic phantoms when compared to mechanical measurements. Finally, the viscoelastic characterization of a fresh porcine cornea was conducted. Preliminary results validate this approach when compared to other methods.


Introduction
Determining mechanical properties of tissue such as elasticity and viscosity is fundamental for better understanding and diagnosing pathological and physiological processes. 1,2In this regard, optical coherence tomography-based elastography (OCE) o®ers the possibility of a noninvasive, high-resolution, and high-contrast measurement of tissue biomechanical properties. 3,4For example, previous studies demonstrate the potential of OCE in assessing mechanical properties of di®erent tissues such as cornea, 5,6 skin, 7 breast, 8 and liver. 9In particular, a subset of dynamic OCE techniques uses short-duration pulses produced by a selected excitation source in order to produce mechanical wave propagation in the tissue being studied. 10Excitation sources include acoustic radiation force (ARF), air-pu® excitation, laser-based thermal expansion, and needles connected to piezoelectric vibrators, to name just a few. 11y tracking the propagating wave, Young's modulus and other biomechanical parameters can be calculated based on the estimation of the wave speed and the selection of the correct wave propagation model dictated by the boundary conditions of the sample. 11nfortunately, in many dispersive lossy tissues, propagation of shear or surface waves will be rapidly damped and distorted, complicating the attempts of estimating wave speed by using typical methodologies such as peak tracking. 12Moreover, the estimation of viscous parameters in addition to the classic elastic modulus is of great interest since it can provide of useful information in the diagnosis of diseases such as the discrimination between malignant and benignant liver tumors 13 and the characterization of glioblastomas in human brain. 146][17][18] Most of them utilize frequency-dependent wave speed measurements to ¯t with theoretical models of dispersion (rheological models), disregarding valuable information given by the wave attenuation process.Therefore, the use of model-independent approaches for the viscoelastic characterization of tissue is an important trending topic in elastography since it can provide an accurate intrinsic information without assumptions of the tissue mechanical behavior, which is in many cases unknown.
Few model-independent approaches have been proposed in ultrasound elastography, 12,[19][20][21][22] and in OCE. 23In ultrasound, Schmitt et al. 19 proposed a method for the calculation of the complex shear modulus by the estimation of shear wave dispersion and attenuation of a sinusoidal continuous plane wave.Nenadic et al. 20 proposed a 2D Fourier approach for the analysis of a cylindrically spreading continuous waves and the estimation of speed dispersion and attenuation.Kazemirad et al. 21estimated the complex shear modulus by ¯tting an analytical model of a continuous cylindrical shear wave with a shear pulse generated using ARF excitation.In OCE, Leartprapun et al. 23 proposed a model-independent method for the reconstruction of complex shear modulus from measurements of continuous hemi-spherical surface waves produce by ARF excitation in viscoelastic media.While these approaches calculate viscoelastic properties of tissue without using a rheological model, some assumptions made may produce biased results.For instance, in Ref. 19, the plane wave assumption can be di±cult to satisfy for most of the excitation methods.Also, Refs.20, 21, and 23, assume that waves will be observed in the asymptotic range far from the excitation source which may not be possible if the wave is highly attenuated by dispersive lossy media.Finally, Parker and Baddour 12 investigated the propagation of a cylindrical axisymmetric Gaussian shear wave in a viscoelastic media by proposing a full analytical model-independent solution that takes a ¯rst-order approach to dispersion.
In this paper, we invert a general wave propagation model following the approach of Parker and Baddour 12 that incorporates space-time propagation, decay, and distortion of pulses in a dispersive medium in order to accurately estimate the elastic and viscous components of such a medium.The model contemplates the initial shape of ARF push in space and time and uses a general ¯rst-order approximation of dispersion, avoiding the use of any particular rheological model of tissue.Experiments will be conducted in elastic and viscoelastic tissuemimicking phantoms by producing a Gaussian push using ARF excitation, and measuring the surface wave propagation using a Fourier domain optical coherence tomography (FD-OCT) system.Results found in the inversion method will be compared to mechanical measurements (MM) for validation.Finally, we will conduct a preliminary experiment in a fresh porcine cornea in order con¯rm the validity of our approach in real viscoelastic tissue.

Theory 2.1. Shear wave pulse in an in¯nite media
An ultrasound (US)-based ARF excitation beam with a Gaussian intensity axisymmetric pro¯le in cylindrical coordinates ðr; ; zÞ is given by in the r-direction, and extending uniformly in the axial (z) direction, where 2 ¼ half variance of the pulse, r ¼ ðx 2 þ y 2 Þ 1=2 , ¼ atanðy=xÞ as shown in Fig. 1.Particle velocity and displacement are set to be zero everywhere as initial conditions in an in¯nite and homogeneous dispersive medium.We assume that, given the ARF direction and extent, the displacements are polarized in the z-direction, and the derivatives with respect to and z are zero.Then, Eq. ( 1) can be written as IðrÞ, and using the notation of Gra®, 24 we can de¯ne the governing equation in a viscoelastic medium as where u z is the displacement of the shear wave in the z-direction; c is the velocity of the shear wave; F z ðrÞ is the applied body force proportional to the ARF beam IðrÞ; and T ðtÞ is the temporal application, which we will take as a rectangular function of duration a delayed a=2, rectð t a À 1 2 Þ.The Laplacian operator in cylindrical coordinates reduces to r 2 ¼ @ 2 @r 2 þ 1 r @ @r .Applying the Hankel transform H in the space, and the nonunitary angular frequency Fourier transform J in time to Eq. ( 2) in cylindrical coordinates yields where Û ð"; !Þ ¼ JfHfu z ðr; tÞgð"; tÞgð"; !Þ, " is the spatial frequency, ! is the temporal frequency, Then, Eq. ( 3) can be written as where k is the wavenumber.In a dispersive medium, the wavenumber is treated as a complex number k ¼ !c À i, where is the absorption coe±cient.Then, we apply the inverse Hankel transform to Eq. (4) as follows where J 0 is the Bessel function of the ¯rst kind.Then, by applying Baddour's theorem 25 and selecting the appropriate solution according to the Sommer¯eld radiation condition we obtain where H ð2Þ 0 ðxÞ is a Hankel function of the second kind.The form of Eq. ( 6) when the temporal application of the force F z ðrÞ is instantaneous (T ðtÞ ¼ ðtÞ) can be found in Parker and Baddour. 12Here, we seek a form when the force F z ðrÞ with a Gaussian beam pattern (Eq.( 1)) is applied for a transient time a in the media.Let where A 0 is the force amplitude, and is related to the curve width.Then, applying the spatial Hankel and temporal Fourier transforms to Eq. ( 7), we have where sincðxÞ ¼ sinðxÞ x .We are interested in ¯nding the analytic solution of Eq. ( 6) for particle velocity vðr; !Þ ¼ i!ûðr; !Þ.Then, using Eqs.( 6) and ( 8), we have which is a closed-form analytical solution in the r À !space that describes the cylindrical spreading of the pulsed wave, attenuation, and dispersion when it propagates through a viscoelastic medium.

Viscoelastic modeling of the medium
As discussed by Carstensen, 26 a variety of loss mechanisms with their own frequency-dependent solutions exist in viscoelastic materials, such as Kelvin-Voigt, Kelvin-Voigt fractional derivative, standard linear solid, and relaxation and hysteresis models.Given the uncertainties in appropriate choice of rheological model, we seek a general solution that is mechanism-independent.Therefore, we solely assume a limited shear wave frequency bandwidth of the excitation in which a Taylor series expansion can be used to express the frequencydependent behavior of the medium.Then, we introduce a ¯rst-order Taylor approximation of the frequency-dependent phase speed c and attenuation of the wave such that where c 0 ) c 1 !.As described by the most conventional loss mechanisms, 27 as !!0, c ! c 0 and !0.Then, we can consider 0 ¼ 0 for a conventional dispersive medium.Therefore, the complex frequency-dependent wavenumber will be described as and this form is used in Eq. ( 9) for k.For further notation simpli¯cation, we call k 0 ð!Þ ¼ !c 0 þc 1 j!j , and k 00 ð!Þ ¼ 1 j!j.In a linear and isotropic viscoelastic material, the complex shear modulus and G 00 ð!Þ is the loss moduli, can be obtained using the real and imaginary parts of the wavenumber in Eq. ( 11) 21 as where is the density of the material.In Figs.2(a) and 2(c), the magnitude of Eq. ( 9) is plotted in the r À !space for the viscoelastic (c 0 6 ¼ 0; c 1 6 ¼ 0, 1 6 ¼ 0), and pure elastic (c 0 6 ¼ 0; c 1 ¼ 0; 1 ¼ 0) cases, respectively.Estimates in r ¼ 0 were disregarded since H ð2Þ 0 ðkrÞ has a singularity at the origin.Moreover, particle velocity pro¯les versus r-axis, vðr; t 0 Þ ¼ J À1 fvðr; !Þgðt 0 Þ, at di®erent instants t 0 for both cases are shown in Figs.2(b) and 2(d).In the pure elastic case, the rapid decease in amplitude follows the asymptotic cylindrical spreading 1= ffiffiffiffiffi kr p ; while in the viscoelastic case, the low-pass smoothing of the particle velocity waveforms is evident, which diminishes the e®ectiveness of peak tracking approaches.

Shear to surface acoustic wave relationship
Surface acoustic waves (SAW) are produced when a vibration is generated at a solid-vacuum or solid-°uid interface. 28,29This work is intended for optical coherence elastography (OCE) applications.Therefore, we use an OCT system to acquire motion frames.Since the penetration depth of this system is typically some millimeters from the interface, SAW are more likely to be scanned than shear waves.If we consider a semi-in¯nite solid medium, the predominant SAW propagating in the solid-vacuum interface are Rayleigh waves. 28The relationship between shear wave and Rayleigh wave phase speed in a linear isotropic medium for a Poisson's ratio % 0:5 is given by 29 Moreover, Rayleigh waves from a point source follow cylindrical spreading 1= ffiffiffiffiffi kr p as described in Ref. 29, which is consistent with the asymptotic form of the Hankel term in Eq. ( 9), jH ð2Þ 0 ðkrÞj ffi ffiffiffiffiffiffi 2 kr

q
, for complex values of k. 30 Therefore, making the adjustment for phase speed as described in Eq. ( 13), the shear wave model of Eq. ( 9) is suitable for describing the Rayleigh wave propagation at a given depth z 0 .

Sample preparation and measurement
Two tissue-mimicking phantoms where used in experiments.A cylindrically-shaped custom shear wave viscoelastic phantom (model No. 16410001, CIRS, Virginia, USA) was used as the dispersive medium (5.4 cm in diameter Â 2.2 cm in height).
A cylindrically-shaped aqueous elastic phantom (Aqua°ex US del pad, Parker Laboratories INC., New Jersey, USA) was selected as the elastic medium (9 cm in diameter Â 2 cm in height).The frequency-dependent Young's modulus of each phantom was measured using a stress-relaxation by compression test.The measurement was conducted using a MTS Q-Test/5 Universal Testing Machine (MTS, Eden Prairie, Minnesota, USA) with a 5 N load cell using a compression rate of 0.5 mm/s, a strain value of 5%, and total measurement time of 600 s.The stress-time plots obtained by the machine were ¯tted to a three-parameter Kelvin-Voigt fractional derivative (KVFD) model 31 for the calculation of frequency-dependent complex Young's For all cases, the input force is a Gaussian pulse as described in Eq. ( 7), with ¼ 0:1 mm, and a ¼ 1 ms.The material properties are selected as c 0 ¼ 4 m/s, and c 1 ¼ 1 Ã 10 À11 m s =Hz. 1 ¼ 0:15 Np m =Hz, and 1 ¼ 0:0015 Np m =Hz for the viscoelastic and elastic cases, respectively.Particle velocity units are shown in radians referring to OCT Doppler phase-shift and can be transformed to m/s using Eq. ( 16).modulus given as where E 0 is the relaxed elastic constant, is the viscous parameter, and is the order of fractional derivative.Measurements were conducted in three samples of each phantom type in order to calculate the mean and standard deviation (error) for each predicted frequency-dependent result.Figure 3 shows the real and imaginary parts of Eq. ( 14) versus frequency for both phantoms.Then, for a nearly incompressible (Poisson's ratio close to 0.5), homogeneous, and isotropic medium, the complex phase velocity of the shear wave for each phantom with a material density of can be calculated 1 as , and the complex wavenumber of Eq. ( 11) can be expressed as where E 0 ¼ RealfE Ã g.For both phantoms, ¼ 1 g/ cm 3 .The three estimated parameters of the KVFD model are detailed in Table 1.We also included results of a compression test on the same phantoms in order to estimate the quasi-static Young's modulus Eð! % 0Þ, and the equivalent shear wave speed cð! % 0Þ.The equivalent frequency for the quasistatic compression test rate was found to be 0.25 Hz which is very low compared to the peak frequency band of the ARF shear waves.Predictions of the real part of the Young's modulus given by the KVFD model at f % 0:25 Hz (¯rst frequency sample in Fig. 3) are similar to results found in the quasistatic compression test which partially supports the validity of the model.Mechanical testing results in phantoms are fundamental for the validation of our proposed approach.Therefore, the selection of the correct rheological model which describes the relationship between the frequency-dependent prediction and the stress-relaxation test is required.Both elastic and viscoelastic phantoms were well-characterized using the KVFD model since a higher r-square ¯tting quality (0.999) was found as compared to other rheological models: Voigt, Kelvin-Voigt, Zener, Standard Linear Solid, and Standard Linear Solid Fractional Derivative models.
For the biological tissue test, we selected a porcine cornea.A fresh porcine eyeball was obtained from an abattoir (Joe's Meat Market, Ontario, NY, USA) with all experiments being performed within 12 h of their collection.For the scanning of the Fig. 3. Log-log of frequency-dependent complex Young's modulus for the elastic and viscoelastic phantoms obtained from mechanical measurements (MM) using a stress relaxation by compression test.!¼ 2f, and f ¼ frequency.Frequencydependent results are predictions of the Kelvin-Voigt fractional derivative model and do not represent independent measurements over the frequency range.

Experimental setup
The experimental setup is shown in Fig. 4(a).A 5 MHz confocal ultrasonic transducer (PIM7550-2inchFL, Dakota Ultrasonics, California, USA) with 5.01 cm of focal length was excited with a 1 ms sinusoidal tone of 5 MHz provided by a function generator (AFG320, Tektronix, USA).The generator was connected to a RF power ampli¯er (25A250, Ampli¯er Research, Souderton, PA, USA) in order to produce an approximate Gaussian, radially symmetric ( ¼ 0:338 mm) focused ARF push within the sample at the air-solid surface interface of the sample.The ultrasonic transducer was coupled to the sample with saline water as shown in Fig. 4(b).A phase-sensitive optical coherence tomography (PhS-OCT) system implemented with a swept source laser (HSL-2100-WR, Santec, Aichi, Japan) of a center wavelength of 1318 nm was used to acquire 3D motion frames of the sample within a region of interest (ROI) of 10 Â 10 mm in the xy-plane, and a maximum depth of 2.5 mm in the z-plane.The OCT acquisition and the excitation of the 5 MHz transducer were triggered by the computer controlling the entire process.The ARF push was focused at a certain (x 0 ; y 0 ) position in the sample's surface [Fig.4(b)], and it produced a localized out-of-plane vertical displacement which generated a cylindrically-shaped Rayleigh wave propagating within the ROI.The ARF vibration amplitude in the sample was adjusted to less than the maximum displacement (u z;max ) that can be detected by the OCT system without unwrapping the Doppler phase signal.For all cases, we used Eq. ( 12) to convert Rayleigh wave speed to shear wave speed.
The PhS-OCT characteristics include a full-width half-maximum (FWHM) bandwidth of 125 nm, and a light source frequency sweep rate of 20 kHz.The source power that entered the OCT interferometer was split by a 10/90 ¯ber coupler into the reference and sample arms, respectively.In the reference arm, a custom Fourier domain optical delay line was used for dispersion compensation.In the sample arm, a collimated light beam diameter of 6.7 mm at 1/e 2 was directed onto a test phantom by a focusing imaging lens (LSM05, Thorlabs Inc., NJ, USA), coupled with a galvanometer scanning mirror placed at the front focal plane of the imaging lens to achieve telecentric scanning.The back-scattered light from the sample was recombined with the light re°ected from the reference mirror with a 50/50 ¯ber coupler.The time-encoded spectral interference signal was detected by a balanced photo-detector (1817-FC, New CA, USA), and then digitized with a 500 Msamples/s, 12-bit-resolution analog-to-digital converter (ATS9350, AlazarTech, QC, CA).The maximum sensitivity of the system was measured to be 112 dB. 32The imaging depth of the system was measured to be 5 mm in air (-10 dB sensitive fall-o®).The optical lateral resolution was approximately 30 m, and the FWHM of the axial point spread function after dispersion compensation was 10 m.
The synchronized control of the galvanometer and the OCT data acquisition was conducted through a LabVIEW platform (National Instruments, Austin, TX, USA) connected to a workstation.The phase stability of the system was calculated as the standard deviation of the temporal °uctuations of the Doppler phase-shift (Á err ) while imaging a static structure. 33Results show Á err ¼ 4:6 mrad when using Loupas' algorithm. 34The displacement sensitivity is measured as the minimum detectable axial particle displacement (u z;min Þ.We found u z;min ¼ 0:358 nm.Finally, the maximum axial displacement supported by the system without unwrapping the phase-shift signal (Á max ¼ ) is u z;max ¼ 0:24 m.

Acquisition and processing approach
Due to speed limitations of the OCT acquisition system, a complete phase front of a single ARF excitation cannot be instantly acquired within a ROI.Therefore, a repeated excitation/acquisition triggered at any single spatial location within the ROI is conducted.This method is also called the M-B mode acquisition protocol, as described in Ref. 10. Methodologically, in a medium with refractive index n, the phase di®erence ÁðzÞ ¼ ðz; t 1 Þ À ðz; t 0 Þ at two consecutive instants t 0 and t 1 (t 0 < t 1 ), for a given (x o , y o ) position, is related to the particle velocity in the axial direction by where is the phase of the A-line signal acquired with the OCT system, n is the refractive index of the medium, 0 is the center wavelength of the laser, and Á M is the time sampling resolution.We select the Loupas' algorithm 34 for the accurate estimation of ÁðzÞ which increased the signal-to-noise ratio compared to other approaches as demonstrated by Zvietcovich et al. 10 The repeated acquisition and ARF excitation is possible since they are synchronized with the same cyclic trigger signal of 7 ms periodicity.The ARF tone is formed by 5000 cycles of a 5-MHz harmonic wave equivalent to 1 ms excitation.The M-B mode approach is used for generating x-axis space-time representations of particle velocity v z ðzÞ at a given z 0 and y 0 position.In this study, we constrained the analysis to the surface of the sample so that z 0 corresponds to the axial position of surface for any (x 0 , y 0 ) location.Each M-B acquisition spans 200 locations in the x-axis (10 mm), and M ¼ 100 time samples (5 ms) as shown in Fig. 5.Then, we repeat this process at each location y 0 in the y-axis for 200 positions (10 mm).
In total, we acquired a 3D matrix volume of 200 Â 200 Â 100 measurements of v z ðz 0 Þ in the x, y, and time axes, respectively, using the described protocol in order to cover the ROI of 10 Â 10 mm at z 0 , and 5 ms in the time frame.Each pro¯le cut of the 3D matrix in the xy-plane corresponds to a motion frame at a frame rate of 20 kHz (time resolution Á M ¼ 50 s).The total acquisition time was 4.6 min.

Inversion approach
In Sec.3.3, a 3D matrix representing particle velocity at the surface of the sample was acquired using the OCT system.Without loss of generality, we can call that data matrix v z 0 ðx; y; tÞ, where x and y represent spatial coordinates on the surface, and t is the time domain.We want to ¯t the analytical model of Eq. ( 9)vðr; !Þto the data in v z 0 ðx; y; tÞ, however, vðr; !Þ is de¯ned in the cylindrical r-axis.Therefore, taking a pro¯le cut of v z 0 ðx; y; tÞ along the xy-plane containing the center of the ARF excitation, and cropping the data to achieve a one-sided propagation matrix, we obtain v z 0 ðr; tÞ.
Taking vz 0 ðr; !Þ ¼ Jfv z 0 ðr; tÞg, and deriving A 0 in Eq. ( 9) by normalizing the data with the model for the maximum peak value, we can solve for the unknown parameters c 0 , c 1 , and 1 .Equation ( 17) is a nonlinear least squares problem and we employ the Levenberg-Marquardt method 35 for generating a solution.We solved Eq. ( 17) in four pro¯le cuts of v z 0 ðx; y; tÞ at angles ½0; =4; =2; 3=4 rad for the posterior calculation of average and standard deviation (error) of the estimated predictions.

Phantom experiments
Figure 6 shows particle velocity frames v z 0 ðx; y; ft 0 ; t 1 ; t 2 gÞ for three di®erent time instants extracted at the surface of each phantom as described in Sec.3.3.We de¯ned the initial time t ¼ 0 s at the falling edge of the ARF pulse.The ARF push location is found to be approximately at the center of the ROI.The cylindrical spreading of the wave when traveling out of the source is evident in both cases as shown in Fig. 6(e).The attenuation of the peak is more pronounced in the viscoelastic case when compared to the elastic case, as expected.
Pro¯le cuts in the x-direction centered to contain the ARF excitation origin were obtained from v z 0 ðx; y; tÞ in the viscoelastic and elastic cases as shown in Figs.7(a) and 7(b), respectively.A Gaussian bell shape is initially observed at the ¯rst instant as described in Sec.2.1.We ¯t F z ðrÞ from Eq. ( 7) to this curve in order to estimate , which is necessary for solving Eq. ( 17).Results are shown in Table 2.In addition, parameter a in the rectangular function of Eq. ( 7) is set to be a ¼ 1 ms as described in Sec.3.3.Figures 7(c the peak toward smaller values of frequency when r increases, which makes evident the presence of a frequency-dependent attenuation attributed to 1 .
In contrast, this behavior is not strong for the elastic case as expected.
The Levenberg-Marquardt method was applied to Eq. ( 17) for di®erent initial values of fc 0 ; c 1 ; 1 g, and we found convergent solutions for the viscoelastic and elastic cases, as reported in Table 2. Equation ( 9) is mesh plotted in the f À r space using the optimized parameters fc Ã 0 ; c Ã 1 ; Ã 1 g in Figs.7(c) and 7(d) for the viscoelastic and elastic cases, respectively.We also included in Figs.7(e) and 7(f), the pro¯le plots extracted at four di®erent r 0 positions from plots in Figs.7(c) and 7(d), respectively, for a better comparison of the experimental data and model predictions.
Given the results of mechanical testing in Table 1, we plotted the real and imaginary parts of the wavenumber as a function of frequency using Eq. ( 15) in Fig. 8, and we compared them with the components of the complex wavenumber obtained using results in Table 2 and Eq. ( 11).We found a good agreement of our model within the boundaries of the mechanical testing results for the viscoelastic and elastic cases.It is also evident that Ã 1 for the viscoelastic case is half an order of magnitude higher compared to the elastic one, which was expected.

Cornea experiment
Figure 9 shows particle velocity frames v z 0 ðx; y; ft 0 ; t 1 ; t 2 gÞ for three di®erent time instants extracted at the surface of a fresh porcine cornea.The ARF push location is found to be approximately at the center of the ROI. Figure 10(a) shows v z 0 ðr; tÞ obtained from a cropped pro¯le cut in the x-direction centered to contain the ARF excitation origin.We compensated r values for the curvature of the cornea.After applying the Fourier transform, jv z 0 ðr; !Þj is shown in Fig. 10  Eq. ( 17) was obtained similarly as in the phantom cases in order to ¯nd fc Ã 0 ; c Ã 1 ; Ã 1 g.Optimized parameters are reported in Table 2.

(b). A solution of
Corneal phase speed c Ã 0 is found to be consistent with other studies developed in porcine cornea. 5,15ttenuation parameter Ã 1 trends similarly to the viscoelastic case when compared to phantom results, which plays an important role in the pulse distortion and damping during wave propagation as demonstrated in Fig. 10.For further viscous parameter comparison with other corneal studies, a rheological model needs to be carefully selected, which will be a matter of future work.In this study, we are not taking into consideration the speci¯c boundary conditions of the cornea and, therefore, the presence of Lamb waves.Thus, speed results in the cornea reported here are solely valid for higher frequencies, when Lamb waves behave similarly to Rayleigh waves as described by Nenadic et al. 36 Concerning attenuation behavior, Lamb waves are spread and damped as described in Sec. 2 according to Schubert et al. 37 Therefore, attenuation results for corneal tissue shown in Table 2 are valid for the frequency-analyzed bandwidth.
Our current model is limited to the study of semi-in¯nite isotropic media which restricts its direct application to thin layers such as the cornea.However, this limitation can be resolved by extending the current general model for Lamb wave propagation in cylindrical coordinates and this will be the objective of future work.The acquisition time reported in this study is large and limits its application to ex vivo studies.The purpose of this research was to validate our method for the viscoelastic characterization of tissue.After validation, a faster and more stable OCT acquisition approach can be implemented for in vivo applications, as demonstrated by Singh et al. 38 and Song et al. 39 In addition, holographic techniques 40,41 demonstrated the fast acquisition of surface wave propagation in tissue, which o®ers an alternative to implement OCE methods in real time.

Conclusion
The application of an analytical general cylindrical wave propagation model for the viscoelastic characterization of dispersive media has been demonstrated.Experimental results in elastic and viscoelastic phantoms support the e®ectiveness of the approach when compared to mechanical testing results.Our proposed model of propagation takes into account (1) the initial force shape of the excitation pulse in the space-time ¯eld, (2) wave speed dispersion, (3) wave attenuation caused by material properties of the sample, (4) wave spreading caused by the outward cylindrical propagation of the wavefronts, and (5) the rheological-independent estimation of the dispersive medium.The model is versatile enough to be tuned with any type of input push-force produced by the desired excitation method.Moreover, in contrast to the majority of approaches that use only the frequency-dependent wave speed data for the calculation of viscoelastic parameters using a rheological model, our approach utilizes the wave attenuation data, which is fundamental for the accurate viscoelastic characterization of the sample without assuming any rheological model.This last advantage is signi¯cant for the study of tissue with unknown mechanical behavior.Finally, wave propagation in porcine cornea was analyzed and compensated for cylindrical spreading using this approach.Viscoelastic parameters of this tissue were calculated and reported.Future work will extend this research to a general Lamb wave model using di®erent rheological models of tissue for the accurate estimation of viscous parameters in porcine and human cornea, and their relationship with ocular diseases and treatments.

Fig. 1 .
Fig. 1.Schematic of a Gaussian-shaped pulse produced in a dispersive medium by an ultrasound ARF transducer.The center of the pulse is the origin for the Cartesian and cylindrical coordinate system.
cornea, the eyeball was immersed in a tank with saline solution at 25 C leaving the corneal tissue in contact with the air.The intra-ocular pressure (IOP) of the eyeball was measured during the experiment using a digital pressure gage (model DPGWB-04, Dwyer Instruments Inc., Michigan City, Indiana, USA).The average IOP reported during experiments was 9.3 mmHg.

Fig. 4 .
Fig. 4. Experimental setup.(a) PhS-OCT system implemented with a swept-source laser for motion detection.(b) Placement of the sample in a water tank, ARF US-transducer, and region of interest (ROI).Motion is produced in the surface of the sample.

Fig. 5 .
Fig. 5. M-B mode for acquisition protocol of motion signals using a PhS-OCT system.

Fig. 6 .
Figure6shows particle velocity frames v z 0 ðx; y; ft 0 ; t 1 ; t 2 gÞ for three di®erent time instants extracted at the surface of each phantom as described in Sec.3.3.We de¯ned the initial time t ¼ 0 s at the falling edge of the ARF pulse.The ARF push location is found to be approximately at the center of the ROI.The cylindrical spreading of the wave when traveling out of the source is evident in both cases as shown in Fig.6(e).The attenuation of the peak is more pronounced in the viscoelastic case when compared to the elastic case, as expected.Pro¯le cuts in the x-direction centered to contain the ARF excitation origin were obtained from v z 0 ðx; y; tÞ in the viscoelastic and elastic cases as shown in Figs.7(a) and 7(b), respectively.A Gaussian bell shape is initially observed at the ¯rst instant as described in Sec.2.1.We ¯t F z ðrÞ from Eq. (7) to this curve in order to estimate , which is necessary for solving Eq. (17).Results are shown in Table2.In addition, parameter a in the rectangular function of Eq. (7) is set to be a ¼ 1 ms as described in Sec.3.3.Figures7(c) and 7(d) show jv z 0 ðr; !Þj plots for the viscoelastic and elastic cases, respectively.The appearance of the sinc function in the frequency axis con¯rms the validity of T ð!Þ in Eq.(8).Moreover, the attenuation of the main lobe in jv z 0 ðr; !Þj for the viscoelastic case tends to move

Fig. 8 .Fig. 9 .Fig. 10 .
Fig. 8. Frequency-dependent wavenumber comparison between mechanical measurements (MM) results using the KVFD model and the application of our general propagation model (using optimized parameters fc Ã 0 ; c Ã 1 ; Ã 1 gÞ.!¼ 2f, and f ¼ frequency.The real and imaginary parts of wavenumber is plotted as a function of frequency for the viscoelastic and elastic phantoms.Frequency-dependent results are the predictions of the KVFD model and our general propagation model, and they do not represent independent measurements over the frequency range.

Table 1 .
Mechanical testing results in elastic and viscoelastic phantoms.Kelvin-Voigt fractional derivative parameters (left col.) and compression test results (right col.) are shown for both media.Quasi-static shear wave speed was calculated using cð! % 0Þ ¼