Nonlinear fourth-order elastic characterization of the cornea using torsional wave elastography

Measuring the mechanical nonlinear properties of the cornea remains challenging due to the lack of consensus in the methodology and in the models that effectively predict its behaviour. This study proposed developing a procedure to reconstruct nonlinear fourth-order elastic properties of the cornea based on a mathematical model derived from the theory of Hamilton et al. and using the torsional wave elastography (TWE) technique. In order to validate its diagnostic capability of simulated pathological conditions, two different groups were studied, non-treated cornea samples (n=7), and ammonium hydroxide (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$NH_4OH$$\end{document}NH4OH) treated samples (n=7). All the samples were measured in-plane by a torsional wave device by increasing IOP from 5 to 25 mmHg with 5 mmHg steps. The results show a nonlinear variation of the shear wave speed with the IOP, with higher values for higher IOPs. Moreover, the shear wave speed values of the control group were higher than those of the treated group. The study also revealed significant differences between the control and treated groups for the Lamé parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μ (25.9–6.52 kPa), third-order elastic constant A (215.09–44.85 kPa), and fourth-order elastic constant D (523.5–129.63 kPa), with p-values of 0.010, 0.024, and 0.032, respectively. These findings demonstrate that the proposed procedure can distinguish between healthy and damaged corneas, making it a promising technique for detecting diseases associated with IOP alteration, such as corneal burns, glaucoma, or ocular hypertension.


Introduction
Maintenance of external ocular morphology and function is related to the biomechanical properties of the cornea [1].Several corneal disorders, such as keratoconus, ectatic disorders, and chemical eye burns, emerge as the most important due to their severity, impact on patient quality of life, and incidence [2][3][4].Alkali burns are the most common emergency related to ocular trauma [5][6][7] (84% of cases [4]), resulting in high visual morbidity [8].Significant corneal damage was determined when the pH was greater than 11.5 [9].Among the most common alkali chemical agents, ammonium hydroxide ( NH 4 OH ) has the fastest penetration rate (<3 min) [10].As a consequence, chemical injuries of the eye produce extensive damage to the ocular surface epithelium, cornea, anterior segment, and limbal stem cells resulting in permanent unilateral or bilateral visual impairment [11].Damage to the corneal and conjunctival epithelium may lead to opacification and neo-vascularization of the cornea [11], causing the appearance of an acute increase in IOP due to shrinkage and contraction.
To provide clinically relevant critical information for the diagnosis of corneal diseases, elastography is proposed as an emerging non-invasive imaging method.Very recently, Optical Coherence Elastography (OCE) was used to detect the propagation of induced elastic waves to biomechanically assess the cornea [12][13][14][15][16], lens [17] and retina [18].The advantages of this and other recent elastography techniques include the microscale sensitivity in motion detection [19][20][21] and the noncontact approach, along with the image's microscale resolution.However, the direction of vibration of the particles, nearly perpendicular to the median plane of the cornea, results in the generation of guided waves motivated by the relationship between wavelength and tissue thickness, which prompts the use of complex guided wave models [22].As additional disadvantages are the low frame rate in 2D imaging, long acquisition imaging times, and exact positioning where slight motion could cause image artifacts [23,24].
Beyond the standard of elasticity maps, a more precise, pressure-and operator-independent interpretation of the measurements might be obtained considering nonlinearity since the dependence of the stiffness modulus with deformation is correlated with IOP pressure.In this sense, nonlinear elastic constants may be much more sensitive to specific diseases, facilitating an early diagnosis [25].In addition, most biological tissues have a nonlinear stress-strain behavior under variable amounts of pressure [26].This nonlinear elasticity can cause biological tissues to become more stressed under increased pressure [27].
The formula proposed by Hamilton et al. [28], a strain energy equation that neglects compression terms, has been applied to the mechanical characterization of soft tissues both with the acoustoelasticity technique [29] and in uniaxial tensile tests [27,30] to determine nonlinear elastic parameters.In acoustoelasticity, Hamilton's formula was used for the breast [26], liver [31,32], or porcine kidney [33] considering only the shear modulus ( ) and the third-order elastic constant (A), and neglecting the fourth-order elastic constant (D) because it was hypothesized that the shear wave displacements were smaller than static compression.According to cornea applications, the influence of the pressure applied in the tissue in determining its stiffness has been studied.In this regard, ex-vivo studies have demonstrated that the stiffness of the cornea and sclera increases with increasing the intraocular pressure [34][35][36][37][38]. Regarding the nonlinear characterization of the cornea, to our knowledge, few studies have been performed.An ex-vivo study has been carried out to characterize nonlinear parameters of the cornea using uniaxial destructive tests [27].However, to date, no work has studied Hamilton-Zabolotskaya's fourth-order parameters in the cornea employing non-destructive tests, with potential application in clinical practice.
Torsional waves are postulated to be a crucial tool, sensitive to the measurements of nonlinear parameters [25,39].Based on the studies presented, there is a knowledge gap about the effectiveness of the different nonlinear analysis methods that can be applied to the existing TWE technique.This work proposes a mathematical model based on Hamilton et al. theory [28] for the porcine cornea nonlinear behavior characterization of two groups, non-treated cornea samples and samples treated with an ammonium hydroxide concentration ( NH 4 OH ), using the TWE tech- nique.The geometry of the wave emitting part and the receiving part are adapted to the curvature of the cornea, an evolution from the original design used in previous studies for the characterization of cervical tissue [40,41].An important focus is the vibration of the particles due to torsion waves.In this case, said vibration is contained in the median plane of the cornea (in-plane propagation).A distinctive feature of this technique was that we obtained information about in-plane shear deformation in the cornea, therefore, the influence of guided waves, likely present in out-of-plane propagation disappears [34], along with the associated complexity, unlike the OCE technique (out-of-plane propagation).
The remainder of this paper is organized as follows.Section 2 describes the preparation of samples, the Torsional Wave Elastography technique, the proposed mathematical model, and the data post-processing.Sections 3 and 4 describe the results and the discussion, respectively.Finally, Section 5 offers the conclusions of this research and suggestions for the future.

