Nearly-incompressible transverse isotropy (NITI) of cornea elasticity: model and experiments with acoustic micro-tapping OCE

The cornea provides the largest refractive power for the human visual system. Its stiffness, along with intraocular pressure (IOP), are linked to several pathologies, including keratoconus and glaucoma. Although mechanical tests can quantify corneal elasticity ex vivo, they cannot be used clinically. Dynamic optical coherence elastography (OCE), which launches and tracks shear waves to estimate stiffness, provides an attractive non-contact probe of corneal elasticity. To date, however, OCE studies report corneal moduli around tens of kPa, orders-of-magnitude less than those (few MPa) obtained by tensile/inflation testing. This large discrepancy impedes OCE’s clinical adoption. Based on corneal microstructure, we introduce and fully characterize a nearly-incompressible transversely isotropic (NITI) model depicting corneal biomechanics. We show that the cornea must be described by at least two shear moduli, contrary to current single-modulus models, decoupling tensile and shear responses. We measure both as a function of IOP in ex vivo porcine cornea, obtaining values consistent with both tensile and shear tests. At pressures above 30 mmHg, the model begins to fail, consistent with non-linear changes in cornea at high IOP.

. Relationships between mechanical forces in the cornea and various mechanical testing methods. Mechanical forces in corneal tissue govern focusing power in (a) healthy and (b) pathological tissue (e.g. keratoconus). (c) Elastic shear waves generated in the cornea with a transient force are tracked using optical coherence tomography to determine corneal mechanical properties without contact. (d) Summary of biomechanical test methods and corneal Young's and shear moduli from the literature. Reported moduli vary by up to four orders of magnitude depending on the loading technique and test configuration, suggesting that isotropic models cannot accurately describe corneal mechanics. In particular, tensile and inflation tests generally agree, as do shear and transient tests. Values listed are for fresh samples (within 2 weeks), and best effort was taken to report moduli in the low-strain/preload region where they generally are near their lowest value.

