Addressing Parameter Variability in Corneal Biomechanical Models: A Stepwise Approach for Parameters’ Optimization

TL;DR

This paper introduces a stepwise method to optimize corneal material model parameters, aiming to reduce variability and improve reliability in biomechanical modeling.

Contribution

A novel stepwise approach is proposed to separately estimate isotropic and anisotropic components in corneal biomechanical models.

Findings

The stepwise approach reduces variability in parameter sets by isolating isotropic and anisotropic components.

The methodology improves the reliability of corneal material model parameters for clinical and biomimetic applications.

Abstract

Biomechanical modeling of the cornea is crucial for understanding the progression of some ocular diseases and optimizing surgical treatments. However, hyperelastic non-linear material models, such as those used for corneal tissue, often yield highly variable parameter sets in the scientific literature, influenced by factors like the chosen optimization intervals and differences between tensile and inflation test curve optimization, both of which are addressed in this study. This variability complicates the understanding of corneal mechanical properties. In this research, the aim is to optimize and calibrate the key parameters of the corneal material model, particularly focusing on c1, c2, k1 and k2, using the Holzapfel–Gasser–Ogden (HGO) hyperelastic model, and a novel methodology is proposed that separately estimates the isotropic and anisotropic components in a stepwise manner, addressing the issue of multiple parameter sets fitting experimental curves similarly. This approach helps to standardize corneal material models and improve the reliability of parameter estimations. Moreover, accurate biomechanical characterization within this framework contributes not only to clinical applications but also to biomimetics, inspiring the design of artificial corneal substitutes and bioengineered materials.

Figures

Full text

1. Introduction

The human eye receives over 90% of external sensory information, highlighting the critical importance of ocular health for quality of life [1]. The cornea, contributing approximately two-thirds of the eye’s refractive power [2], plays a vital role in this process. This avascular structure, with an average thickness of 550 μm centrally and 750 μm peripherally, is primarily composed of stroma, which accounts for 90% of its thickness and consists of collagen fibers embedded in a proteoglycan matrix and arranged in a region-specific orientation [3].

The cornea’s mechanical behavior is commonly modeled using anisotropic hyperelastic models [4,5], essential to optimize surgical outcomes and understanding diseases like ectasia, characterized by structural weakening of collagen fibers [6]. Hyperelastic models, such as the Gasser–Holzapfel–Ogden (GHO) and Holzapfel–Gasser–Ogden (HGO) models, have become standard in corneal studies due to their capacity to describe anisotropic properties and fiber dispersion [7,8,9,10]. These models assume near-incompressibility and provide parameters representing isotropic stiffness, fiber stiffness, and non-linear deformation. However, parameter variability across studies, influenced by differences in experimental techniques, tissue properties, and mathematical simplifications, complicates direct comparisons and the reliability of conclusions between studies [11,12,13,14].

Tensile and inflation tests are the primary experimental methods to estimate corneal parameters. However, the two approaches yield inconsistent results due to differences in testing conditions and the mathematical models applied. For example, studies using tensile tests, such as those by [11,12,13], have reported wide variability in parameter values, while inflation tests, as demonstrated by [14,15], highlight the sensitivity of material parameters to experimental conditions.

Experimental comparisons have shown that strip extensometry and inflation produce markedly different responses [16] and that multiple loading scenarios are critical for accurate parameter identification [17]. Moreover, recent advances such as Brillouin microscopy have enabled non-invasive biomechanical characterization [18], while age- and IOP-related collagen recruitment studies further emphasize the variability in corneal mechanical response [19]. Recent ex vivo and animal studies also demonstrate discrepancies between uniaxial and biaxial testing outcomes [20], highlighting the difficulty of establishing universally valid parameter sets.

Variations in geometric measurements, non-uniformity of corneal properties, and model-specific assumptions contribute to these discrepancies. Furthermore, high correlations between experimental data and model predictions (high correlation coefficients $[eqn]$ ) do not always ensure unique or physically valid parameter sets. This issue arises when parameter compensations occur, leading to multiple sets that fit experimental data similarly. Consequently, it becomes necessary to establish robust methods for parameter estimation that minimize variability and improve reliability.