Tissue preparation
Due to the similarity of the mechanical behavior between porcine samples and the human ones [42], in this study, porcine eye globes were used to characterize the nonlinear behavior of the cornea.Fourteen porcine eye globes were enucleated postmortem from a local abattoir and placed immediately in phosphate-buffered saline solution (PBS, pH 7.4) to maintain tissue hydration and to prevent the cornea from hardening until Torsional Wave Elastography measurements were performed.The solution was prepared using di-Sodium Hydrogen Phosphate anhydrous (Reag.Ph.Eur.99%), Potassium di-Hydrogen Phosphate (Reag.Ph.Eur.99% purity), and Sodium Chloride (USP, BP, Ph.Eur.JP 99%) from Panreac AppliChem.All samples were tested at room temperature within eight hours postmortem.Exposures to household cleaning products commonly involve the eyes.This exposure occurred in 8.4% according to a study carried out in United States [43].Ammonium hydroxide ( NH 4 OH ) solution (EMSURE ACS, Reag.Ph Eur 28-30%) is one of the typical cleaning products.For that reason, two different groups were studied, nontreated cornea samples (n=7) and NH 4 OH-treated samples (n=7).Ammonium hydroxide was mixed with distilled water in a similar concentration (3 mM NH 4 OH at 10% v/v) to that of cleaning products [44].The exposure time for the treated samples was 5 min, while no treatment was applied to the control group.After that, all the treated samples were washed in PBS before examination with TWE.