Results nearly incompressible transversely isotropic (niti) model of cornea. Corneal stroma contains
hundreds of vertically stacked collagen lamellae, each 0.2-2.5 μm thick 34 with a preferred collagen orientation. While some in-plane anisotropy has been reported 28,[35][36][37][38] , its magnitude at low intraocular pressure (IOP) suggests that macroscopic behavior can be treated as isotropic in-plane 22,38,39 . Recent second harmonic generation imaging studies show that lamellar orientations are more random than previously suggested 32,40,41 , further supporting the assumption of in-plane isotropy. One interpretation of this structure, based on a fiber-composite model 42 , is that collagen fiber mechanical properties govern in-plane behavior, while those of the connective tissue matrix govern out-of-plane behavior.
A transversely isotropic (TI) model is the most appropriate given an isotropy plane. It contains five independent elastic constants (C 11 , C 12 , C 13 , C 33 , C 44 ) rather than the two (Lamé constants) of isotropic materials. To simplify notation, we adopt the shorthand: C 11 = λ + 2μ, C 12 = λ, and C 44 = G. We also assume that the cornea, like most soft tissue, is nearly-incompressible. That is, the medium's internal pressure remains finite as λ → ∞ and the dilatation approaches zero. When applied to a transversely isotropic solid, this assumption leads to a set of conditions on the longitudinal terms of the stiffness matrix (C 11 , C 12 , C 13 , C 33 ). These conditions do not reduce the number of independent constants, but they do provide limiting relationships that define the nearlyincompressible condition 43 .
Of note, any transversely isotropic solid for which C 13 and C 33 are asymptotically equal to λ will behave as a nearly-incompressible solid, provided the stiffness matrix remains invertible (Supplementary Note 2). This last condition must be treated carefully. When the longitudinal part of the stiffness matrix is isotropic (C 13 = C 12 and C 11 = C 22 = C 33 ), the stiffness matrix is always invertible. This condition corresponds to weak anisotropy, an assumption supported by Brillouin microscopy measurements 44 . We note that more exact estimates of C 13 and C 33 can be independently measured with Brillouin microscopy; however, these terms contribute little to predicting mechanical (shear, surface, guided) wave behavior in nearly-incompressible solids, i.e. have very little effect on shear moduli inversion from dynamic elastography measurements (Supplementary Note 4).
Thus, we take C 13 = λ and C 33 = λ + 2μ and finally define the stiffness tensor of a nearly-incompressible transversely isotropic (NITI) material as: www.nature.com/scientificreports/ The constants λ and μ mimic those in an isotropic solid, with λ enforcing incompressibility and μ defining in-plane shear, tensile, and compressive behavior. Similar to an isotropic material, the Young's modulus is simply related to μ, E = 3μ. An additional shear constant G governs out-of-plane shear and is completely decoupled from E.
Uniaxial tensile and inflation tests yield Young's modulus estimates related only to μ. However, shear torsional tests depend only on G. In addition, the speed of vertically-polarized bulk shear waves is a function of G. This decoupling of normal and shear deformations helps explain the discrepancy between tensile/inflation test modulus estimates (on the order of MPa) and shear/transient estimates (on the order of kPa). Supplementary Note 2 shows how to obtain NITI parameters μ and G from tensile, inflation, shear, and transient mechanical tests.
Wave behavior in a bulk NITI medium. Unlike isotropic materials, which support only two bulk waves (longitudinal and shear), transversely isotropic materials support three-quasi-longitudinal, quasi-shear, and shear. Soft tissue is nearly-incompressible (λ ≫ μ), implying that the quasi-longitudinal wave speed is nearly constant over all directions of a NITI medium. The anisotropy of longitudinal wave propagation in cornea has been recently shown to be very weak, with angular variations in the quasi-longitudinal wave speed < 5% 44 . In contrast, quasi-shear and shear wave speeds vary dramatically with angle and depend on both G and μ. This has important implications for OCE measurements.
Many dynamic OCE methods track mechanical waves propagating along the air-cornea interface, ignoring liquid loading on the cornea's posterior surface 28,36,45,46 . In other words, the cornea is considered semi-infinite with a simple Rayleigh wave propagating along the surface. The in-plane Rayleigh wave speed can be obtained numerically from the Stroh formalism [47][48][49][50][51] . For materials with G < μ, such as we expect for cornea, the Rayleigh wave speed varies from c R = √ G/ρ in the highly anisotropic limit (G ≪ μ) to c R = 0.9553 √ G/ρ in the isotropic limit (G = μ) (Supplementary Note 3 , Fig. S3.2). Even for varying degrees of anisotropy, Rayleigh wave speed is primarily governed by G and only slightly influenced by μ. This leads to an important conclusion. If the cornea could be considered a semi-infinite NITI medium, it would be extremely difficult to determine μ from OCE measurements. For reference, Supplementary Note 3 derives bulk and Rayleigh wave speeds in a NITI medium.
In reality, the cornea's thickness is typically on the order of the wavelength of propagating mechanical waves considered in dynamic OCE 22,26,52 . As such, corneal thickness and boundary conditions lead to guided waves. As shown below, guided waves provide additional information to help extract both G and μ from dynamic OCE measurements.
Guided wave behavior in a bounded niti layer. The cornea's bounded geometry produces dispersive guided waves, which must be analyzed in frequency-wavenumber (ω-k) space to quantify elasticity. Partial wave analysis assuming the cornea as a flat isotropic layer bounded by air and water (Fig. 2a) leads to a secular equation describing guided modes 23,26 . Here, we introduce the dispersion relation for a NITI layer bounded by air and water, derived from partial wave solutions to the elastic wave equations satisfying corneal boundary conditions. Supplementary Note 4 provides a full derivation, and functions solving this dispersion relation are provided in Supplementary Software.
While the solution includes an infinite number of quasi-antisymmetric (A) and quasi-symmetric modes (S), only the first two (A 0 and S 0 ) are typically captured in OCE. Figure 2 shows some examples for NITI materials with varying levels of anisotropy (G = 20 kPa, h = 0.55 mm, varying μ). As μ increases, the S 0 mode propagates with higher phase velocity at low frequencies (Fig. 2c). In contrast, the high frequency asymptote of the A 0 mode is primarily governed by G, with μ controlling both the rate at which the A 0 mode approaches its asymptote and weakly affecting the asymptotic value (< 5%) (Fig. 2d). Thus, the cornea's bounded structure provides a potential avenue to determine both G and μ using OCE. finite element model of guided waves in a bounded niti medium. Although an analytical solution exists, it is not valid near the mechanical wave excitation source (near field) and does not describe how energy is distributed among guided modes in a NITI layer. To address this, we developed a finite element model using OnScale (OnScale, Redwood City, CA). It contains a NITI layer bounded above by air and below by water with dimensions similar to the cornea. Elastic waves were excited with a spatio-temporally short pressure applied to the top surface, mimicking AμT. Supplementary Note 5 presents a detailed model description, and the OnScale input file is available in Supplementary Software.
We examined temporal changes in the surface velocity field over a range of lateral positions. Unitless surface vibrations illustrate the spatio-temporal behavior of the surface wave field and mimic an OCE measurement (Fig. 3a,b). This is referred to as an XT plot. A two-dimensional Fourier transform was then applied to analyze guided wave dispersion. Figure 3 shows XT plots and corresponding 2D Fourier spectra for increasing μ. In the isotropic limit (μ = G), the XT plot (Fig. 3c) clearly shows multiple interacting guided waves, with A 0 and S 0 modes both visible in the spectrum (Fig. 3g). As anisotropy increases (μ > G), XT plots ( Fig. 3d-f) change dramatically. Interference between A 0 and S 0 modes decreases at μ = 2G and is nearly absent at μ = 5G. Furthermore, as μ increases, the shapes of A 0 and S 0 modes shift. Energy in the S 0 mode also decreases, disappearing almost www.nature.com/scientificreports/ entirely for μ ≥ 5G (Fig. 3i,j). These results strongly suggest that if the NITI model is valid, then only one mode, A 0 , should be visible in cornea. This is in contrast to an isotropic layer of the same thickness, where multiple high order modes can be observed.

