Experimental Determination of Corneal Elastic Constants and Their Use in Biomechanical Modeling

Abstract

Corneal biomechanics aims to establish the physico-mathematical bases that allow for predicting the corneal response to physiological and pathological situations by creating models of tissue behavior. Determining the characteristic parameters of these models is a formidable challenge in the biomechanical modeling process. To contribute to corneal tissue characterization, an experimental set-up was designed, built and tested to study corneal behavior by applying changes in pressure. The elastic constants of porcine corneas were determined, and a Young’s modulus of 0.188 MPa and 26.22% hysteresis were obtained. A computational cornea model was developed to analyze the influence of different factors. Minor variations in the applied conditions were found for apical displacement and pachymetry, and the corneal behavior was reproduced. However, the optical power behavior was affected by variations in the applied conditions, and the experimentally obtained data could not be reproduced. Despite its importance, this parameter has not been analyzed in-depth by other studies, which shows that the quality of a biomechanical cornea model should not be evaluated only by apical displacement.

1. Introduction

Nowadays, there is a growing interest in knowing the biomechanics of biological soft tissues because these materials perform nonlinear behavior that depends on the time and state of the previous load [1]. However, the complexity and diversity of its structure at both the macroscopic and microscopic levels present a challenge in developing numerical simulation models, given the desire to simulate a real structure to improve accuracy.

In clinical practice, in silico mechanical models have been developed to simulate surgical procedures in ocular tissues [2,3], muscles [4,5,6] or arteries and veins [7,8]. In all these models, the properties of materials are obtained from experimental adjustment data. All these studies coincide with the need to simplify these models while maintaining a certain level of complexity to reduce computational costs [9].

The human cornea is a transparent structure that acts as a convergent lens. Its refractive power represents about two thirds of the optical power needed to focus light on the retina. This is why different techniques seek to modify its curvature to obtain a new refractive state that would ideally remain stable over time [10].

In order to predict corneal behavior, it is extremely important to develop a numerical model capable of simulating this tissue and being personalized to patients’ specific characteristics. To develop this model, it is not only important to measure the real corneal geometry, but it is also necessary to know both the boundary conditions and corneal material elastic properties.

Nowadays, it is possible to very accurately determine both the geometry and boundary conditions from regular clinical measurements, but establishing the material properties is a harder task.

Truly representative corneal properties should take into account the inhomogeneity and viscoelastic characteristics of this tissue. However, variations in these properties are found not only among individuals but also due to changes in age, temperature, humidity, some diseases, nourishment and many other factors [11].

The need for such measurements and lack of modern equipment were the motivation to build an experimental set-up to measure the mechanical properties of a cornea submitted to different pressures and to simultaneously determine its geometry by an anterior segment scanner (Pentacam). The results of these measurements were used by a biophysical model of the human eye to predict its mechanico-optical responses [12,13,14,15]. The main limitation of this set-up is that it only works with ex vivo corneas.

2. Materials and Methods

In order to carry out inflation tests, it was necessary to design and build an experimental facility to be used in future human cornea tests.

2.1. Experimental Set-up Description

The main set-up’s components are shown in Scheme 1. The previously prepared corneal sample (1) was placed inside a pressurization chamber (2), which was kept motionless by a fixation instrument (3) in the right position to obtain the corneal topography with a Pentacam (Pentacam Scheimpflug imaging system, Oculus) (4). One of the ducts (5) from the pressurization chamber was connected to a syringe (6), which was used to control the applied pressure, while the other (7) was connected to a pressure sensor (8). The pressure sensor was connected to an Arduino UNO (9), which was controlled by a laptop (10) that also registered pressure changes.

2.1.1. Sample Collection and Preservation

Seven porcine corneas were used in the experiments. They were taken from a local slaughterhouse immediately after the animals’ deaths and preserved in a cold solution of sodium chloride (1%) for 10 h until the tests began. An ophthalmologist surgeon checked for tissue abnormalities and extracted the corneas using a section ring. The selected corneas were left at room temperature for 30 min to thermally stabilize them. Then, the corneas were fixed inside the pressurization chamber, and the tests began. The corneas were submitted to a preconditioning process of three increasing and lowering pressure cycles. The pressure was gradually built up to the maximum value to be used in the experiments. The initial pressure to start testing was maintained for 10 min in the last cycle.