Torsional wave elastography
The same torsional wave device was used as in the study carried out by Torres et al. [34], where the sensor (receiving ring -4 mm base) and excitation (emitting diskexternal and internal diameters were 13 mm and 9.6 mm respectively) components were assembled.Torsional waves were generated by rotating the emitter in contact with the cornea, which are a type of shear waves that propagate axisymmetrically until they reach the receiver.The displacement peaked at around 2 m.The geometry of the sensor was optimized so that the volumetric displacement components were negligible compared to the shear displacement components, and this is achieved using a receiving ring whose center coincides with the center of the emitting disk [45,46].That configuration makes it possible to obtain, after the analysis of the received signal, the shear mechanical parameters in-plane that interrogate the properties of the cornea.In addition, the receiving ring had an internal curvature that was adapted to the geometry of the cornea [47].
The diagram of the experimental configuration is shown in Fig. 1.A multichannel AD/DA converter with a 192 kHz sampling rate was employed to emit the signal and receive it after the propagation through the cornea.All the samples were measured with a single sinusoidal pulse of 1000 Hz, according to previous references in the literature [40,48].An amplifier was used to output a 25V peak-to-peak signal.The received signal by the TWE probe was preamplified (40 dB gain) before reaching the digital to analog converter.In order to reduce random and high-frequency noise, an average of 16 signals and a 5 kHz low-pass filter were applied.The total measurement time was 3.2 s, composed of 16 intervals of 200 ms.A calibration measurement was then subtracted from the received cornea signal to counterbalance crosstalk effects.All elements were computer-controlled using high-speed communication ports and a Matlab environment (R2018b, The MathWorks Inc., Natick, MA, USA).Based on previous results, no guided waves were produced inplane shear deformation caused by TWE measurement [34].Therefore, the phase speed was calculated by subtracting a quarter of the period (inverse of the received signal central frequency) from the first signal peak.Due to the small propagation distance between the emitter and the receiver (2.8 mm), the curvature of the cornea was approximated as a straight line [49].
Each eyeball was placed in a custom-made holder and the axis of the emitter was aligned with the optical axis of the eyeball.A needle was inserted into the anterior chamber of the eye, which was connected to a saline solution reservoir.The IOP was modulated by adjusting the height of the reservoir [50].The tests increased IOP from 5 to 25 mmHg with steps of 5 mmHg.The IOP was held constant at each IOP level for 2 min before TWE measurements.All corneas were measured three times by repositioning the probe.The pressure applied with the TWE probe was small enough so that the received wave had the higher possible amplitude to avoid slippery conditions [41] and did not influence the IOP controlled by the saline solution reservoir.

Mathematical model
Considering the deformation due to the intraocular pressure and that due to shear wave propagation (see Fig. 2), the deformation gradient tensor is given by: where is the stretch (defined as = 1 + , is the strain) in the two circumferential directions and a strain of −2 in the radial direction, which ensures incompressibility ( J = det(F) = 1 ).In addition, torsional waves generate shear strain in the cornea, which is given by (x 12 , t) ( x 12 shear strain direction -Fig.2, t time), considering that the strain due to the torsional wave is much smaller than the strain induced by changing the IOP (i.e.,  << ).
The Green-Lagrange strain tensor is given by: The strain energy density equation proposed by Hamilton et al. [28] is, Where I 2 and I 3 are the invariants of the strain tensor, is the Lamé constant, A is the third-order elastic constant, and D is the fourth-order elastic constant.These terms are sufficient for describing shear deformation when the energy stored in compression is comparatively insignificant.The invariants are calculated as, ( 1) The second Piola-Kirchhoff (PK) stress tensor may be computed as S = W∕ E .Using this equation, equation 3, and taking into account that  <<  [34, 51] we obtain (see Appendix (Sect.6) for more details), Considering the porcine cornea as a spherical shell of radius ( R = 8.45mm ) [47], and a typical thickness of a porcine cor- nea ( = 1mm ), the equilibrium equation is, where IOP is the intraocular pressure in mm of Hg, and is the Cauchy stress tensor.
The relation between the Cauchy stress tensor and the second Piola-Kirchhoff stress tensor is, ( 4) Fig. 2 Scheme of the stress and strain field.A differential element of the cornea is represented in blue.Subfigures a and b show the normal stresses to which the cornea is subjected in three directions according to the reference system adopted ( e 1 , e 2 , e 3 ).The strain field is representated in Subfigure c.Finally, subfigure d shows the in-plane strain generated by the torsion wave To obtain a relationship between the torsional wave speed ( c s ) and the stretch, the Green-Lagrange strain tensor has been decomposed into static ( E s ) and dynamic ( E d ) strain as follows, The dynamic strain has been obtained considering the displacements associated with an in-plane ( e 1 e 2 ) torsional wavefront.
Taking into account the following shear wave equation, We arrive at a nonlinear relationship between the shear wave speed ( c s ) and the parameters , A,and D: For a given value of IOP and the , A, and D parameters, we can obtain using equation 10.Afterwards, the corresponding torsional wave speed ( c s ) and the value of are intro- duced to minimize the equation 14 and obtain the parameters (10) (11) , A and D. It is iterated until the error of the three aforementioned parameters is minimized for each sample of each group with genetic algorithms (Matlab -R2018b, The Math-Works Inc., Natick, MA, USA).