Elastic modulus estimates with AμT-driven OCE.
A spectral-domain OCT system with a 46.5 kHz effective frame rate, as detailed previously 53 , tracked guided waves in isotropic polyvinyl alcohol (PVA) cryogels and porcine cornea, providing experimental measurements to compare with theoretical predictions. A cylindrically-focused 1 MHz air-coupled ultrasound transducer (AμT) provided a spatio-temporally sharp push to the surface of a thin isotropic PVA phantom or freshly excised porcine cornea (n = 6) with IOP incrementally increasing from 5 to 35 mmHg, generating mechanical waves with bandwidths up to 4 kHz ("Methods" section). Figure 4 compares OCE-measured surface velocity fields obtained for an isotropic PVA phantom (Fig. 4a) and porcine cornea at an IOP of 5 mmHg (Fig. 4d). Guided waves are apparent in the isotropic phantom (Fig. 4b).
The 2D spectrum of this XT wave field clearly has two guided modes, and fitting an isotropic dispersion relation 26 yields a shear modulus estimate, μ PVA = 14.4 kPa (Fig. 4c, yellow curves).
Porcine cornea displays very different behavior, with the wave energy concentrated in a single dispersive mode (Fig. 4e). Two-dimensional spectra highlight differences between PVA and cornea wave fields (Fig. 4c,f). Clearly, only the A 0 mode is present in porcine cornea (Fig. 4f), as predicted by numerical simulations (Fig. 3).
In porcine cornea, we fit only the A 0 mode to both isotropic and NITI dispersion relations. Isotropic fits produced an estimate of the isotropic shear modulus whereas NITI fits produced estimates of both G and μ (μ = E/3). Using a simplex optimization method (Methods, Supplementary Software), we found dispersion curves that most closely matched the mode structure in 2D spectra. The isotropic model provided a poor fit (Fig. 4f, yellow curve). In contrast, the NITI model closely followed the A 0 mode (Fig. 4f, pink curve). This trend was consistent for IOP ranging from 5-20 mmHg (Fig. 5).
Moduli estimates from the NITI model for all corneas and IOPs are summarized in Fig. 6. We observe a multiple order-of-magnitude difference between estimated Young's modulus (E TI = 3μ TI ) and the shear modulus G TI , consistent with literature values for static tests and OCE measurements. Both moduli increase with IOP over the observed range; however, the orders-of-magnitude difference in the moduli remains consistent across all IOP. www.nature.com/scientificreports/