2.1.2. Pressure Control

A 2.5-mL syringe was used to manually increase and decrease the pressure. To control the pressure, a differential pressure sensor (MPX 5010, 10-kPa maximum pressure value) was connected to an Arduino UNO, which was programmed to save the measurements every second on a laptop.

2.1.3. Corneal Geometry Measurements

For each deformed state corresponding to an applied pressure, geometrical data (e.g., curvature radius and pachymetry) were taken by the Pentacam system.

2.2. Biomechanical Modeling

2.2.1. Geometry

Spherical and biconical surfaces were used to model both the anterior and posterior corneal surfaces by keeping the same central thickness for each type of model employed. The Pentacam data obtained directly from the tests were used.

2.2.2. Boundary Conditions

A normal constant pressure was applied to the posterior corneal surface. Different kinds of supports were analyzed in the limbus zone, which are shown in

Table 1. Boundary conditions.

Condition Type	Scheme 1	Observations
Fixed border		Both the side and bottom edge are fixed.
Scleral ring		A 2-mm wide scleral ring embedded at its end is used.
Scleral globe (30°)		Half the scleral balloon embedded at this end is used. The joint between both materials has a 30° inclination in relation to the horizontal plane.
Scleral globe (45°)		Half the scleral balloon embedded at this end is used. The joint between both materials has a 45° inclination in relation to the horizontal plane.

¹t = corneal thickness, Ra = corneal anterior radius, Rp = corneal posterior radius and C.D = corneal diameter.

.

The scleral ring and globe radii, as well as the material assigned to the sclera, were obtained from a previous study [15]. The aim of these applied conditions was to show the influence of the constraints and geometric models on the system’s behavior. A very wide variety of options can be studied that depend on the model to be reproduced.

2.2.3. Material Model Selection and Parameter Determination

Previous works [12,15] have shown that, for hyperelastic models, a 5-parameter Mooney–Rivlin model better reproduces experimental results. The viscoelastic and anisotropic cornea properties were not analyzed because additional experimental data were lacking. However, hyperelastic models have shown good results in reproducing experimental tests [15,16,17,18].

The Mooney–Rivlin model is based on deformation invariants, and its general expression is as follows [12]:

$${\displaystyle \sum }_{\mathrm{i}+\mathrm{j}=1}^{\mathrm{N}}{\mathrm{C}}_{\mathrm{ij}}{\left({\mathrm{I}}_1-3\right)}^{\mathrm{i}}{\left({\mathrm{I}}_2-3\right)}^{\mathrm{j}}+{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{N}}\frac{1}{{\mathrm{d}}_{\mathrm{k}}}{\left(\mathrm{J}-1\right)}^{2\mathrm{k}}$$

(1)

where C_ij represents the specific material parameters, I₁ and I₂ are the Cauchy tensor invariants, d_k is a Lagrange multiplier that represents material incompressibility and J is the elastic deformation gradient determinant.

The 5-parameter Mooney–Rivlin model could be obtained using Equation (1) by having N = 2. The model parameters could be obtained by fitting the stress-strain curves in Ansys 18.1.

2.3. Parameter Definitions

2.3.1. Optical Power

The corneal optical power is one of the most important parameters of the eye functioning as an optical system. However, its behavior has not been discussed in most of the consulted references. The average optical power can be determined according to Equation (2) [19]:

$${\mathrm{P}}_{\mathrm{m}}={\mathrm{P}}_{\mathrm{m} \mathrm{ant}}+{\mathrm{P}}_{\mathrm{m} \mathrm{post}}-\frac{\mathrm{t}}{{\mathrm{n}}_{\mathrm{cornea}}}{\ast \mathrm{P}}_{\mathrm{m} \mathrm{ant}}{\ast \mathrm{P}}_{\mathrm{m} \mathrm{post}}$$