This study aims to address these challenges by (1) developing a method to characterize corneal material parameters of the HGO model through a novel stepwise optimization approach that separates the estimation of the isotropic ground-matrix components from the anisotropic fiber parameters to reduce variability and ensure physical consistency in the estimated parameters; (2) proposing parameter bounds and optimization algorithms to demonstrate the high variability inherent in the characterization problem, driven by the large number of parameter combinations that yield similarly high correlation with the experimental curve.

Beyond its clinical relevance, this research also fits within a biomimetic perspective. The cornea, combining transparency, anisotropy, and mechanical resilience, serves as an inspiring model for bioengineered materials. Advances in contact lenses, tissue substitutes, and artificial corneal constructs seek to replicate these properties, highlighting how biomechanical modeling can guide the design of biomimetic solutions that emulate their structural and functional efficiency. In this context, an accurate characterization of the cornea’s hyperelastic properties is essential to ensure that such designs achieve physiological-like behavior.

2. Materials and Methodology

In this research, a detailed methodology is implemented to optimize the material parameters of corneal tissue ( $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ ) using tensile and inflation tests within the Holzapfel–Gasser–Ogden (HGO) model framework [7]. This model is adopted because it combines physiological realism with simplicity: it separates the isotropic ground matrix ( $[eqn]$ , $[eqn]$ ) from collagen–fiber reinforcement and reproduces the cornea’s characteristic exponential stiffening ( $[eqn]$ , $[eqn]$ ). Other models, such as the five-parameter Mooney–Rivlin formulation [21], can also capture pronounced nonlinearity; however, they make it more difficult to distinguish which parameters govern the initial, matrix-dominated segment of the curve and which govern the later, fiber-driven non-linear response.

Both tensile and inflation tests are mathematically defined to simulate the cornea’s mechanical behavior under physiological conditions. The methodology is summarized in the scheme shown in Figure 1.

Python version 3.11.9 (Python Software Foundation) is used for data analysis, parameter optimization, and curve fitting. Core libraries were NumPy for array algebra, SciPy (differential_evolution, best/1/bin) for parameter optimization, and Matplotlib for plotting.

In the optimization process, a differential evolution algorithm is employed to iteratively adjust the material parameters ( $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ ) to minimize the difference between the predicted stress values ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) and the experimental data.

This approach, known as least-squares minimization, aims to reduce the total error in the predicted stress values. Penalties are also introduced to avoid physically meaningless solutions, such as negative stress values ( $[eqn]$ ), ensuring that the parameters remain physically valid throughout the optimization process. The coefficient of determination, $[eqn]$ is used to evaluate the quality of the fit between the model’s predictions and the experimental data; it measures how much of the variance in the experimental data is captured by the model, with a value close to 1 indicating a strong fit.

The optimization analysis is conducted by implementing three different parameter bounds to demonstrate that various parameter sets can achieve high experimental fits with $[eqn]$ values close to 1, while still maintaining physical validity in order to explain the variability of parameter sets commonly reported in the literature. The goal is to analyze which parameter interval is preferable, ranging from highly restrictive, moderately restrictive, to non-restrictive bounds. Only different intervals are set for $[eqn]$ and $[eqn]$ , while maintaining for $[eqn]$ , and $[eqn]$ . The optimization process is iterated within these intervals until the set that minimizes the error is found.

Experimental data for the tensile and inflation tests’ curves are obtained from published investigations. To mitigate inconsistencies and enhance parameter stability, a stepwise optimization approach is proposed in Section 3.3. First, $[eqn]$ is optimized for small deformations dominated by isotropic behavior. Once $[eqn]$ is established, $[eqn]$ and $[eqn]$ are optimized to capture the anisotropic non-linear behavior at larger deformations, while $[eqn]$ is neglected. This sequential approach simplifies the optimization process and reduces parameter variability, addressing a major limitation in prior studies. By first fitting $[eqn]$ on the small-strain, matrix-dominated segment and then fitting $[eqn]$ and $[eqn]$ on the non-linear, fiber-dominated segment, the estimation is split into lower-dimensional subproblems, which removes many compensating parameter combinations and thereby reduces parameter variability even when $[eqn]$ remains high.

