**Purpose**:
The purpose of this study was to demonstrate accurate measurement of corneal elastic moduli in vivo with noncontact and noninvasive optical coherence elastography.

**Methods**:
Elastic properties (in-plane Young's modulus, *E*, and both in-plane, μ, and out-of-plane, *G*, shear moduli) of rabbit cornea were quantified in vivo using noncontact dynamic acoustic micro-tapping optical coherence elastography (AµT-OCE). The intraocular pressure (IOP)-dependence of measured mechanical properties was explored in extracted whole globes following in vivo measurement. A nearly incompressible transverse isotropic (NITI) model was used to reconstruct moduli from AµT-OCE data. Independently, cornea elastic moduli were also measured ex vivo with traditional, destructive mechanical tests (tensile extensometry and shear rheometry).

**Results**:
Our study demonstrates strong anisotropy of corneal elasticity in rabbits. The in-plane Young's modulus, computed as *E* = 3μ, was in the range of 20 MPa to 44 MPa, whereas the out-of-plane shear modulus was in the range of 34 kPa to 261 kPa. Both pressure-dependent ex vivo OCE and destructive mechanical tests performed on the same samples within an hour of euthanasia strongly support the results of AµT-OCE measurements.

**Conclusions**:
Noncontact AµT-OCE can noninvasively quantify cornea anisotropic elastic properties in vivo.

**Translational Relevance**:
As optical coherence tomography (OCT) is broadly accepted in ophthalmology, these results suggest the potential for rapid translation of AµT-OCE into clinical practice. In addition, AµT-OCE can likely improve diagnostic criteria of ectatic corneal diseases, leading to early diagnosis, reduced complications, customized surgical treatment, and personalized biomechanical models of the eye.

^{1}

^{–}

^{3}compared to a tensile stress along the same plane.

^{4}

^{–}

^{7}To account for the anisotropic deformation response, a nearly incompressible transverse isotropic (NITI) material model based on collagen fiber macro-structure was developed. It was shown that the cornea must be described by at least two shear moduli (in-plane tensile and out-of-plane shear moduli),

^{8}a departure from commonly used single-modulus models. Additionally, it was shown that elastic waves in dynamic optical coherence elastography (OCE) studies of the cornea were accurately described by the NITI model. Noncontact acoustic micro-tapping OCE (AµT-OCE) was used to measure anisotropic mechanical properties in ex vivo porcine cornea, where the cornea's bounded structure forced guided wave propagation and enabled reconstruction of in-plane, μ (where

*E*= 3μ), and out-of-plane,

*G*, shear moduli by fitting experimentally obtained dispersion curves in the wavenumber-frequency domain to the theoretical (NITI) model.

^{8}

^{7}Specifically, the Young's modulus,

*E*, and out-of-plane shear modulus,

*G*, of fresh porcine whole globes inflated to controlled IOPs were measured with AµT-OCE before being cut into corneal buttons and tested using a parallel-plate rheometer and then stripped and tested via in-plane extension loading. Results were consistent with an order of magnitude difference between in-plane,

*E*, and out-of-plane,

*G*, elastic moduli for all testing methods. Although there were some differences in corneal condition during testing (such as curvature, boundary conditions, loading type, preconditioning, etc.), it was shown that moduli quantified from OCE data analyzed with the NITI model were accurate and provided a noncontact, nondestructive path to measure corneal anisotropic biomechanical properties.

^{9}

^{–}

^{11}Thus, to date, there are no studies demonstrating noncontact, noninvasive measurements of cornea shear elastic anisotropy in vivo.

*G*) and tensile extensometer (for independent determination of in-plane Young's modulus,

*E*). The results demonstrated that AµT-OCE can reliably generate and track elastic waves in rabbit cornea in vivo, enabling accurate reconstruction of corneal elasticity. In vivo results were compared with both controlled ex vivo and destructive mechanical tests that currently serve as the gold standard in biomechanical testing.

*G*, via rotational rheometry. Following rheometry, each sample was cut into strips for tensile testing of the in-plane Young's modulus,

*E*, using an extensometer. All tests were performed within 6 hours of euthanasia and all samples were transported using cloths dampened with BSS.

^{12}

^{,}

^{13}) contributes to a stress-strain relationship that can be approximated using four independent elastic constants λ, δ,

*G*, and μ (assuming small-strain deformation).

^{7}

^{,}