Results
Figure 3 shows the shear wave speed of the TWE technique for both groups, control and NH 4 OH , and for each IOP (5-25 mmHg with 5 mmHg steps).We observed that the values of the control group were higher than those of the treated group.Furthermore, the shear wave speed varied nonlinearly with the IOP, with an increase in the speed observed as the IOP increased.
The parameters reconstructed from the proposed model (Eqs.10, and 14) for the control and NH 4 OH groups are displayed in Fig. 4. The boxplots of the Lamé parameter , third-order elastic constant A, and fourth-order elastic constant D showed higher values in the control group.The Mann-Whitney-Wilcoxon test was used to study the significance level of the reconstructed parameters for both groups.A value of p < 0.05 was considered statistically significant for the difference testing.We observed that the difference between the Lamé parameter , third-order elastic constant A, and fourth-order elastic constant D for the two groups was significant, p = 0.010 , p = 0.024 , and p = 0.032 respectively.Fig. 5 displays the representative fit (comparing shear wave speed to IOP) of the experimental data for both a control sample and an NH 4 OH-treated sample using our proposed mathematical model.The reconstructed parameters for the control fit were = 25.54 kPa, A = 767 kPa, D = 303 kPa, while for the NH 4 OH fit were = 7.33 kPa, A = 238 kPa, D = 152 kPa.Across all samples from both groups, the coefficient of determination ( R 2 ) demonstrated a strong fit with a minimum value of 0.9.
The relationship between stretch and IOP was represented for control and NH 4 OH cornea samples using the Eq. 10 for a given , A, and D parameters (see Fig. 6).It is clear that the IOP varies nonlinearly with the stretch and that the slope of the curve, which is determined by the nonlinear parameters A and D, is higher in the control group.