2.1. Biomechanical Modeling of Corneal Material

In hyperelastic material modeling, the local deformation is described by the deformation gradient $[eqn]$ , which maps a material point from its reference position $[eqn]$ to its current deformed position $[eqn]$ . Its determinant $[eqn]$ is the local volume ratio; for incompressible response, J = 1 (volume preservation). The Cauchy stress tensor is written as

$$[eqn]$$

where $[eqn]$ is the Cauchy stress tensor, $[eqn]$ is a Lagrange multiplier that enforces incompressibility and represents the hydrostatic/volumetric part of the stress, and I is the identity tensor (unit tensor) that produces an equal normal contribution in every direction. The symbol $[eqn]$ denotes the strain–energy density function. The tensor C = $[eqn]$ is the right Cauchy–Green deformation tensor, which measures strain in the reference configuration. The derivative $[eqn]$ is the tensor of partial derivatives of the energy with respect to C. In practice, ψ is often split into a volumetric term and a deviatoric (shape-changing) term. In this work, incompressibility is enforced by the Lagrange multiplier $[eqn]$ , so the volumetric contribution is represented by the hydrostatic term $[eqn]$ in Equation (1); therefore, ψ in Equation (1) denotes only the deviatoric strain-energy density. The deviatoric component of $[eqn]$ in the HGO model is as follows:

$$[eqn]$$

where $[eqn]$ + $[eqn]$ defines shear modulus, $[eqn]$ describes fiber stiffness, and $[eqn]$ accounts for nonlinearity in large strains. The invariants are $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ , where $[eqn]$ and $[eqn]$ are initial fiber directions.

Assuming incompressibility and following the theory of continuum mechanics, the Cauchy stress tensor can be expressed in terms of the invariants and the left Cauchy–Green deformation tensor $[eqn]$ through the following expression:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the orientation vectors in the preferred directions in the deformed configuration, and $[eqn]$ is the identity matrix.

Tensile Test

During the tensile test, the tissue is clamped at both ends, and a tensile force is applied along the length in the 1-direction, resulting in elongation characterized by the stretch ratio $[eqn]$ > 1. To maintain volume due to incompressibility, the tissue contracts in directions 2 and 3 perpendicularly. The deformation gradient tensor is given by

$$[eqn]$$

where $[eqn]$ ensures incompressibility. For a tensile load in the 1-direction, with the tissue confined to the 1–2 plane, $[eqn]$ and $[eqn]$ . The Cauchy stress tensor is expressed as follows:

$$[eqn]$$

For uniaxial tensile tests, stress is applied only in the 1-direction, so $[eqn]$   $[eqn]$ ; $[eqn]$ . Figure 2 illustrates the virtual tensile setup where the specimen is clamped at both ends and loaded along direction 1.

Tensile strips are typically cut along the inferior–superior or nasal–temporal anatomical axes; these are widely considered the principal corneal fiber directions in the central stroma of the cornea; to generalize across published datasets, this orientation is assumed for all cases. Accordingly, two preferred in-plane fiber families are defined and taken orthogonal in the reference configuration: $[eqn]$ aligned with the loading direction and $[eqn]$ aligned transversely. The local axis {1, 2, 3} is also indicated in the figure.

The schematic Figure 3 represents the stress–stretch curves used in this research for three different corneal specimens under tensile testing obtained from [13,14], highlighting varying degrees of non-linearity in their mechanical responses. Additionally, it illustrates the collagen fiber alignment during the stretching process in three stages.

2.2. Inflation Test