^{14}The elasticity matrix in the cornea (in Voigt notation) can thus be described as:

_{ij}denotes engineering stress, ε

_{ij}denotes engineering strain, τ

_{ij}denotes shear stress, γ

_{ij}= 2ε

_{ij}denotes shear strain, and the subscripts

*x*,

*y*, and

*z*denote standard Cartesian coordinates. Because the cornea is nearly incompressible, the longitudinal modulus λ does not define deformation, and the corneal strain response to an applied stress can be fully defined by the out-of-plane and in-plane shear modulus (

*G*and μ, respectively), and an additional term, δ.

^{7}However, corneal structure constrains the latter to the range −2μ < δ < 0, which restricts the Young's modulus to the range:

*E*suggests that μ and

_{T}*G*can provide very close approximations of full cornea deformations assuming small deformation. Whereas the exact relationship between the in-plane Young's modulus (

*E*) and δ has not yet been determined, we assume tensile isotropy here so that

_{T}*E*=

*E*≅ 3μ. The approximation of δ = 0 assumes corneal tensile isotropy and likely leads to a slight overestimation of

_{T}*E*.

_{T}^{7}Briefly, elastic waves were generated using a cylindrically focused air-coupled piezoelectric transducer (AµT) driven with a 100 µs-long chirped (1 MHz to 1.1 MHz) waveform providing a temporally localized and spatially focused acoustic “push.” The resulting elastic wave was measured using a phase-stable Michelson-type fiber-optic interferometer where a broadband superluminescent diode (SLD1018P, Thorlabs, Newton, NJ) with central wavelength 1310 nm (45 nm full-width-half-maximum bandwidth) was coupled into polarization maintaining fibers and components for depth encoded (1.5 mm effective imaging range) OCT imaging and motion detection.

^{15}). Each MB-scan provided a 3D data set with 1024 depth × 256 lateral × 256 temporal dimensions. The scan rate of the system (determined by the line-scan camera) was 90 kHz, corresponding with a total scan time of around 750 milliseconds. For ex vivo scans at low IOP, 512 temporal scans were taken to allow the elastic wave to propagate across the full 10 mm field of view.

^{16}where the measured displacement sensitivity of the system was approximately 1 nm in water. The log-compressed real-part of the OCT signal was used to reconstruct corneal structure, and the surface of the cornea was detected using an automatic segmentation algorithm. The vibration velocity along the surface of the cornea was determined using a weighted-average (one half of a Gaussian window, HWHM = 90 µm, weight decreasing with depth) along the anterior 183 µm of the cornea, providing raw space-time (

*x-t*) maps of the vertical displacement from propagating guided elastic waves detected along the air-tissue boundary.

*x-t*data where vibration frequencies below 50 Hz and above 4 kHz were removed. Additionally, randomly propagating diffuse wavefields

^{17}were limited by a directional filter, where a two-dimensional Fast Fourier Transform (2D FFT) was performed on the

*x-t*plot and a mask applied to the second and fourth quadrant of the

*k-f*space, followed by an inverse 2D FFT.

^{18}Finally, unwanted reflections and forward propagating waveforms from sidelobes in the acoustic excitation, or remaining diffuse propagations, were removed by applying a moving temporal window centered on the peak of the vertical velocity in the

*x-t*plot. The moving temporal window utilized a super-Gaussian (

*SG*) function (Equation 3) that followed the maximum vibration velocity of the wavefield \(t_m^{wf}( x )\) at each discrete position

*x*:

_{t}= 0.5 ms. The resulting surface

*x-t*information was used to reconstruct elastic moduli assuming a NITI material.

*x-t*) plots (Fig. 2b) were subject to a 2D FFT to display the waveform in wavenumber-frequency (

*k-f*) space. An inversion method based on the solution to guided elastic waves in a NITI material quantified out-of-plane shear modulus,

*G*, as well as in-plane Young's modulus,

*E*, (assuming

*E*=

_{L}*E*).

_{T}^{7}

^{,}

^{8}As the resolution limit (determined in guided materials by the wavelength and corneal thickness

^{19}) was approximately 1 mm, the corneal surface area available for imaging (approximately 10 mm) was sufficient for accurate reconstruction of elastic moduli.

_{0}-mode is expected in the cornea for the frequency ranges produced by AµT, the A

_{0}-mode solution (detailed previously

^{8}) was solved in