Discussion
Several studies have suggested that nonlinearity could enhance the current diagnosis in elastography since the operator dependency is decoupled by correlating the stiffness and the tissue deformation [25].In this study, a mathematical model based on Hamilton et al. theory [28] is proposed for the porcine cornea nonlinear behavior characterization of two groups, seven non-treated cornea samples and seven samples treated with NH 4 OH , using the TWE technique.
Recent work has focused on corneal mechanical characterization.Weng et al. [52] used a high-frequency ultrasound elastography system based on an external vibrator and ultrafast ultrasound imaging to estimate the viscoelasticity of porcine cornea ex-vivo by using a Lamb wave model.A study carried out by Torres et al. [34] developed an elastography technique based on torsional waves (TWE) adapted to the specificities of the cornea to characterize the viscoelasticity of two groups, a control group and one group of alkali burn treatment ( NH 4 OH ).The TWE technique reflected mechanical properties changes after treatment, showing a high potential for clinical diagnosis due to its rapid performance time.
According to the nonlinear characterization of the cornea, few studies have been performed.By using a similar experimental setup than [52], Zhang J. et al. [53] assessed the nonlinear elastic properties of the cornea and ciliary body using the hyperelastic Blatz model.However, they did not compare the nonlinear response of the healthy with other model corneas.Ashofteh Y. et al. [27] evaluated which hyperelastic model (Hamilton-Zabolotskaya model, Ogden model, and Mooney-Rivlin model) could best describe the nonlinear behavior of the cornea to discriminate structural changes in a damaged cornea (NaOH treatment) with quasi-static uniaxial tests.However, they did not study the nonlinearity of the damaged cornea with one of the most common alkali chemical agents and the fastest penetration rate (ammonium hydroxide -NH 4 OH ).Another disadvantage is that the uniaxial tensile test is a destructive procedure, unlike our study.Thus, the cornea samples were characterized by high deformation ratios.
The correlation between the shear wave speed and the (IOP) was studied for both groups, the control, and NH 4 OH group, and for each IOP 5-25 mmHg with 5 mmHg steps, see Fig. 3. Shear wave speed values for the control group are higher than those of the treated group.This was also evidenced in the work carried out by Torres et al. [34], which is motivated by significant changes in the composition of the extracellular matrix of the cornea after a corneal chemical injury [54].Collagen fibrils behave nonlinearly for high intraocular pressure values [55].Results from Fig. 3 show how the shear wave speed depends on the IOP nonlinearly, which is consistent with previous literature studies [34,55].
A complex microstructure governs the nonlinear behavior of the cornea.It is composed of a hydrated proteoglycan matrix with collagen fibrils embedded in it [56].The collagen fibrils are packed and arranged in the circumferential direction [57,58].The TWE technique used in this work interrogates the properties of the cornea in the plane e1-e2 (in-plane -deformations produced by the torsion wave in that plane, see Fig. 2), which correspond with the collagen fibrils of the cornea, unlike the OCE and ultrasound techniques which measured in the plane e1-e3 or e2-e3 (out-ofplane), the matrix of the cornea.Given this, our technique's advantage lies in the wave propagation, not generating Lamb waves, which was demonstrated in the work carried out by Torres et al. [34].However, Lamb waves are generated by measuring out of the plane, and a more complex study is required [52,59].Figure 3 shows how the reconstructed values of shear wave speed for the untreated cornea are higher than previous studies present in the literature [52,53,53] for different values of IOP, which may be motivated by in-plane (plane e1-e2) measurements.
In recent years, there has been a growing interest in the mechanical characterization of the cornea.Researchers, such as Weng et al., have employed a high-frequency ultrasound elastography system to estimate the viscoelasticity of the porcine cornea ex-vivo.Similarly, Torres et al. developed an elastography technique based on torsional waves that was adapted to the specificities of the cornea.The TWE technique demonstrated its potential for clinical diagnosis due to its rapid performance time and ability to reflect mechanical property changes after treatment.While some studies have focused on the nonlinear characterization of the cornea, such as Zhang J. et al.'s assessment of the nonlinear elastic properties of the cornea and ciliary body using the hyperelastic Blatz model, few have compared the nonlinear response of healthy corneas to other model corneas.Ashofteh Y. et al. evaluated which hyperelastic model could best describe the nonlinear behavior of the cornea to discriminate structural changes in a damaged cornea, yet they did not study the nonlinearity of the damaged cornea with one of the most common alkali chemical agents and the fastest penetration rate.Furthermore, the uniaxial tensile test is a destructive procedure, which sets our study apart as the cornea samples were characterized by high deformation ratios.Our study proposes a nonlinear fourthorder elastic model that can differentiate between healthy and damaged corneas.We found significant differences in elastic constants between the control and NH4OH groups.The TWE technique has the potential to measure nonlinear parameters and detect diseases related to intraocular pressure, making it a promising concept with diagnostic potential for diseases like glaucoma.The correlation between the shear wave speed and the intraocular pressure was studied for both groups, the control, and NH4OH group.Shear wave speed values for the control group are higher than those of the treated group.This was also evidenced in the work carried out by Torres et al., who were motivated by significant changes in the composition of the extracellular matrix of the cornea after a corneal chemical injury.Results show how the shear wave speed depends on the intraocular pressure nonlinearly, which is consistent with previous literature studies.The nonlinear behavior of the cornea is governed by a complex microstructure.In this study, the TWE technique interrogated the mechanical properties of the cornea in the plane e1-e2, which corresponded with the collagen fibrils of the cornea.This is in contrast to other techniques such as OCE and ultrasound, which measured in the matrix of the cornea.Our technique's advantage lies in the wave propagation, which reconstructs values of shear wave speed for the untreated cornea that are higher than previous studies present in the literature for different values of intraocular pressure.This may be motivated by in-plane measurements.
The aforementioned nonlinear fourth-order elastic model is based on the strain energy density equation by Hamilton et al. [28], which allows a nonlinear fourth-order elastic characterization of the cornea using Torsional Wave Elastography.The model is governed by the Lamé constant ( ), the third-order elastic constant (A), and the fourth-order elastic constant (D).These terms are sufficient for describing shear deformation when the energy stored in compression is comparatively insignificant.This is achieved because torsional waves present the advantage of isolating a pure shear motion, thus, eliminating the generation of spurious compressional waves due to the configuration emitterreceiver of the probe [45,60].Genetic algorithms (Matlab -R2018b, The MathWorks Inc., Natick, MA, USA) were employed to minimize the discrepancy between the pairwise measurements of torsional wave speed ( c s ) and IOP, equa- tions described in detail in the Mathematical model section (Eqs.10, and 14). Figure 4 shows the boxplots of the values the Lamé parameter , third-order elastic constant A, and fourth-order elastic constant D for both groups, control, and NH 4 OH group.
The reconstructed parameters for the control group fit were = 25.9 ± 6.85 kPa, A = 215.09± 133.58 kPa, D = 523.50± 161.35 kPa, while for the NH 4 OH group fit = 6.52 ± 3.94 kPa, A = 44.85± 13.30 kPa, D = 129.63 ± 76.09 kPa.A minimum value of 0.9 for the coefficient of determination ( R 2 ) was obtained for both groups' samples.Generally, the reconstructed parameters for the control group are notably higher than those for the NH 4 OH group.This discrepancy is particularly evident in the parameters A and D when observed on the Stretch-IOP curve (Fig. 6).The curve shows a steeper slope as the strain in the control corneas increases, providing superior contrast when calculating A due to the smaller data deviation compared to the fourth-order elastic parameter D. Alireza et al. [27] employed the same nonlinear model, using measurements from uniaxial tensile tests for a control group and a NaOH-treated group of porcine corneas.However, the parameters , A, and D values that we obtained for the control group differ from those in their study.This variation could be due to the differences in the stretch analyzed; the strain in their study is higher, and it is known that the shear modulus ( ) increases with the applied strain [55].However, measurements by Bernal et al. [26] of the nonlinear parameter A in tissues align with the values obtained in our study, as the order of magnitude of the shear deformation is similar.To compare the two groups in our study, we applied the Mann-Whitney-Wilcoxon test to assess the significance level of the reconstructed parameters.The differences between the Lamé parameter , the third-order elastic constant A, and the fourth-order elastic constant D for the two groups were significant (p = 0.010, p = 0.024, and p = 0.032, respectively).The threshold for statistical significance was set at p < 0.05 .Despite studies in the literature often disregarding the fourthorder parameter D for small deformations [26,29], our results highlight the potential of parameter D as a valuable biomarker for distinguishing between healthy and alkali-burned porcine eyes within the studied deformation range.Figure 6 demonstrates the relationship between the stretch and IOP for both groups, using equation 10 after reconstructing the three nonlinear parameters.This figure depicts the nonlinear relationship between IOP and stretch and shows that the slope of the curve, determined by the nonlinear parametersA and D, is steeper in the control group.
In Fig. 6, the relationship between the stretch and the IOP was represented for both groups, using equation 10 after reconstructing the three nonlinear parameters.We can observe the nonlinear behavior between IOP and the stretch, and that the slope of the curve, which is determined by the nonlinear parameters A and D, is higher in the control group.
The proposed method presents some limitations that need to be exposed.The TWE technique only receives a signal; therefore, no 2D image is obtained that could provide us with additional information, such as focal changes in tissue consistency.Besides, it is difficult to change the IOP in conventional examination procedures.However, according to daily IOP fluctuations of the human corneas [61], this technique could be an emerging concept with the diagnostic potential to measure the nonlinear parameters and to detect diseases associated with intraocular pressure such as corneal burns, glaucoma, or ocular hypertension.Ocular inspection requires that the patient does not make involuntary movements that cancel the measurement, which is why the proposed technique has the advantage of measuring in a short period of time.The tonometry technique, currently used in the diagnosis of glaucoma, consists of flattening the cornea to calculate IOP, which can be uncomfortable for the patient [62].However, the TWE technique only requires a low contact pressure between the probe and cornea to perform the measurement.In addition, such a technique could diagnose diseases related to IOP by correlating it with shear modulus.In future work, the study of the correlation between IOP and shear stiffness for diseases such as glaucoma is proposed using the nonlinear method.