The inflation test is critical to understand the cornea’s biomechanics, as it closely replicates in vivo conditions by applying pressure similar to intraocular pressure (IOP). Unlike uniaxial tensile tests, it provides insights into the realistic mechanical response of corneal tissue under physiological conditions. Note that an axisymmetric thin-membrane assumption is adopted to match the experimental inflation setup and as a test simplification, implying equibiaxial loading ( $[eqn]$ ). The cornea, however, is not a perfect spherical cap and exhibits regional variation and fiber dispersion; these aspects are acknowledged in Section Limitations of the Study. The deformation gradient tensor for the inflation test is

$$[eqn]$$

where λ_1_, λ_2_, and λ_3_ are the stretches in the meridional $[eqn]$ , circumferential $[eqn]$ and radial $[eqn]$ directions, respectively. Due to incompressibility, $[eqn]$ = $[eqn]$ and $[eqn]$ = $[eqn]$ , with $[eqn]$ .

Assuming a thin-membrane model, the cornea is treated as sufficiently thin to neglect through-thickness stresses and bending effects in the radial direction $[eqn]$ , resulting in a biaxial stress state in the meridional $[eqn]$ and circumferential $[eqn]$ directions. The stress tensor is as follows:

$$[eqn]$$

where $[eqn]$ represents the stress in both directions. This setup allows studying corneal behavior under applied pressure, as shown in Figure 4.

The equations to convert pressure–apical displacement data into stress–stretch relationships are detailed in Appendix A according to [22]. These include the conversion of pressure units, calculation of the radius of curvature and stretch, and the application of Laplace’s law to determine corneal stress. The equations account for geometric changes in the cornea during deformation, assuming it behaves as a spherical shell. Figure 5 shows the transformation of the curves [10] before and after applying the normalization process.

2.3. Parameter Optimization for Tensile and Inflation Tests

The four parameters ( $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ ) are fitted simultaneously using the differential evolution algorithm (strategy best/1/bin), implemented in Python via SciPy (from scipy.optimize import differential_evolution). This algorithm efficiently handles complex, non-linear problems typical of hyperelastic material modeling. The optimizer minimizes a least-squares error between model predicted and experimental stresses, augmented with large penalties to enforce test physics: penalties avoid negative stress, and non-monotonic axial stress as stretch increases.

The algorithm runs up to 10,000 iterations with a tolerance of 1 × 10^−6^ on the objective. Parameter bounds are implemented as box constraints; their restrictiveness is varied to demonstrate that different intervals can yield different parameter sets, with further discussion provided in Section 3. Goodness of fit is reported with the coefficient of determination $[eqn]$ (where $[eqn]$ is the residual sum of squares (predicted vs. experimental), and $[eqn]$ is the total variance in experimental data).

Section 3.3 introduces a novel stepwise approach seeking to stabilize parameter identification by separating the estimation of the isotropic part ( $[eqn]$ ) from the anisotropic part ( $[eqn]$ and $[eqn]$ ).

3. Results

3.1. Parameter Optimization: Tensile Test

In the optimization process, three bound intervals are tested (Interval 1 (slightly restrictive), Interval 2 (non-restrictive), and Interval 3 (highly restrictive)) to demonstrate how bound choices affect the resulting parameter sets.

The rationale for the bounds selection is explanatory (to show how bounds affect identifiability), keeping the fiber terms effectively unrestrictive ( $[eqn]$ in [0, 100] MPa and $[eqn]$ in [0, 1000]) while varying only the isotropic bounds. Interval 1 (“slightly restrictive”) uses positive-only ranges for $[eqn]$ and $[eqn]$ (e.g., [0, 50] MPa); whether the upper limit is 50 or 100 MPa is immaterial, provided the range is sufficiently large to calculate without restriction. Interval 2 (“non-restrictive”) allows $[eqn]$ and $[eqn]$ in [−100, 100] MPa, admitting negative $[eqn]$ values, sometimes reported to compensate $[eqn]$ ; any solutions with $[eqn]$ + $[eqn]$ ≤ 0 are deemed physically inadmissible even if the statistical fit is high. Interval 3 (“highly restrictive”) adopts narrow literature-based ranges around typical $[eqn]$ and $[eqn]$ values to illustrate how tight priors can over-constrain the fit and promote parameter compensation in $[eqn]$ and $[eqn]$ .