(2)

where P_{m ant} and P_{m post} are the average powers of the anterior and posterior corneal surfaces, respectively, t is corneal thickness (given in meters) and n_cornea (1.376) is the corneal refraction index.

The average powers of the anterior and posterior corneal surfaces can be determined according to [19]:

$${ \mathrm{P}}_{\mathrm{m} \mathrm{ant}}=\frac{{\mathrm{n}}_{\mathrm{cornea}}-{\mathrm{n}}_{\mathrm{air}}}{{\mathrm{R}}_{\mathrm{m} \mathrm{ant}}}$$

(3)

$${\mathrm{P}}_{\mathrm{m} \mathrm{post}}=\frac{{\mathrm{n}}_{\mathrm{aqueous}}-{\mathrm{n}}_{\mathrm{cornea}}}{{\mathrm{R}}_{\mathrm{m} \mathrm{post}}}$$

(4)

where n_air (1.000) and n_aqueous (1.336) are the refractive indices of air and the aqueous humor, respectively, and R_{m ant} and R_{m post} are the mean curvature radii of the anterior and posterior corneal surfaces, respectively (given in meters).

2.3.2. Apical Displacement

Corneal vertex displacement has been one of the most studied parameters of corneal deformation in both experimental and biomechanical modeling terms. In this work, its direct determination was not possible, which was why the shell theory was applied [20].

For a spherical shell of radius r and thickness h subjected to pressure p, it was possible to calculate the surface normal displacement w using Equation (5) [20]:

$$\mathrm{w}=\frac{{\mathrm{p}\ast \mathrm{r}}^2}{2\ast \mathrm{E}\ast \mathrm{h}}\left(1-\mathsf{\mu }\right)$$

(5)

where E is Young’s modulus and μ is Poisson’s ratio.

2.3.3. Stress-Strain

Stress-strain relations were obtained from the shell theory [20] and applied to the experimental results:

$$\mathsf{\sigma }=\frac{\mathrm{p}\ast \mathrm{r}}{2\ast \mathrm{t}}$$

(6)

$$\mathsf{\epsilon }=\frac{ \mathrm{r}\prime -\mathrm{r}}{\mathrm{r}}$$

(7)

where p is the applied pressure, t is the corneal thickness and r and r′ are the mean corneal radii before and after each pressure increment, respectively.

3. Results and Discussion

Different parameters were tested to compare the experimental and modeling results. Experimental data can be found in Supplementary materials (Tables S1–S7).

3.1. Experimental Results

3.1.1. Curvature Radii

The principal meridian radii of the anterior and posterior corneal surfaces were obtained from the Pentacam data. The variations in the average radius of each meridian on both corneal surfaces with the changes in pressure are shown in Figure 1. The variations were taken from an initial common pressure of 2.47 kPa.

Increases in pressure led to higher radii values, and the cornea became flatter with each increment in pressure and steeper with the lowering pressure values. One important fact is that the final average anterior surface radius approximately equaled the initial one, while the posterior surface had a smaller final average radius compared with the initial one. This had consequences for the optical power of the entire cornea.

Apical displacement (corneal vertex displacement) is normally used as a corneal deformation measure [21,22,23,24,25]. The authors of these papers calculated the curvature radius from the apical displacement measures and geometric simplifications and obtained lower values when the pressure increased and higher values when the pressure dropped. However, the studies that applied the Scheimpflug technology for inflation tests [26,27,28] reported similar results to those obtained herein.

3.1.2. Pachymetry

The mean pachymetry value obtained for the initially applied pressure was 885 μm (σ = 158 μm). Figure 2 shows that the pachymetry value lowered for both increasing and decreasing pressures. This behavior was not found in those studies that employed apical displacement as a corneal deformation measure, where the pachymetry was calculated with the corneal incompressibility and curvature radii (determined as previously explained). This was why those authors obtained lower corneal thickness values when the pressure increased and higher values when the applied pressure dropped.