*k-f*space for a broad range of input shear moduli (

*G*and μ) with a fixed corneal thickness. The central cornea thickness was measured in each scan using the OCT image assuming a refractive index of 1.38.

^{20}

*G*, respectively), where the best-fit theoretical dispersion relation was performed by maximizing the following objective function using simplex optimization (

*fminsearch*, MATLAB, MathWorks, Natick, MA):

*k-f*space is \(\hat v\), the A

_{0}mode solution for a NITI material is decribed by

*w*(

*f*,

*k*; μ,

*G*), and the nearly incompressible assumption is defined by β (set to 1 based on an L-curve analysis).

^{2}The value of λ was updated at each iteration according to \(\lambda = \rho c_L^2 - 2\mu .{\rm{\;}}\)The corneal density was assumed to be

*ρ*= 1000 kg/m

^{3}and corneal longitudinal wave speed

*c*= 1540 m/s. The cornea was assumed bounded from below by water with a density of 1000 kg/m

_{L}^{3}and longitudinal wave speed of 1540 m/s. To avoid convergence to a local (as opposed to global) maximum in Equation 4, five independent fits were performed, with quasi-random (within a reasonable range of expected values) initial values of

*G*

_{0}and µ

_{0}. The final output of μ and

*G*were set to those corresponding to the highest value in Equation 4. An example of the resulting best-fit A

_{0}-mode superimposed on top of the

*k-f*spectrum in a rabbit cornea is shown in Figure 2c.

_{0}mode that most closely matches experimental data, it does not provide information on the fit quality. Although the best-fit line generally follows the spectral maxima (see Fig. 2c), the degree to which the theoretical A

_{0}mode follows the measured spectral maxima is not directly determined by the fitting regime. By normalizing the best-fit line to the unconstrained global maximum of the spectral peaks at each spectral bin in

*k-f*space (Φ

_{max}), the quality of the fit can be estimated:

*g*

_{NITI}describes the fraction of the maximum mode energy captured by the best-fit A

_{0}dispersion curve. A value of

*g*

_{NITI}= 1 would indicate that the experimental dispersion curve has the exact match with the NITI model (Φ

_{NITI }captures all the measured mode's maximum energy). Because low values of

*g*

_{NITI}correspond with a failure to converge to a solution that captures most of the mode energy, it can be used to exclude unreliable scans (due to misalignment, rabbit motion, model inaccuracies, etc.). The output values of

*G*and μ were considered inaccurate for any fit that corresponded to

*g*

_{NITI}below 0.87, and 0.92, respectively. Details on determining exclusion parameters can be found in Supplementary Methods.

*G*and μ were varied independently around the optimum values of Φ

_{NITI}. For each combination of

*G*and μ, Equation 4 was used to calculate Φ (μ,

*G*), which was then used to determine:

*g*

_{NITI}) was determined independently for at least five repeat scans (in both in vivo and ex vivo data sets). The standard deviation of

*g*

_{NITI}was then used to determine the cutoff value in ψ for each set of scans.

*g*

_{NITI}(

*g*

_{NITI}= ψ

_{max}) was 0.98 with a standard deviation of 0.003 (or 0.3%) across 5 independent scans. The variation in

*g*

_{NITI}determined the uncertainty interval using the corresponding

*G*and μ values for ψ at 0.3% below the maximum (i.e. ψ = 0.977 for

*G*and μ). This method produced uneven error bars due to the shape of ψ. Note that the absolute value of

*g*

_{NITI}provides an estimate of model error. When

*g*

_{NITI}is reduced, the shape of ψ widens for both moduli and, therefore, model uncertainty increases.

*G*and μ, and their respective uncertainty ranges. The mean value of

*G*and μ for all repeat fits was considered for each cornea (at each IOP). Uncertainty ranges were calculated as root mean squared (RMS) values of lower and upper uncertainty limits in independent scans. The process provided a value for the OCE-measured out-of-plane shear modulus,

*G*, and in-plane Young's modulus,

*E*(assuming

*E*= 3μ), as well as error bars associated with the uncertainty of reconstruction. Uncertainty intervals for each independent OCE measured value can be found in the Supplementary Material.

*G*. A rheometer (Anton Paar MCR 301 Physica) assessed the frequency-dependent shear behavior (storage,

*G*′(ω), and loss,

*G*″(

*ω*)) of corneal buttons over a range of 0.16 Hz to 16 Hz. A 5 N compressive preload was applied and the peak shear strain was approximately 0.1%. The test was performed twice, and the mean of each run was used for the final value.