Interval 1 (Slightly Restrictive) provides the most reliable results, with parameter bounds of [0, 50] for $[eqn]$ and $[eqn]$ , and [0, 100] and [0, 1000] for $[eqn]$ and $[eqn]$ , respectively. This interval ensures physical validity, avoiding negative shear stress values ( $[eqn]$ + $[eqn]$ ). The algorithm consistently assigns $[eqn]$ near zero, simplifying the isotropic component to a single parameter ( $[eqn]$ ), which effectively captures the material’s initial stiffness. This interval results in high $[eqn]$ values while maintaining interpretability, with parameter values reflecting the stiffness and fiber behavior across different non-linearities of the corneal tissue.

Interval 2 (Non-Restrictive), with bounds of [−100, 100] for $[eqn]$ and $[eqn]$ , leads to unrealistic results, such as negative shear stress, rendering the parameters physically invalid despite high $[eqn]$ values. This demonstrates that overly broad bounds can cause overfitting, producing results that are statistically accurate but lack meaningful representation of material properties.

Interval 3 (Highly Restrictive) tightly constrains $[eqn]$ and $[eqn]$ to [0.005, 0.05] and [−0.05, −0.005], respectively, while maintaining the same bounds for $[eqn]$ and $[eqn]$ as in Interval 1. Although this interval produces physically valid results, it forces some parameters to their boundary values, leading to compensatory adjustments in $[eqn]$ and $[eqn]$ . While useful in specific scenarios where parameter ranges are well-defined, such tight constraints may limit the optimization process and risk over-constraining the system.

Results of the estimated HGO material parameters for the optimization intervals are shown in Table 1.

Figure 6 presents the stress–stretch curves for Interval 1 (Slightly Restrictive). Plot (a) corresponds to the highly non-linear curve, (b) represents the moderately non-linear curve, and (c) shows the slightly non-linear curve. It illustrates the correlation between the experimental stress and the predicted stress. The high $[eqn]$ values (ranging from 0.97 to 0.99) demonstrate an excellent fit between the experimental and calculated curve, confirming the accuracy of the optimization process.

The three datasets fitted in this section correspond to the literature tensile curves shown previously in Figure 3 obtained from [13,14]. Stresses in directions 2 ( $[eqn]$ ) and 3 ( $[eqn]$ ) are expected to remain zero, as required by the boundary conditions of a uniaxial tensile test, confirming that the parameter set obtained from Interval 1 produces physically meaningful and reliable results across all directions.

3.2. Parameter Optimization: Inflation Test

Results from the inflation test are visually represented in Figure 5b, where three corneas (Cornea 1, Cornea 2, and Cornea 3) are displayed. These curves reflect a typical stress–stretch behavior, showing increasing stiffness as stretch progresses. It can be observed that the sudden non-linear behavior, or “knee,” occurs at different points along the curves, with some corneas showing an earlier rapid increase in stiffness, while others present a more gradual transition.

Cornea 1 exhibits a gradual increase in stress at lower strains, followed by a sharper rise later on. This indicates that Cornea 1 shows delayed fiber stiffening, as it is dominated by a lower initial stiffness, influenced mainly by a lower $[eqn]$ . Cornea 2, on the other hand, displays a more intermediate response, with a smoother balance between early-stage stiffness and later non-linear behavior. Cornea 3 presents the sharpest increase in stress almost from the beginning of the test, which suggests that Cornea 3 has a higher initial stiffness (likely due to a higher $[eqn]$ ) and reaches the non-linear fiber-dominated region earlier than the others.

As in Section 3.1, the parameter optimization for the inflation test is carried out for three different intervals, from which the material parameters are optimized using the differential evolution algorithm. A detailed analysis of the set of parameters obtained for each interval and the strong correlation with the experiments are shown in Table 2, with $[eqn]$ values close to 1.