The data pachymetry found in experiments using the Scheimpflug technology indicated two different behaviors. Lombardo et al. (2014) [27] reported narrowing pachymetry variations for both increasing and decreasing pressure values, which is similar to what we found. Kling et al. [26,28] revealed a different behavior. Changes in pachymetry were significantly less marked, with a slight increase in pachymetry when the pressure dropped, and the pachymetry had even higher values than the initial ones. Although this behavior is not completely clear, it might be related to the previous treatment to which the corneal samples were submitted. Further experiments must be carried out on this subject.

3.1.3. Optical Power

Figure 3 shows the variations in the mean corneal optical power versus the changes in pressure obtained from the experimental results. We can see how the optical power decreased when the pressure rose and increased when the pressure dropped. This value was not reported in the consulted bibliography.

The observed behavior was directly related to the changes observed in the pachymetry and the curvature radii of both the anterior and posterior corneal surfaces, as discussed above. Obviously, the changes in the hydration state and other internal changes that could affect the corneal refraction index were not considered in this model and will be studied in a future work.

3.1.4. Apical Displacement

Given the difficulties in achieving a common initial pressure for all the taken measurements, only the corneas to which it was possible to apply an average initial pressure of 1.61 kPa (12.09 mmHg) (with a 0.06-kPa standard deviation) were taken into account to obtain the pressure versus apical displacement and stress-strain relations (Corneas 1–5). Figure 4 depicts the behavior of the average obtained cornea. This behavior presented an approximately linear region between 2.6 and 5.0 kPa. In accordance with this value, the cornea rigidity slightly increased.

A wide range of reported results has been observed for this parameter [13,22,24,25,28] using porcine corneas extracted from eyeballs. The results obtained in our experiments fell within the range of the results reported in the consulted bibliography.

3.1.5. Stress-Strain

Figure 5 shows the stress-strain relations for each analyzed cornea, calculated from the average radius of each one. The stress-strain curves obtained for an increasing pressure are depicted by the solid colored lines, and those corresponding to a decreasing pressure are represented by dashed lines. For Cornea 3, only the increase is shown, as the drop in pressure could not be carried out. In all cases, greater deformation was observed with increasing stress, but it lowered when the internal pressure decreased.

The area enclosed between both curves is related to the energy lost during the experimental process, which confirmed the corneal viscoelastic character and provided us with a tissue hysteresis measure.

Figure 6 depicts the stress-strain curves obtained from the average curvature radii when the applied pressure was increased. These curves show a linear behavior between approximately 6 and 15 kPa (this range corresponds to pressures between approximately 1.60 and 3.33 kPa). The stiffness increased up to 50 kPa, and the anterior corneal surface was stiffer than the posterior one. Different authors have reported similar behavior [20,21,22,23,24].

3.1.6. Young’s Modulus Determination

Young’s modulus characterizes the material stiffness by relating the strains and stresses with Hooke’s law [20]. Although we could not use a Young’s modulus value to describe the nonlinear behavior of the corneal material, it was possible to fit the linear region of each curve shown in Figure 6 (obtained for stresses under 15 kPa) to determine this parameter for the selected range. The fitted values are found in

Table 2. Young’s modulus values for the anterior, posterior and mean corneal surfaces.

Surfaces	E (MPa)	Random Uncertainty (MPa)	R2
Anterior	0.22475	0.00152	0.9996
Posterior	0.14084	0.00428	0.9908
Mean	0.18804	0.00143	0.9995

. They are the mean values calculated to estimate the average material stiffness of the studied corneas for pressures between approximately 1.60 and 3.33 kPa.

The regression coefficients were above 0.99, which means that the linear model adequately reproduced the experimental data. The anterior surface was around 37% stiffer than the posterior one.

A wide range of Young’s modulus values have been reported, which can be explained by considering the different experimental techniques used to measure it. These results also depend on the corneal origin (e.g., human or porcine), donors’ ages, applied strain ranges, method followed to determine Young’s modulus and experimental conditions, among other factors.

Table 3. The Young’s modulus values reported for the inflation tests with porcine corneas.