*E*, was defined as the tangential slope of the stress-strain curve. Extension testing provided a value for the strain-dependent Young's modulus,

*E*, up to 10% strain.

^{21}

^{,}

^{22}a correction factor was applied to all in vivo tonometry measurements to facilitate comparison between in vivo and ex vivo measurements. The correction factor in this study was determined using a direct comparison of Tono-Pen measurements made on five ex vivo whole globe samples inflated to known IOPs.

*IOP*(corresponding to a coarse estimate of the “actual” IOP) was found for each Tono-Pen determined value (

_{Control}*IOP*). The result of this test suggests that a Tono-Pen measured value of 10 mmHg, for example, corresponds to an actual value of 13 mmHg (with uncertainty associated with the standard deviation of the values at each pressure between approximately 6 mmHg and 21 mmHg). The mean difference between the Tono-Pen and actual values was approximately 6.5 mmHg across all pressures for the sample size used. As such, error bars corresponding to ±6.5 mmHg are included in the adjusted in vivo pressure values displayed in Figure 4. All error bars were cut off at 0 mmHg for display.

_{Tonopen}*G*(blue squares) and

*E*(red squares) are plotted at the corrected pressure (horizontal error bars are associated with the range of possible “actual” pressure values). Ex vivo moduli (triangles) were measured and displayed at the actual IOP value. The vertical error bars in OCE measured values correspond with the standard deviation of the mean in each cluster at each IOP. Note that of the 10 samples tested, one had a fit quality below the exclusion criteria and was omitted. Because the goodness of fit generally decreased with increasing IOP for both moduli, there are fewer data-points at pressures greater than 17 mmHg (Supplementary Methods).

*G*was 34 kPa to 261 kPa, and for

*E*was 20 MPa to 44 MPa. For both in vivo and ex vivo whole globes measured via OCE, the stiffness generally increased with increasing pressure. In ex vivo whole globes,

*G*increased from 31 kPa (±15 kPa) in the IOP range from 3 mmHg to 5 mmHg, to 98 kPa (±65 kPa) in the IOP range from 11 mmHg to 13 mmHg, and

*E*increased for the same IOP ranges, respectively, from 27 MPa (±9 MPa) to 47 MPa (±13 MPa). The values are not reported at higher pressures due to decreased sample size with reliable quality of fit. The mean and standard deviation of rheometry measured values of

*G*′(ω) (see Fig. 5a) at 16 Hz was 75 kPa (±43 kPa). Note that shear rheometry was performed over a lower range of frequencies than OCE, which presumably should lead to a slightly lower estimate of

*G*, as shown. The tensile modulus,

*E*, increased with strain from 2.8 MPa (±1.1 MPa) at 1% to 32 MPa (±20 MPa) at 10%. The mean and standard deviation of the high strain value is displayed in Figure 5c. Corresponding frequency-dependent values of

*G*′(ω) and

*G*″(

*ω*)) (measured in all corneal buttons with shear rheometry), as well as the strain-dependent value of the Young's Modulus,

*E*, along with OCE measured values, can be found in Supplementary Methods for each sample included in the analysis. As can be seen in Figure 5, the range of values measured with all methods (for comparable IOP in OCE) were in close agreement.

^{23}have shown (in their case for anesthetized rodents) that the frequency of respiratory motion and the heartbeat do no exceed 5 Hz and that the associated root mean square (RMS) motion amplitude of both cornea and retina does not exceed 1.5 µm (sub-pixel in our case). A typical time for blood pressure change-induced axial motion is about 1 second. On the other hand, the time between each M-scan in our study is about 3 ms (approximately 256/90e3). As such, we can infer that recordings are not highly influenced by axial motion of the eye. In addition, the low frequency filtering induced by windowing the

*x-t*plot removes the residual low frequency component. Finally, because the top-surface displacement is in fact an averaged and weighted displacement from the first 200 µm below the surface of the cornea, small motion effects can be disregarded. Bulk motion in awake animal and human models would likely require faster scan times or surface correction algorithms to account for micrometer and millimeter scale movements.

^{19}

^{,}

^{24}Dense-scanning of propagating elastic waves creates high-resolution images of local group velocity,

^{25}but these images do not map anisotropic elastic moduli at the same resolution. In reconstructing quantitative elastic moduli using guided waves, the resolution is determined by (i) the elastic wavelength of the generated mechanical wave