Discussion
In dynamic OCE, the cornea is typically considered flat, semi-infinite, and isotropic. These approximations may lead to inaccurate interpretations for real corneal geometry and anatomy. For example, the cornea's finite thickness and bounded structure were ignored until recently, when comprehensive simulations and measurements showed that guided modes cannot be ignored 52 .
A flat isotropic layer is also approximate because the cornea is curved. To see the influence on guided mode structure, we performed numerical simulations in flat and spherical bounded layers. The numerical model was similar to the flat layer one ("Methods" section), except the domain was a curved isotropic layer of thickness h = 0.55 mm with outer radius R 0 = 6.5 mm. The axisymmetric solution convolved over a line approximates the AμT source. The wave field and Fourier spectrum along the midline of the propagating wave (measured by AμT-OCE) showed little difference between flat and curved models, suggesting that a flat layer can be a reasonable www.nature.com/scientificreports/  www.nature.com/scientificreports/ approximation for cornea. This result is unsurprising as curvature adds propagation speed variations of order h/R 0 (< 5%) 54,55 . Considering the dramatic effect of corneal anisotropy on measured spectra, curvature can be ignored to first order. Supplementary Note 6 provides a detailed description of this analysis, and the OnScale input file and MATLAB processing functions to reproduce it are available in Supplementary Software. Corneal microstructure and an array of biomechanical studies strongly suggest that the cornea is transversely isotropic rather than purely isotropic. Recently, Brillouin microscopy was used to probe mechanical anisotropy in the cornea based on angle-dependent measurements of the longitudinal wave speed 44,56 . Using a transversely isotropic model, moduli C 11 , C 13 , and C 33 were estimated for porcine and human corneas, with values on the order of GPa 44 . This is a promising technology for measuring the longitudinal stiffness coefficients in transversely isotropic tissues, which can provide complementary information to our NITI model and shear-wave dynamic OCE method. Unfortunately, Brillouin microscopy is insensitive to the shear moduli G = C 44 and µ = C 66 , which dominate the deformational behavior of NITI tissues. Anisotropy in the longitudinal terms is a second-order effect (Supplementary Note 4), and accurate estimates of G and µ will likely be essential in driving both clinical measurements of ocular biomechanics and computational models of corneal deformations.
Here we have shown that a NITI model more accurately characterizes elastic waves in dynamic OCE studies of the cornea. In particular, elastic waves measured in cornea and isotropic PVA phantoms produce markedly different wave fields and spectra (Fig. 4), demonstrating that an isotropic model is not appropriate for cornea.
The NITI model is defined by two shear moduli (G and μ), decoupling tensile/inflation responses from shear responses commonly monitored in torsional tests and dynamic OCE measurements. Based on existing literature, the Young's modulus for cornea is expected to be on the order of MPa, while the shear modulus is on the order of kPa. This is not physically possible for isotropic materials.
The Rayleigh wave speed in a NITI material is almost entirely defined by the modulus G (see Supplementary Note 3). However, the cornea's finite thickness produces guided waves depending on both G and μ, allowing both parameters to be estimated. Theoretically, the optimal way to determine μ is from the phase velocity spectrum of the S 0 mode, which is largely defined by μ (Fig. 2c). Unfortunately, numerical simulations and OCE measurements show that excitation at the air/cornea interface transfers little energy to the S 0 mode, making it nearly impossible to detect. The absence of an S 0 mode is strong evidence of anisotropy, but also means that the A 0 mode alone must be used to evaluate both G and μ.
Using the proposed NITI model, we estimated decoupled Young's (E TI = 3μ TI ) and shear (G TI ) moduli from measurements of elastic wave propagation in porcine cornea. Results agree well with literature values, with www.nature.com/scientificreports/ E TI > 500 kPa and G TI in the range of 6-200 kPa depending on the IOP (Fig. 6), and generally support the observed orders-of-magnitude difference between the two moduli. We note that μ shows greater variance relative to G. For high anisotropy (μ ≫ G), the A 0 mode is increasingly insensitive to μ. Measurement noise inside the elastic wave bandwidth has a greater effect on μ estimates, and uncertainty increases with increasing μ. Thus, while the multiple order-of-magnitude difference between G and μ is accurate, determining the true value of μ requires increasingly higher signal-to-noise ratio as anisotropy increases. This produces large confidence intervals for μ in our measurements. The choice of regularization (used to stabilize the fitting routine, see "Methods" section) can also affect the exact estimate of μ. Consequently, there is far more uncertainty in estimated μ than G values, especially for low SNR measurements. Because of this uncertainty and the need to choose a suitable regularization parameter, future work is needed to refine these estimates. Increasing the OCE measurement SNR should improve sensitivity to μ, as described above. Ideally, new methods designed to excite and accurately measure the A 0 and S 0 modes together could provide more accurate and sensitive estimates of μ.
It is important to note that the NITI model has limitations. As IOP increases, guided modes change noticeably, particularly at very high IOP. Figure 7 shows two-dimensional spectra and best-fit dispersion curves for a porcine eye measured over a larger IOP range (5,15,25,and 35 mmHg). For lower IOP (5 and 15 mmHg), Fourier spectral peaks follow the general A 0 mode shape, and the NITI model closely fits the data (Fig. 7a,b). However, at higher IOP (25 and 35 mmHg), mode shape changes dramatically, and the NITI fit no longer describes it well (Fig. 7c,d).
Two factors must be considered in this high IOP regime-nonlinearity and complex anisotropy. The cornea exhibits nonlinear elasticity that changes at approximately 30 mmHg due to a two-stage deformation process 57 . Deformation is governed by a net matrix response at low strain and collagen stiffness at higher strain. Changes in load-bearing characteristics at high IOP may dramatically change wave propagation and/or induce more complex anisotropic behavior.
Stress-strain testing and polarization-sensitive imaging of collagen alignment suggest that the cornea exhibits a relatively symmetric tensile response at low strain 38 . As strain increases, fiber orientation changes. Transient elastography studies have also observed in-plane anisotropy in the cornea starting at 15-20 mmHg and increasing with IOP 22,28,35,36,39 . This suggests that the cornea becomes more nonlinear and anisotropic at high IOP, and the NITI model no longer adequately describes it. Complex anisotropy models, such as orthotropic or fibril-based models, may be required at high IOP.
Further studies are needed to evaluate the NITI model's clinical utility. As human and porcine cornea differ slightly, it must be tested using OCE measurements on humans. More complex models may also be considered to better estimate corneal stiffness at high IOP. At physiologically relevant IOP, however, we expect the NITI model to perform similarly for human and porcine cornea.
The model's relative simplicity should also facilitate future clinical trials, as it requires a single non-contact measurement, obtained within seconds, to estimate G and μ. While other clinical tools, such as the Ocular Response Analyzer (ORA), have been used clinically to infer both corneal stiffness and IOP 58, sufficient measurement error 59, in addition to patient discomfort 60, hinder widespread use. Further, the relationship between in vivo applanation and corneal Young's Modulus remains unclear 61,62 . The non-contact nature and consistency of AμT-driven OCE strongly support its potential as a practical clinical tool to evaluate corneal elasticity, monitor glaucoma, and study corneal response to ophthalmic interventions.
Non-contact clinical assessment of cornea biomechanics may enable in vivo clinical trials that can provide further insight into the role of biomechanics in cornea function. AμT-OCE can produce accurate maps of mechanical properties, providing reliable, non-contact assessment of corneal biomechanics. Such features make www.nature.com/scientificreports/ it a potentially valuable tool to evaluate cornea in vivo and to develop future procedures to improve vision by modifying tissue microstructure.