Surfaces	E (MPa)	Pressure (kPa)
Elsheikh 2008 [22]	0.1449–0.2209	0–2.67
Asejczyk-Widlicka 2008 [29]	0.05–0.24	1.60–3.33
Boschetti 2012 [21]	0.0945–0.2873	0.40–4.00
Matteoli 2016 [30]	0.247	2.00–4.00
Mean cornea of this work	0.18804	1.60–3.33

shows the Young’s modulus values reported for the inflation tests with porcine corneas, which were calculated with the physiological pressure values. The variations were significantly narrower than those indicated in other experiments. The results obtained in this work fell within the range of the values reported for similar tests.

3.1.7. Corneal Hysteresis

Viscoelastic corneal hysteresis is characterized by different behaviors under load and unload conditions. It is a measure of the energy absorbed during the cycle and represents a combination of this tissue’s elastic and viscous properties. It is useful for understanding corneal resistance to intraocular pressure [31].

Hysteresis can be determined as the ratio of the area contained in the cycle divided by the area under the load curve. The corneal hysteresis values obtained in this work varied between 15.64% and 44.21%, with an average value of 26.22% and a standard deviation of 13.36%.

3.2. Biomechanical Modeling

The geometric and material corneal computational model characteristics were obtained from the experimental results as explained before. The used data can be found in Appendix A and Appendix B (

Table A1. Corneal geometrical parameters.

Type of Surfaces	Rfa (mm)	Rsa (mm)	Qfa	Qsa	Rfp (mm)	Rsp (mm)	t (µm)	C.D (mm)
Biconic	8.03	7.49	0	−0.36	6.71	6.02	816	12.43
Spheric	7.76		0		6.37		816	12.43

Rfa = flattest radius of the anterior corneal surface; Rsa = steepest radius of the anterior corneal surface; Qfa = flattest asphericity parameter of the anterior corneal surface; Qsa = steepest asphericity parameter of the anterior corneal surface; Rsp = steepest radius of the posterior corneal surface; Rfp = flattest radius of the posterior corneal surface; t = corneal thickness; C.D = corneal diameter; The Pentacam provides only the asphericity parameters of the anterior corneal surface.

and

Table A2. The applied Mooney–Rivlin material model parameters.

Structure	C10 (MPa)	C01 (MPa)	C20 (MPa)	C11 (MPa)	C02 (MPa)
Cornea	1.19	−1.11	17.84	−26.83	10.59
Sclera	34.73	−32.43	75,029	−153,636	78,853

). The geometry influence and comparisons between the experimental and computational results are shown for several parameters.

3.2.1. Pachymetry

Influence of the Geometry

In Figure 7, a reduction in the pachymetry is observed when the applied pressure was increased in both models (spherical and biconical geometries). This behavior corresponded to that obtained experimentally. There is a coincidence in the values of these models, except in the points closest to the maximum pressure (8.00 kPa), but they were still similar.

Influence of the Boundary Conditions

The pachymetry variation versus the pressure for the different applied boundary conditions is shown in Figure 8.

Small differences were observed in the results obtained among the used models, and the bigger ones were less than 10 μm. Pachymetry is a parameter that is barely affected by both the use of different geometric models and the analyzed boundary conditions.

Model–Experiment Comparison

It was necessary to obtain a stress-free cornea configuration [12] to compare the experimental and computational results, in which the displacements could be measured by starting at the initially applied pressure (corresponding to the initial geometry). The results are provided in Figure 9. Similar behaviors were noted, and the differences were less than 10 μm, representing <1.5%.

3.2.2. Apical Displacement

Influence of the Geometry

Figure 10 shows the apical displacements obtained for the two geometric models. Increasing the applied pressure led to greater apical displacement, as observed in the experiments. A complete coincidence between both models occurred, which means that the geometry used in these cases did not seem to affect this parameter’s behavior.

Influence of the Boundary Conditions

Figure 11 shows the apical displacement behavior due to increases in pressure for different supports.