^{24}and (ii) the corneal thickness.

^{19}For AµT excitation, the theoretical elastic wavelengths generated can yield spatially resolved elastic moduli maps with mm-scale resolution. Demonstration of local anisotropic elastic moduli remains an area of future focus.

*E*, is on average multiple orders of magnitude larger than the out-of-plane shear modulus,

*G*. The NITI model used in this work introduced two independent shear moduli, separating in-plane from out-of-plane moduli. It helps to explain the order(s) of magnitude difference in corneal stiffness estimates extracted from shear- and tensile-based mechanical measurements. The results of this study show that AµT-OCE can be used to quantify anisotropic mechanical properties in living corneal tissue. As reported above, anisotropic mechanical moduli in rabbit corneas measured via OCE under both in vivo and ex vivo conditions were in close agreement with the values measured via destructive ex vivo testing.

*G*, only. As such, the shear modulus, μ, cannot be determined in a semi-infinite or bulk material. However, guided elastic wave propagation in the cornea plane also creates a weak dependence on μ through boundary conditions.

^{8}As shown here and in our recent work,

^{7}

^{,}

^{8}a change in μ results in a small change in the low-kHz region of the A

_{0}-mode dispersion. Thus, a deviation of the NITI model from the actual cornea may further worsen the goodness of fit and result in a large confidence interval for reconstructed modulus μ. A recent study,

^{26}for example, suggests that the wave dispersion can also be affected by corneal prestress induced by IOP and its in-plane boundary conditions, which are not taken into account in our study. Although the reconstruction of

*G*seems to be quite stable to small inaccuracies in the model used here, this model can be further refined for better quantification of μ. Therefore, the reconstruction of μ from experimental data is less accurate compared to the reconstruction of

*G*. Under these circumstances, it is important to evaluate whether the measured elastic waveform can be appropriately described using the NITI model. To determine whether μ can be accurately determined assuming a NITI model, a standardized method to determine data inclusion and model appropriateness is very important. In this work, we developed a method to quantify data quality based on what we refer to as the goodness of fit (between measured and theoretical dispersion curves). Using a statistical analysis of experimental data, we have shown that the reconstruction of

*G*and μ can be considered correct for any fit when

*g*

_{NITI}is above 0.87, and 0.92, respectively. Details on the determination of exclusion parameters can be found in Supplementary Methods.

*E*=

_{T}*E*= 3μ). Most likely, this assumption requires further refinement of the NITI model, particularly in the definition of the proportionality coefficient between

_{L}*E*and μ. Because corneal shear anisotropy is extremely strong, we expect this relationship to be closer to

_{T}*E*= 2μ. Measuring the exact relationship on a large population of cornea samples will be a part of our future studies. Here, however, we prefer to keep the relationship

_{T}*E*=

*E*=

_{T}*E*= 3μ until it is accurately measured.

_{L}*G*and

*E*increased with IOP for in vivo and ex vivo measurements. Additionally, the stress-strain curves in Supplementary Materials show that tensile Young's Modulus

*E*is strain dependent. Clearly, the pre-strain condition will play a role in measured corneal stiffness. As an aside, it has also been shown that increased stress plays a role in the propagation behavior of elastic waves, even for linear elastic materials.

^{27}

^{,}

^{28}Different techniques have shown the possibility of disentangling the effect of IOP induced pre-stress from the change of stiffness for linearly isotropic solids.

^{26}

^{,}

^{28}

^{,}

^{29}Decoupling nonlinear mechanical responses from strain-induced changes in wave speed for anisotropic materials remains an area of future interest.

*E*and

*G*are real moduli (i.e. the loss modulus related to tissue viscosity is assumed to be negligibly small). As can be seen in the rheometry plots shown in the Supplementary Methods,

*G*has a non-zero loss modulus and the frequency response is consistent with a measurable viscosity. Methods to estimate corneal viscosity (such as by tracking elastic wave energy) are currently under investigation. However, we notice here that cornea viscosity mostly affects the higher frequency range of the dispersion relationship for the

*A*

_{0}mode,

^{30}

^{,}

^{31}and thus has very little effect on values estimated from the real part of moduli. In other words, the cornea loss modulus is not measured in this study, but introducing complex elastic moduli into the NITI model would produce small changes in the reconstructed elastic (real) part of moduli.

^{21}

^{,}