Interval 1 (Slightly Restrictive) emerges as the most suitable for parameter estimation due to its ability to provide physically meaningful values without encountering negative or unrealistic results. However, as shown in Table 2, all curves in this interval exhibit zero values for both $[eqn]$ and $[eqn]$ , effectively neglecting the isotropic contribution to the cornea’s material behavior. This simplification suggests that the model primarily captures anisotropic stiffening due to fiber alignment, which may not fully represent the cornea’s initial stiffness, as the isotropic matrix plays a critical role in early deformation responses. This limitation highlights the need for further refinement in optimizing $[eqn]$ , as addressed in Section 3.3.

In contrast, Interval 2 (Non-Restrictive) allows for a broader exploration of parameter values but yields physically invalid results. For example, Cornea 2 shows $[eqn]$ = 2.905 MPa and $[eqn]$ = −3.792 MPa, resulting in a negative shear modulus ( $[eqn]$ + $[eqn]$ ), which is inconsistent with the physical interpretation of the corneal material. Despite achieving high $[eqn]$ values, these results lack physical coherence and cannot reliably describe the cornea’s mechanical properties.

Interval 3 (Highly Restrictive), while yielding physically valid results, imposes excessively tight constraints that limit the optimization process. Parameters like $[eqn]$ are consistently forced to boundary values, such as $[eqn]$ = 0.005 MPa, restricting the model’s ability to explore the cornea’s true mechanical behavior. Although this interval may be useful in specific scenarios, its restrictive nature makes it less suitable for general applications.

Figure 7 illustrates the stress–stretch curves for the inflation test using parameters from Interval 1, with results presented for highly non-linear (a), moderately non-linear (b), and slightly non-linear (c) cases. The model accurately captures the anisotropic mechanical behavior of the cornea under inflation conditions. In plots (a, b, c), the predicted stress aligns closely with the experimental data, achieving high $[eqn]$ values ranging from 0.9967 to 0.9993. The circumferential and meridional stresses are equivalent, consistent with the equibiaxial tension assumption, while the radial stress remains near zero, as expected under thin membrane conditions. These findings confirm that the parameter set derived from Interval 1 provides a reliable and consistent representation of the cornea’s mechanical response.

3.3. Novel Methodology for Parameter Optimization in Corneal Inflation Tests

During the optimization process for the inflation test, it became clear that the values of $[eqn]$ were unrealistically low, being calculated as zero in some configurations. This is physically incorrect because $[eqn]$ is responsible for the isotropic stiffness of the corneal matrix, which plays a crucial role in the initial response to deformation. To address this issue, a strategy is implemented to refine the model and ensure more accurate values for $[eqn]$ .

The initial portion of the stress–stretch curves, where the response is approximately linear and dominated by the isotropic matrix, is isolated and optimized separately. This step allowed the model to focus on adjusting $[eqn]$ , ensuring that the early stiffness of the cornea is accurately captured. The exact cutoff for this linear segment can vary between curves, but a practical guideline is that collagen fibers begin to bear a significant share of the load at roughly 2% strain [20]. Therefore, the 0–2% strain range can be treated as approximately linear and dominated by the isotropic matrix.

After isolating the initial phase of the curves, strict limits close to zero are imposed on $[eqn]$ , and $[eqn]$ to refine the optimization of $[eqn]$ . It is important to remember that $[eqn]$ can be neglected to simplify the material model, reducing it to a system with only three parameters to optimize. The final values for $[eqn]$ are in Table 3, demonstrating a clear increasing trend from Cornea 1 to Cornea 3. This validates the method, as the graphs show a corresponding increase in initial stiffness, aligning with the expected mechanical behavior.

Once $[eqn]$ is fixed, the anisotropic parameters $[eqn]$ and $[eqn]$ are optimized to capture the non-linear behavior of the cornea. The resulting values reflecting the influence of fiber alignment and stiffening in the later stages of the stress–strain response are shown in Table 4:

Upon analyzing $[eqn]$ , there appears to be no clear trend across the curves (it increases from Curve 1 to Curve 2 but then decreases in Curve 3). This pattern suggests that $[eqn]$ does not solely govern the non-linear behavior of the cornea and likely interacts with both $[eqn]$ and $[eqn]$ , complicating its interpretation. In contrast, $[eqn]$ shows a clear increasing trend from Curve 1 to Curve 3, playing a more dominant role in the non-linear stiffening behavior observed at higher strains. This makes sense, as $[eqn]$ is closely linked to fiber stiffening. The overlap of $[eqn]$ ’s effects with both $[eqn]$ and $[eqn]$ complicates the optimization, reinforcing the importance of first optimizing $[eqn]$ to separate isotropic behavior from fiber-driven non-linearity. Once $[eqn]$ is fixed, $[eqn]$ and $[eqn]$ can be optimized more effectively.

In summary, $[eqn]$ increases steadily from Curve 1 to Curve 3, reflecting greater initial stiffness. $[eqn]$ follows a non-monotonic pattern, with no clear trend, while $[eqn]$ consistently rises, capturing the sharp increase in stiffness due to fiber alignment at larger strains, aligning with typical non-linear stiffening seen in the cornea during inflation tests.

4. Discussion

The significant variability observed in corneal biomechanical studies, particularly in tensile and inflation tests, arises from various causes. In this research, the authors highlight that these discrepancies mainly occur because different parameter sets can still make physical sense and achieve strong correlations with experimental data. This introduces challenges in defining a universally accepted set of material parameters for corneal tissues.

In tensile tests, a wide range of parameter values is reported across different studies (see Table 5). The parameters $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ show considerable variability, reflecting the anisotropic stiffening behavior of the cornea under uniaxial loading conditions. For example, $[eqn]$ values range from 0.0025 MPa [13] to 15.11 MPa [10] (Hoeltzel et al., 1992), illustrating a great difference in the isotropic stiffness parameter. Similarly, $[eqn]$ spans from 0.0338 MPa [15] to 103.5 MPa [12] (M. Nambiar et al., 2022), showing a wide range in fiber alignment’s contribution to corneal stiffness. $[eqn]$ values vary from 5.74966 [13] to 1280.07 [15], highlighting a massive disparity in scale in the non-linear stiffening behavior of collagen fibers across the studies. These variabilities may also contribute to geometrical differences in the corneal specimens, testing different corneal regions (central vs. peripheral), variations in experimental conditions, and the methods used to model the anisotropic behavior of the tissue.

To aid comparison across sources, Table 5 (tensile) and Table 6 (inflation) use (−) to indicate values not reported. The differences across studies reflect model choices. Some studies use an HGO formulation with a two-parameter Mooney–Rivlin isotropic base ( $[eqn]$ and $[eqn]$ ), others report only a single isotropic parameter ( $[eqn]$ ), and some fit their curves with a pure two-parameter Mooney–Rivlin model without fiber terms $[eqn]$ and $[eqn]$ . For Table 6, note that [15] present a re-fit in which $[eqn]$ and $[eqn]$ are recalculated after fixing $[eqn]$ to the mean across specimens (0.004 MPa).

Inflation tests present distinct challenges. Variability in inflation studies is also very high (see Table 6). Apart from the previously mentioned possible sources of variability, an additional source arises from converting pressure–apical displacement curves into stress–stretch curves, which relies heavily on geometrical measurements such as the corneal base radius, thickness, and apex height. If these are measured inaccurately, they can drastically alter the resulting curves. Consequently, variation in material parameters in inflation tests stems from differences in these geometric inputs, in addition to model assumptions and experimental conditions. For instance, $[eqn]$ values in inflation tests range from 0.00291 MPa [15] to 0.32 MPa [10]. Similarly, $[eqn]$ values vary between 0.0217 MPa [15] and 0.055 MPa [10], and $[eqn]$ ranges from 291.28 [15] to 750 [10].