No differences were observed between the model with a fixed border and the scleral ring one, but the results differed for the other conditions. The scleral globe model had the smallest apical displacements. This could be due to the sclera assuming part of the deformation caused by the applied pressure. The scleral globe model with a 30° inclination on the surface that separated the cornea and sclera presented fewer deformations than the scleral globe model at 45°. This means that employing such a model requires analyzing the interaction between both tissues if the intention is to model corneal behavior inside the eyeball.

Model–Experiment Comparison

Figure 12 illustrates the apical displacements obtained for the computational model (black squares) and the experimental results (red circles).

Although real experimental conditions are hard to reproduce with a numerical model, the obtained results indicated a marked resemblance. Studying the different boundary conditions and material properties could lead to a model that accurately reproduced reality.

3.2.3. Optical Power

Influence of the Geometry

The optical power variations with pressure are shown in Figure 13 for the models that employed spherical and biconical surfaces.

The optical power behavior obtained from computational modeling presented significant differences compared with the experimental results (see Figure 3). In the former, an increase in optical power for the pressures took place, which rose to 4.00 kPa, plus lower values for an ulterior pressure increase. In the latter, a decreasing monotone behavior was noted.

Major differences were observed between the spherical and biconical models, which were around 2 D within the studied pressure range.

Influence of the Boundary Conditions

Figure 14 shows the optical power variation with the pressure for the different employed boundary conditions.

Similar behavior occurred to that obtained in Figure 13. The maximum optical power value obtained was at 4.67 kPa in the model with a fixed border and 4.00 kPa in the remaining cases. The biggest differences were observed for 4.00 kPa. The models that used the eyeball came closer to one another and presented bigger differences, falling between 3.33 and 4.67 kPa. The fixed border model presented the highest optical power values and the widest variations, which was related to the applied support, and forced the material to undergo more marked deformations. These two deformation types were lesser for the scleral ring model and the eyeball model.

In these cases, an optical power analysis showed that the differences found among the models led to a much more in-depth study of the real cornea conditions, which was essential if the cornea behavior was to be reproduced. It was also very interesting to see how the used models very accurately reproduced the experimental behavior of some parameters, such as apical displacement and pachymetry, but that of the optical power of the cornea differed. This result indicates that determining the apical displacement or pachymetry is insufficient to verify the quality with which a computational model reproduces the experimental results.

Reported Results

The optical power variation from computational modeling has only been reported in five of the consulted references [32,33,34,35,36]. Their authors employed the material and geometrical parameters taken from the experimental results reported by other authors and did not directly determine the curvature radii. According to these results, the change in the optical power monotony versus the pressure curves obtained in our simulations would be the expected result, but the causes of this behavior and the conditions that could influence it were unclear. A more in-depth study into the influence of mathematical modeling is necessary to characterize the material properties and to analyze different boundary conditions. To better understand this tissue behavior, more experimental results are necessary, in which the curvature radii could be measured for each pressure increment.

Our laboratory lacks an automatic pump with the appropriate characteristics to carry out such tests, which is why we resorted to manual insertion. However, tests were carried out on numerous occasions to establish a robust sensitive protocol with the experiment to be run.

4. Conclusions

This study obtained the elastic properties of porcine corneas from an experimental set-up designed and built to be used together with topographic equipment. Geometric cornea characterization was performed under controlled pressure conditions. The Young’s modulus (0.188 MPa) and hysteresis (26.22%) were determined, as was a five-parameter material Mooney–Rivlin model. Variations in pachymetry, apical displacement and optical power were analyzed for both the experimental and modeling results. The effect of the different types of geometries and boundary conditions were studied, with very little influence on the variation in the pachymetry and apical displacement but significant differences in optical power. Although the pachymetry and apical displacement results were similar in both the experimental and computational cases, the optical power showed significant differences in them. Very little data about optical power have been reported, and those that exist all show a similar behavior to that obtained from the modeling in the present work. A more in-depth study about the conditions that influence this parameter’s behavior is necessary. We conclude that apical displacement and pachymetry determination are not enough to test the quality of a computational model for reproducing the experimental results.

For future research, more complex material models and other optical and mechanical parameters can be analyzed for greater accuracy to help surgeons in their clinical practice.