^{22}the mean tonometry-measured IOP generally underestimated the actual IOP at high pressures, with very large individual variability between samples. One limitation in the present study was that only five samples were used to create the correction factor. Due to large individual variations, a more accurate correction factor would require a much larger population sample size. Because contact tonometers (such as the Tono-Pen) use a simple model to estimate internal pressure based on displacement of the cornea, the final value it measures is inherently linked to not only IOP, but also to cornea biomechanics. Independent measurements of corneal elastic moduli may help to develop better mechanical models that can potentially lead to more accurate estimates of IOP.

*G*, at different frequencies. Whereas rheometry operates in the Hz-range, OCE operates in the multiple-kHz regime. As can be seen in the individual plots shown in Supplementary Methods, the apparent modulus,

*G*, increases with increasing frequency, suggesting that the OCE results reported in Figure 5 would be higher than the expected shear modulus of the cornea under lower frequency shear strain, as well as under typical biological shear strain rates (due to eye rubbing, for example). However, the results obtained with all methods were still in relatively close agreement, which supports the NITI model for reconstruction of cornea moduli from OCE data.

^{32}rabbit cornea may not be accurately described by the NITI model at high IOP (above 15 mmHg), and an orthotropic model (and, therefore, a larger number of elastic moduli) may be required.

**Data Availability:**The authors declare that all data from this study are available within the article and its Supplementary Information. Raw data for the individual measurements are available upon reasonable request.

**M.A. Kirby**, None;

**G. Regnault**, None;

**I. Pelivanov**, None;

**M. O'Donnell**, None;

**R.K. Wang**, None;

**T.T. Shen**, None

*Investig Ophthalmol Vis Sci*. 2012; 53: 873–880. [CrossRef]

*J Biomech*. 2014; 47: 723–728. [CrossRef] [PubMed]

*Investig Ophthalmol Vis Sci*. 2014; 55: 7919–7924. [CrossRef]

*J Biomech*. 2001; 34: 533–537. [CrossRef] [PubMed]

*J Cataract Refract Surg*. 2003; 29: 1780–1785. [CrossRef] [PubMed]

*J Biomech Eng*. 2012; 134: 031003. [CrossRef] [PubMed]

*Ophthalmol Sci*. 2021; 1(4): 100058. [CrossRef] [PubMed]

*Sci Rep*. 2020; 10: 12983. [CrossRef] [PubMed]

*Biomed Opt Express.*2019; 10: 6272–6285. [CrossRef] [PubMed]

*Biomed Opt Express.*2022; 13: 2644–2654. [CrossRef] [PubMed]

*Prog Retin Eye Res*. 2015; 49: 1–16. [CrossRef] [PubMed]

*Exp Eye Res*. 2015; 133: 81–99. [CrossRef] [PubMed]

*IEEE Trans Ultrason Ferroelectr Freq Control*. 2004; 51: 322–328. [CrossRef] [PubMed]

*J Biomed Opt*. 2013; 18: 121509. [PubMed]

*Appl Phys Lett*. 2007; 90: 164105. [CrossRef]

*J Biomed Opt.*2017; 21: 126013. [CrossRef]

*IEEE Transacctons Ultrason Ferroelectr Freq Control*. 2011; 58: 2032–2035. [CrossRef]

*Biomed Opt Express*. 2022; 13: 4851–4869. [CrossRef] [PubMed]

*Optics Letters*. 2004; 29(1): 83–85. [CrossRef] [PubMed]

*Sci Rep*. 2016; 6: 35187. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci.*2005; 46(7): 2419–2423. [CrossRef] [PubMed]

*Exp Eye Res*. 2010; 91(1): 63–68. [CrossRef] [PubMed]

*J Biomed Opt*. 2019; 24: 1–16. [CrossRef]

*Sci Rep*. 2016; 6: 38967. [CrossRef] [PubMed]

*J Mech Behav Biomed Mater*. 2022; 128: 105100. [CrossRef] [PubMed]

*Commun Phys*. 2022; 5: 1–7.

*J Eng Sci Med Diagn Ther*. 2023; 6(1): 011006. [PubMed]

*J Acoust Soc Amer*. 2022; 151(4): 2403–2413. [CrossRef]

*Opt Express*. 2019; 27: 16635. [CrossRef] [PubMed]

*Sci Rep*. 2020; 10: 17366. [CrossRef] [PubMed]

*Opt Elastography Tissue Biomech III*. 2016; 9710: 97100T.