Tensile and inflation test parameters are not directly comparable. The mechanical behavior of the cornea under uniaxial tensile loading (tensile test) is characterized by higher stresses and strains compared to the inflation test, which simulates physiological conditions with lower stresses but captures the biaxial response of the cornea. Tensile test curves reach stress values up to 7.5 MPa and stretches up to 1.18, indicating a stronger response to uniaxial deformation. In contrast, inflation test curves display much lower stress levels, peaking around 0.05 MPa, with stretches reaching only 1.045. These differences explain why optimized parameters obtained from tensile and inflation tests are not interchangeable. For example, $[eqn]$ values in inflation tests are typically much larger than those in tensile tests.

Limitations of the Study

This study is subject to some limitations that must be acknowledged. First, the datasets used were limited in number and derived from previously published experimental studies, which may not fully capture the biological variability of corneal tissue. In addition, the study relied on representative smoothed stress–stretch curves from literature instead of raw experimental data, which limits the ability to evaluate the robustness of the optimization process against real experimental noise and variability.

Regarding the modeling assumptions, the inflation test was treated as a spherical shell with uniform thickness, neglecting regional variations in geometry and thickness distribution. While this simplification is common in corneal biomechanics, it inevitably introduces discrepancies with physiological reality.

A particularly relevant limitation concerns the transformation of pressure–apical displacement curves into stress–stretch curves in inflation experiments. This conversion is highly sensitive to the assumed geometric input parameters, such as corneal thickness, apex height, and base radius. Small variations or inaccuracies in these values can significantly alter the computed stress–stretch relationships, thereby affecting the optimized material parameters.

The present analysis also adopts a simplified, general model with fixed geometric inputs rather than subject-specific anatomy; this choice facilitates comparison across datasets but does not capture inter-individual variability in contour, stiffness, and thickness. Finally, the stepwise sequential approach was applied and validated only for the HGO model; although it may extend to models that clearly separate isotropic and anisotropic contributions with few parameters, its applicability to highly coupled multi-parameter laws (e.g., five-parameter Mooney–Rivlin) is limited.

To mitigate these limitations, future work may use raw experimental data rather than smoothed curves; include subject-specific geometry from imaging (apex height, base radius, thickness maps); replace the axisymmetric membrane with a biaxial, non-axisymmetric formulation; add fiber-dispersion and allow non-orthogonal preferred directions; and validate the estimated parameters on independent datasets and in finite-element simulations under physiological intraocular pressure.

5. Conclusions

The objectives of this study are to optimize the mechanical parameters of the cornea using both tensile and inflation tests and to address the variability in parameters reported across studies. It is established that different testing methods (tensile and inflation) yield significantly different results in the HGO model parameters $[eqn]$ , $[eqn]$ , and $[eqn]$ , excluding $[eqn]$ to simplify the model. To achieve this, the parameters were optimized across different optimization ranges for three tensile curves with varying degrees of nonlinearity and three inflation curves.

A key factor contributing to the variability across studies includes differences in the specific corneal zones that are tested, and initial geometric measurements like base, height, and radius in the case of inflation tests. Mathematical simplifications in virtual tests also play a role in this variability. However, it is identified and demonstrated as the primary source of variation that multiple sets of parameters can highly correlate with experimental stress–stretch curves while still making physical sense. This presents a challenge in reducing parameter variability.

This study proposes a stepwise approach to mitigate this variability. First, $[eqn]$ should be optimized based on initial deformations from tensile tests, where isotropic stiffness dominates. Once $[eqn]$ is fixed, the anisotropic parameters $[eqn]$ and $[eqn]$ , which capture fiber alignment and non-linear stiffening, may be optimized. This method simplifies the optimization process and reduces variability in parameter estimation.

Through a comparative analysis, the study demonstrates that tensile and inflation test data are not interchangeable. Each test yields different parameter values, as the stress–stretch scales differ significantly.