Conclusions
In our study, we propose a mathematical model utilizing the torsional wave elastography (TWE) technique to characterize the cornea's nonlinear fourth-order elastic properties based on the theory developed by Hamilton et al. [28].We conducted our analysis on two groups of cornea samples, one untreated (n=7) and the other treated with ammonium hydroxide ( NH 4 OH ) (n=7).Using a torsional wave device to measure in-plane, we incrementally increased the intraocular pressure (IOP) from 5 to 25 mmHg in 5 mmHg increments, and observed that the untreated group had higher shear wave speed values when compared to the treated group.Furthermore, the shear wave speed exhibited nonlinear variation with IOP, with an increase observed as intraocular pressure increased.Our analysis of reconstructed parameters revealed significant differences between the untreated and treated group for the Lamé parameter , third-order elastic constant A, and fourth-order elastic constant D. Additionally, we demonstrated that the stretch exhibited nonlinear variation with IOP.Through the correlation of stiffness and tissue deformation, this technique has the potential to detect diseases associated with IOP alteration, such as corneal burns, glaucoma, or ocular hypertension, by decoupling operator dependency and measuring nonlinear parameters.

Appendix: Mathematical model
The second Piola-Kirchhoff (PK) stress tensor may be computed as S = W∕ E.
The strain energy density equation proposed by Hamilton et al. [28] is, where I 2 and I 3 are the invariants the strain tensor, is the Lamé constant, A the third-order elastic constant, and D the fourth-order elastic constant.The invariants are calculated as, Therefore, Assuming incompressibility, the deformation gradient tensor is given by: where is the stretch (defined as = 1 + , is the strain) in the two circumferential directions and a strain of −2 in the radial direction, which ensures incompressibility ( J = det(F) = 1 ).In addition, torsional waves generate shear strain in the cornea, which is given by (x 12 , t) ( x 12 shear strain direction -Fig.2, t time), considering that the strain due to the torsional wave is much smaller than the strain induced by changing the IOP (i.e.,  << ).
The Green-Lagrange strain tensor is given by: (15) The invariants and their respective derivatives to obtain S 11 are: Therefore, Considering the porcine cornea as a spherical shell of radius ( R = 8.45mm ) [47], and a typical thickness of a porcine cor- nea ( = 1mm ), the equilibrium equation is, where IOP is the intraocular pressure in mm of Hg, and is the Cauchy stress tensor.The relation between the Cauchy stress tensor and the second Piola-Kirchhoff stress tensor is, We obtain a relationship between the intraocular pressure of the cornea and the parameters of the nonlinear model , A, and D.
To obtain a relationship between the torsional wave speed ( c s ) and the stretch, the Green-Lagrange strain tensor has been decomposed into static ( E s ) and dynamic ( E d ) strain as follows, The dynamic strain tensor has been obtained considering the displacements associated with an in-plane ( 1 e 2 ) torsional wavefront.
Taking into account the following wave equation, Next, we calculate the S 12,1 term of the second Piola-Kirch- hoff stress tensor.Now, the invariants and their respective derivatives to obtain S 12,1 : (31) The relationship between the derivative of the dynamic strain tensor terms and the displacements associated with the shear was obtained with the equation 32.
Deriving and simplifying the expression for S Funding Funding for open access publishing: Universidad de Granada / CBUA.

Fig. 1
Fig. 1 Experimental setup for the nonlinear mechanical characterization of the cornea using TWE

Fig. 3
Fig.3 Comparison of the shear wave speed (in m/s) from TWE (Torsional Wave Elastography) for both groups (control -black solid line boxplots and NH 4 OH -red solid line boxplots), and intraocular pressure (IOP in mmHg)