Analytic solution.
To derive an analytic solution for guided waves in the cornea, we consider an infinite NITI layer of thickness h and density ρ bounded above by air and below by water. The stiffness tensor contains material constants λ, μ, and G. We assume a plane strain state, consistent with the pseudo-line source generated in our AμT experiments. Here, we briefly overview the derivation for the guided wave solution. Supplementary Note 4 provides a complete derivation. For a NITI material, the elastic wave equations take the following dimensionless form: where u and v are the x-and z-components of the displacement, respectively, and subscripts denote partial differentiation. Assuming harmonic plane wave solutions for displacements in the elastic wave equations for the tissue layer and harmonic acoustic wave solutions for the bounding fluid leads to a secular equation that can be solved for the guided wave frequency-wavenumber dispersion relation (Eq. S4.25 , Supplementary Note 4).
numerical simulation. We developed a finite element numerical model of guided wave propagation in a NITI layer using OnScale (OnScale, Redwood City, CA) 63 . Supplementary Note 5 provides a full description of the model. Briefly, we model the cornea as a thin elastic layer of thickness h = 0.55 mm and density ρ = 1,000 kg/ m 3 , bounded above by air (free surface condition) and below by a layer of water (modeled as an isotropic solid with density 1,000 kg/m 3 , shear wave speed 0 m/s, and a longitudinal wave speed roughly matched to the solid layer). The outer boundaries of the computational domain were set to absorbing conditions. A pressure load was applied to the air-tissue interface with a Gaussian profile in space and super-Gaussian profile in time. The spatial full-width-at-half-max (FWHM) was 600 μm. It was measured with a needle hydrophone (HNC-1000, Onda, Sunnyvale, CA, USA) in air sampling along a 45° line through the AμT transducer focus. The temporal FWHM was 100 μs, also chosen to match AμT experiments. The computational domain was discretized using linear finite elements on a regular rectangular grid with at least 40 elements per elastic wavelength. Simulations were solved using explicit time stepping, and the vertical velocity component was extracted for analysis, similar to OCE experiments where only this component is available. Velocity data were directional and band-pass filtered using the same processing as OCE experiments to remove reverberations from faster wave components.
AμT-OCE to track mechanical waves. To generate elastic waves, we excited samples with acoustic micro-tapping (AμT), a technique using a cylindrically focused, air-coupled ultrasound transducer to induce a localized radiation force at the sample surface 22,64 . The AμT transducer effectively applied a line load to the surface over a wide region relative to the propagation distance of interest, resulting in approximately planar elastic waves (normal to the OCT image plane). The transducer's full-width-at-half-maximum lateral focus was measured as 420 µm. In the phantom experiment, the spatial push width was approximately 600 µm due to its tilt angle relative to the sample 22 . Due to corneal geometry, the axial focus in the porcine experiment was closer to the theoretically measured axial width of 420 µm 53 .
For AμT (as for laser excitation of ultrasound) 65, when an elastic wave is excited by an infinitesimally short (in time) push, the spectral characteristics of the wave are defined by the spatial width of the push. In practice, this pressure confinement can be realized when the push duration is shorter than the shear wave propagation time across the excitation zone. Taking this into account, we utilized a 100 µs pulse duration to generate broadband (up to 4 kHz) mechanical waves and induce tissue displacements on the order of hundreds of nanometers. Because the shear wave excitation is reflection-based, most of the acoustic energy (99.9%) is reflected from the sample surface, leading to acoustic exposures well below the limits for ophthalmic applications 22,66 .
The axial particle vibration velocity of propagating mechanical waves was detected using a phase-sensitive frequency-domain OCT (PhS-OCT) system, which has been described in a previous study 53,[67][68][69] . The M-B mode PhS-OCT system includes a broadband super-luminescent diode with a 1,310 nm center wavelength and 86 nm spectral BW (Denselight Ltd., Singapore), a 90/10 beam splitter, a stationary reference arm, a sample arm integrated with a set of galvo-scanners, and a high-speed spectrometer. The incident power on the sample at focus was approximately 8 mW, providing an OCT SNR of 40 dB. This incident power is well within the safety limits for ophthalmic imaging 70 . The sampling rate of the 1,024-pixel line-scan InGaAs array was 46.5 kHz, determining the A-line rate of the system (temporal resolution). The optical resolution was approximately 15 µm axially and 24 µm laterally.
To track mechanical wave propagation on the sample surface, an external TTL trigger synchronized the PhS-OCT system with wave excitation for each M-scan. All data were collected in an M-B format in which 512 A-scans are repeated in the same location (M-scan) at 256 different horizontal locations (B-scan) across the imaging plane (dx = 54.7 µm), forming a complete M-B scan (1,024 depth × 256 lateral locations × 512 temporal frames) with an effective imaging range of 1.5 mm × 10 mm (axial × lateral). One full M-B scan took 3.66 s. www.nature.com/scientificreports/ The resulting three-dimensional dataset was then used to reconstruct the propagating wave based on the OCT-measured local particle vibration velocity. The axial vibration velocity at a given location ( v z (x, z, t)) was obtained from the optical phase difference �ϕ opt (x, z, t) between two consecutive A-line scans at each location using the following equation 56 : where − was the center wavelength of the broadband light source, − n was the refractive index of the medium, and f s was the sampling frequency. The refractive index in the cornea and phantom were assumed to be 1.38. The system was able to reliably detect displacements greater than ≈ 5 nm.
The surface signal was measured by automatic detection of the sample surface using an edge detection algorithm. Phase data in a 183 μm window below the surface were then extracted and averaged, weighted using one half of a Gaussian window (HWHM = 90 μm, weight decreasing with depth).
fitting experimental data with the niti model. Quantitative moduli estimates in bounded materials require a method to determine the dispersion relation most closely matching observed guided wave modes. Here, we performed this analysis in the frequency-wavenumber domain using a simplex optimization method (fminsearch, MATLAB, MathWorks, Natick, MA). A complete description of the fitting procedure, including rationale for fitting in the frequency-wavenumber domain, is included in Supplementary Note 7.
The theoretical solution presented in Supplementary Note 4 acted as the forward model for optimization. A number of physical parameters were considered fixed, including the corneal density (1,000 kg/m 3 ), corneal longitudinal wave speed (1,540 m/s), and mean corneal thickness (measured from B-mode OCT images). The cornea was bounded from below by water with a density of 1,000 kg/m 3 and longitudinal wave speed of 1,480 m/s. Because we did not observe the S 0 mode in corneal measurements, we extracted only the A 0 mode from the forward model using a mode-tracing routine (similar to Pavlakovic et al. 71 ).
A two-dimensional Fourier transform was applied to OCE-measured surface velocity data to generate a normalized power spectrum. An optimization routine based on the simplex method estimated both shear moduli, G and μ, by fitting the experimentally obtained 2D spectra with the analytic solution (Eq. S4.25 , Supplementary Note 4). For each sample, the thickness (h) was measured using automated edge detection algorithms applied to the OCT structural image collected at the start of each MB scan sequence. At each iteration, a dispersion relation for the A 0 mode was computed for the current iterate (G i , μ i ) based on the forward model. The average power within a small 7-point Gaussian window centered on this dispersion curve was computed, and the algorithm updated iterates of G and μ to maximize this quantity. A regularization term is included to ensure that the ratio μ/λ remains small, consistent with the nearly-incompressible assumption. To compare isotropic and NITI model fits, a goodness-of-fit (GOF) metric was defined based on the maximum energy at each given frequency. When this goodness-of-fit metric approaches a value of 1, the model's dispersion curve accurately captures all of the energy of the A 0 mode (see Supplementary Note 7 for details).
isotropic phantom preparation. A homogenous, isotropic, elastic phantom with controllable mechanical properties was created to experimentally measure wave behavior in a thin plate model. It was fabricated using a similar protocol to that described by Kharine et al. 72 Briefly, polyvinyl alcohol (PVA) (146-186 kDa, > 99% hydrolyzed, CAS: 9002-89-5, Sigma-Aldrich Corp., St. Louis, MO, USA) was added to a 4:1 mixture of dimethylsulfoxide (DMSO, CAS: 67-68-5, EMD Millipore Corp) and water at a concentration of 4 wt%. To tune the phantom's optical properties, we added 0.025 wt% titanium dioxide nanoparticles. The solution was covered and stirred at a temperature of 95 °C for approximately 1 h until the PVA was completely dissolved. The solution was degassed in a vacuum chamber to remove any air bubbles before casting in a round mold with a radius of 10 cm. Phantom thickness (0.6 mm, measured by OCT) was controlled by the amount of PVA solution poured into the mold and allowed to settle. The mold was stored at − 20 °C for at least 12 h, or until the phantom was completely frozen. The phantom was then thawed at room temperature, completing one freeze-thaw cycle. After casting, phantoms were removed from their molds and placed in a water bath for at least 48 h to allow the DMSO to diffuse out. Prior to imaging, the PVA phantom was suspended on top of water to force asymmetric boundary conditions similar to those of the cornea 52 . porcine cornea samples. Porcine eyes were enucleated immediately after death and stored in physiological saline until imaging. All OCE measurements were performed within 1 h of euthanization. The whole porcine eyeball was placed into a custom-built holder with a hemispherical cup filled with saline-moisturized cotton to provide an in situ environment. The eye globe was oriented cornea side up with the optic axis vertical and aligned with the OCE scanning beam. A 23-gauge needle connected to an infusion reservoir was inserted through the sclera to control intraocular pressure (IOP). The reservoir height was adjusted to maintain IOP between 5 and 35 mmHg. The mean thickness of the cornea samples was 0.71 ± 0.11 mm.
All studies were carried out in accordance with institutional guidelines and regulations for tissue studies. All experimental protocols followed standard operating procedures established by the University of Washington for the use of animal tissue acquired from an abattoir in research studies.

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 on reasonable request. In addition, we have www.nature.com/scientificreports/ included a Supplementary Software Library containing the MATLAB scripts and functions used in this study, as well as the OnScale finite element input files. A detailed description of the functions and scripts is provided in the Supplementary Software Documentation. We also include three example MATLAB data files: (1) Example OCE data from one porcine cornea measurement, (2) example OnScale results for the NITI guided wave model, and (3) example OnScale results for the spherical layer model.