Individualized Characterization of the Distribution of Collagen Fibril Dispersion Using Optical Aberrations of the Cornea for Biomechanical Models

Abstract

Purpose: The spatial distribution of collagen fibril dispersion has a significant impact on both corneal biomechanical and optical behaviors. The goal of this study was to demonstrate a novel method to characterize collagen fibril dispersion using intraocular pressure (IOP)-induced changes in corneal optical aberrations for individualized finite-element (FE) modeling.

Methods: The method was tested through both numerical simulations and ex vivo experiments. Inflation tests were simulated in FE models with three assumed patterns of collagen fibril dispersion and experimentally on three rhesus monkey corneas. Geometry, matrix stiffness, and the IOP-induced changes in wavefront aberrations were measured, and the collagen fibril dispersion was characterized. An individualized corneal model with customized collagen fibril dispersion was developed, and the estimated optical aberrations were compared with the measured data.

Results: For the theoretical investigations, three assumed distributions of fibril dispersion were all successfully characterized. The estimated optical aberrations closely matched the measured data, with average root-mean-square (RMS) differences of 0.29, 0.24, and 0.10 µm for the three patterns, respectively. The overall features of the IOP-induced changes in optical aberrations were estimated for two ex vivo monkey corneas, with average RMS differences of 0.57 and 0.43 µm. Characterization of the fibril dispersion in the third cornea might have been affected by corneal hydration, resulting in an increased RMS difference, 0.8 µm.

Conclusions: A more advanced corneal model with individualized distribution of collagen fibril dispersion can be developed and used to improve our ability to understand both biomechanical and optical behaviors of the cornea.

Finite-element (FE) modeling offers the potential to provide insight on corneal biomechanics, including the impact on postoperative optical outcomes.^1–3 The biomechanical behavior of the cornea is mainly determined by the material structure of the stroma, which is composed of hundreds of layers of collagen fibrils (lamellae) embedded in a matrix of proteoglycans. Stiffness of the extrafibrillar matrix varies through the thickness, with anterior stroma being stiffer than posterior.^4,5 It has been shown that this depth dependence has a significant impact on the biomechanical behavior of intact and surgically treated corneas.^6,7 The depth-dependent matrix stiffness of individual corneas can be measured through ex vivo mechanical tests, such as the shear test.^4,5 In vivo measurement is always more challenging. Although recent studies of Brillouin microscopy have shown some potentials for in vivo measurement,⁸ the measurement result strongly depends on cornea hydration,⁹ and the correlation between the measured Brillouin frequency shift and the real mechanical modulus of the cornea still requires further investigation. 

Compared to the extrafibrillar matrix, the microstructure of the collagen fibrils is more complex. Previous x-ray diffraction studies showed that variations exist in the collagen fibril orientation among different species.^10–12 In the human cornea, the collagen fibrils are regularly packed in collagen lamella. The stroma consists of hundreds of layers of collagen lamellae, with collagen fibrils lying in all directions within the cornea plane.^13,14 Further quantitative analyses have shown that collagen fibrils are more aligned in two orthogonal preferential orientations—nasal–temporal (NT) and superior–inferior (SI)—at the center of the cornea and one circumferential preferential orientation at the limbus.^15–17 In addition, the degree of fibril dispersion exhibits a spatial variation over the corneal surface,^15,18 affecting the biomechanical behaviors of the cornea, including the degree of anisotropy, apical displacement, and stress distribution.³ Our previous sensitivity study identified fibril dispersion as an important material parameter that influences corneal aberrations induced by IOP elevation.¹⁹ The spatial distribution of fibril dispersion may also vary among individual healthy corneas and between healthy and diseased corneas, such as keratoconus.^15,20,21 Although biomechanical models accounting for the details of collagen fibril structure have been proposed in several studies,^2,3,22 inter-subject variations in fibril dispersion have not been studied yet. Capturing this variation in individualized models has been challenging due to limitations in the currently available techniques. Several different imaging techniques, such as electron microscopy,^13,14 second harmonic generation imaging,^23,24 and x-ray diffraction,^18,25,26 have been used to investigate the microstructure of corneal collagen fibrils; however, there is still a lack of methods to directly quantify the mechanical properties of the collagen fibers from the imaging outputs and transfer these to individualized modeling. Ex vivo mechanical tests such strip extensometry tests^27–30 and inflation tests^31–33 can measure the highly nonlinear stress–strain behavior of the cornea, but independent measurement of the modulus of extrafibrillar matrix and collagen fibers is very difficult because the collagen fibrils in corneal stroma are always embedded in the extrafibrillar matrix. To develop an individualized biomechanical model of the cornea that improves the prediction of the optical outcomes, a new method to characterize the regional distribution of the collagen fibril dispersion is needed. 

Our preliminary study suggested that aberration changes are approximately correlated with degree of fibril dispersion, as well as the IOP elevation (see Discussion for more details). Based on this finding, we proposed a novel method in the present study to characterize the spatial distribution of collagen fibril dispersion of individual corneas using the intraocular pressure (IOP)-induced changes in corneal optical aberrations. The method was tested through both numerical simulations and ex vivo experiments on rhesus monkey corneas. First, the changes in wavefront aberrations induced by IOP elevation were obtained for both simulated and experimentally measured corneas. Then, with additional information regarding measurable biomechanical parameters, corneal geometry, and depth-dependent matrix stiffness, the distribution of the collagen fibril dispersion was estimated based on proposed mathematical equations. Finally, three-dimensional (3D) individualized FE models of the cornea with customized geometry, depth-dependent matrix stiffness, and collagen fibril dispersion were developed for each cornea, and the final predictions of corneal optical aberrations were compared with the experimental data. 

Methods

FE Model of the Cornea

The 3D anisotropic hyperelastic corneal models in this study were all developed using FE software (Abaqus Unified FEA; Dassault Systèmes Simulia Corp., Johnston, RI, USA) as described in our previous study¹⁹ and shown in Figure 1. An improved strain energy function from the original Gasser–Holzapfel–Ogden model³⁴ with pre-integrated fibril distribution was adopted to describe the microstructure of the cornea. To be consistent with the boundary condition of the cornea in the ex vivo experiment, the peripheral edge of the corneal model (limbus) was fixed in displacement. Four main material parameters describe the material behavior of the cornea: matrix stiffness, C₁₀; fiber stiffness, k₁; fiber nonlinearity, k₂; and fibril dispersion, κ(ρ, θ), which is a function of the distance from the optical axis, ρ, and the angle of any meridian with the horizontal axis, θ. Two fiber families are considered, with preferential orientations progressively varying from two orthogonal NT and SI orientations at the center of the cornea to one circumferential orientation at the limbus.^15,16 The parameter κ(ρ, θ) represents a continuous distribution of the collagen fibrils over the entire cornea. It has a minimum value of 0 where fibrils are perfectly aligned along their preferential orientations and a maximum value of 1/3 where fibrils are fully dispersed and the cornea behaves nearly isotropically. 

Figure 1.

 

Cross-sectional view of an example of the 3D FE model of the cornea, where the color map represents the spatially varied fibril dispersion parameter κ(ρ, θ). The range of κ(ρ, θ) is shown in the color bar; smaller values indicate more aligned fibrils, and larger values indicate more dispersed fibrils.

Correlation Between Fibril Dispersion and IOP-Induced Changes in Corneal Wavefront Aberrations

We proposed mathematical equations to characterize the fibril dispersion parameter κ(ρ, θ) in our 3D anisotropic corneal model from the IOP-induced changes in corneal wavefront aberrations. Aberrations of the cornea can be represented using Zernike circle polynomials, which are also functions of normalized radial distance ρ and azimuthal angle θ, defined as W(ρ, θ). The total changes in corneal wavefront aberrations induced by IOP elevation, W(ρ, θ)\(_{IO{P_2} - IO{P_1}}\), have five contributing biomechanical parameters, as indicated in Equation 1:  

\begin{eqnarray}W\!\left({{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{IO{P_2} - IO{P_1}}} &=& \big[ W{{\left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}_{geometry}} \nonumber\\ &&+ W{{\left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}_{{C_{10}}}} + W{{\left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}_{\kappa \left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}} \nonumber\\ &&+ W{{\left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}_{{k_1}}} + W{{\left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}_{{k_2}}} \big]_{IO{P_2} - IO{P_1}}\end{eqnarray}

(1)

where the subscripts geometry, C₁₀, κ(ρ, θ), k₁, and k₂ represent the aberration changes due to those parameters. Because variations in fiber stiffness (k₁) and fiber nonlinearity (k₂) between the eyes have nearly negligible impact on the optical aberrations,¹⁹ Equation 1 can be approximately simplified as:  

\begin{eqnarray}W\!\left({{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{IO{P_2} - IO{P_1}}} &=& \big[W\left( {{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{geometry}} + W\left( {{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{{C_{10}}}}\nonumber\\ &&+ W\left( {{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{\kappa \left( {{\rm{\rho }},{\rm{\theta }}} \right)}}\big]_{IO{P_2} - IO{P_1}}\end{eqnarray}

(2)

 

From Equation 2, the relationship between fibril dispersion κ(ρ, θ) and the changes in corneal wavefront aberrations can be expressed as follows:  

\begin{eqnarray}{\left[W\!\left( {{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{\kappa \left( {{\rm{\rho }},{\rm{\;\theta }}} \right)}}\right]_{IO{P_2} - IO{P_1}}} &=& W\left( {{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{IO{P_2} - IO{P_1}}} \nonumber\\ &&- \big[W\left( {{\rm{\rho }},{\rm{\;\theta }}} \right){{\rm{\;}}_{geometry}} \nonumber\\ &&+ W{\left( {{\rm{\rho }},{\rm{\;\theta }}} \right)_{{C_{10}}}}\big]_{IO{P_2} - IO{P_1}}\end{eqnarray}

(3)

 

Finally, we hypothesized that there is an approximately linear correlation between fibril dispersion κ(ρ, θ) and \({[W{( {{\rm{\rho }},{\rm{\theta }}} )_{\kappa ( {\rho ,\theta } )}}]_{IO{P_2} - IO{P_1}}}\) based on our preliminary results. To estimate the magnitudes of κ(ρ, θ), we needed a scaling factor A (with units of µm^–1) to first convert the wavefront into the full range of κ(ρ, θ) (0–0.333), and then we used an optimization coefficient (C) to find the optimal range of κ(ρ, θ) by varying only the maximum value of κ(ρ, θ) and fixing the minimum value to be 0:  

\begin{equation}{\rm{\kappa }}\left( {{\rm{\rho }},{\rm{\theta }}} \right) = C \cdot \left\{ {A \cdot \left\{ {{{\left[W{{\left( {{\rm{\rho }},{\rm{\theta }}} \right)}_{{\rm{\kappa }}\left( {{\rm{\rho }},{\rm{\theta }}} \right)}}\right]}_{IO{P_2} - IO{P_1}}}} \right\}} \right\}\end{equation}

(4)

 

Because there are techniques available for individual measurement of corneal geometry,^35,36 as well as the depth-dependent matrix stiffness (C₁₀),^4,5 the term \({[W( {{\rm{\rho }},{\rm{\theta }}} ){{\rm{\;}}_{geometry}} + W{( {{\rm{\rho }},{\rm{\theta }}} )_{{C_{10}}}}]_{IO{P_2} - IO{P_1}}}\) in Equation 3 can be quantified by fitting the measured geometry and matrix stiffness (C₁₀) into the FE model of the cornea. To characterize fibril dispersion κ(ρ, θ), the total changes in corneal wavefront aberrations \(W( {{\rm{\rho }},{\rm{\theta }}} ){{\rm{\;}}_{IO{P_2} - IO{P_1}}}\) need to be measured. Therefore, numerical simulations and ex vivo inflation tests were carried out to measure \(W( {{\rm{\rho }},{\rm{\theta }}} ){{\rm{\;}}_{IO{P_2} - IO{P_1}}}\) for individual corneas for the quantification of κ(ρ, θ). 

Numerical Simulations

We first demonstrated the feasibility of the proposed method with numerical simulations in our 3D FE corneal models. The FE models of three corneas were developed for a simulated inflation test. To focus on the characterization of the distribution of fibril dispersion κ(ρ, θ), we did not vary other biomechanical parameters and used average values from previous literature.^4,31,37 The same simplified axially symmetric corneal geometry (0.572-mm central thickness, 0.779-mm peripheral thickness, and 7.97-mm radius of curvature of anterior surface) was used for the three models.³¹ The preferential orientations of the collagen fibrils were also assumed to be known and modeled based on previous studies.^15,16 For depth-dependent matrix stiffness C₁₀, five material layers with uniformly distributed central thickness from anterior to posterior corneal surface (C_{10_a} = 23.1 kPa, C_{10_b} = 14.0 kPa, C_{10_c} = 5.5 kPa, C_{10_d} = 2.8 kPa, C_{10_e} = 1.8 kPa) were assigned in each model based on ex vivo experimental data of corneal shear modulus.⁴ Material parameter k₁ and k₂ were set to average values of 91.55 kPa and 785.68, respectively.³⁷ 

A blind test with three steps was designed to test the characterization of the distribution of fibril dispersion of each cornea. Specifically, the value of fibril dispersion at each node position in corneal plane was quantified. The minimum value of fibril dispersion value across a single cornea was always 0, representing strongly aligned collagen fibrils at limbus, whereas the maximum value varied depending on the specific collagen fibril structure of the cornea. In step 1, investigator A selected three different patterns of fibril dispersion κ(ρ, θ) to generate widely varying corneal wavefront aberrations (Fig. 2). To represent physiological loading, IOP elevations of 20, 30, and 40 mm Hg were applied uniformly to the posterior corneal surface and used to generate the simulated experimental data of changes in corneal wavefront aberrations. The 3D anterior and posterior corneal surface profiles were exported from each corneal model at each level of IOP elevation to compute the changes in wavefront aberrations for the central 6-mm corneal diameter using a custom-developed MATLAB (MathWorks, Natick MA, USA) program. To estimate the fibril dispersion κ(ρ, θ) of the entire cornea, the wavefront data from the central region were then expanded to a full corneal diameter of 11.25 mm using a customized linear extrapolation algorithm. In step 2, the information of corneal geometry, material parameter C₁₀, k₁, k₂, and the simulated results of the IOP-induced changes in wavefront aberrations were passed to investigator B as known parameters to estimate the specific κ(ρ, θ) of each cornea using Equations 1 to 4. The accuracy of the predicted values was evaluated in step 3. By fitting the estimated κ(ρ, θ) into the individualized FE model developed by investigator B, the predicted changes in corneal wavefront aberrations for the central 6-mm cornea induced by IOP elevation were compared with the simulated experimental data obtained from investigator A. 

Figure 2.

 

Three different assumed fibril dispersion patterns that generated different optical aberrations in numerical simulations. (a) Pattern 1 in the FE model (left); wavefront map of aberrations for the central 6-mm cornea composed of the Zernike aberration terms defocus (Z₄), astigmatism (Z₅), spherical aberration (Z₁₂), and quadrafoil (Z₁₄) (right). (b) Pattern 2 in the FE model (left); wavefront map of aberrations for the central 6-mm cornea composed of Zernike aberration terms defocus (Z₄), astigmatism (Z₅), trefoil (Z₆ and Z₉), coma (Z₇ and Z₈), spherical aberration (Z₁₂), and quadrafoil (Z₁₄) (right). (c) Pattern 3 in the FE model (left); wavefront map of aberrations for the central 6-mm cornea composed of Zernike aberration terms defocus (Z₄), astigmatism (Z₅), trefoil (Z₆ and Z₉), coma (Z₇ and Z₈), spherical aberration (Z₁₂), and quadrafoil (Z₁₄) (right).

Ex Vivo Experiment

Measurement of Optical Aberrations

After the initial verification through numerical simulations, the characterization method of fibril dispersion was further tested through ex vivo experiments. The IOP-induced changes in corneal wavefront aberrations were measured on six rhesus monkey corneas (male; five of them were 156 days old, one was 10 years old) in an ex vivo inflation test. The measurements on monkey corneas were conducted based on the University of Rochester and University of Houston's Institutional Review Board/Ethics Committee-approved protocols. Corneas were harvested immediately after animal death with approximately 3 to 4 mm of scleral rim attached and the NT orientation marked (Fig. 3a). They were stored in corneal storage medium (Optisol-GS; Bausch & Lomb, Bridgewater, NJ, USA) at 37°F until the inflation tests, which were all performed within 24 to 36 hours of harvest. 

Figure 3.

 

(a) Corneal specimen with scleral rim attached and the marker for NT orientation. (b) Illustration of the ex vivo experimental setup for the corneal inflation test and wavefront measurement.

The experimental setup for measuring corneal wavefront aberrations with increasing IOP is shown in Figure 3b. Each cornea was mounted in a specially designed three-part wet cell with aligned orientations. A middle clamp was used to provide watertight fixation along the sclera rim. One side of the bottom chamber was connected to a water column filled with Optisol-GS to increase the IOP and the other side of the bottom chamber was connected to a digital pressure gauge (Dwyer Instruments, Inc., Michigan City, IN, USA) to monitor the variation of IOP. The top chamber of the wet cell was also fully filled with a small amount of Optisol-GS during the test to control the hydration state of the cornea and prevent evaporative effects on the measurements. For optical measurements, both the top and the bottom chambers had a central optical window (12.5 mm in diameter) that allowed passage of the laser beam, and a Shack–Hartmann wavefront sensor was used to measure wavefront aberrations of the cornea. From the baseline IOP, where the shape of the cornea was fully maintained, IOP elevations of 5, 15, and 25 mm Hg were applied. Before collection of the wavefront data, a 10-minute preloading process with an IOP elevation of 25 mm Hg was carried out for each cornea. The entire inflation test time was 30 to 40 minutes for one cornea. The measurement temperature was controlled at lab room temperature (73°F) set by a centralized system at the university. The wavefront aberrations were all analyzed for the central 6-mm diameter region of the cornea and then extrapolated to a full corneal diameter size using the same method as numerical simulations. 

Measurement of Corneal Geometry

For the measurement of individual corneal geometry, the anterior corneal surface topography was first measured in vivo using a corneal topographer (EyeSys Vision, Inc., Houston, TX, USA) before animal sacrifice to obtain the steepest and flattest corneal radius of curvature and the corresponding meridional angles. After wavefront measurement, the top chamber of the wet cell was removed, and the IOP was lowered back to the baseline. The geometry of each cornea was measured again using a commercial optical coherence tomography (OCT) system (Spectralis HRA+OCT; Heidelberg Engineering, Heidelberg, Germany). During the measurement process, Optisol-GS was added on top of the cornea every 2 to 3 minutes; the tip of a glass pipette was gently guided over the corneal surface to create a uniform Optisol layer to keep it hydrated and to avoid any artifacts in the imaged corneal surface profile. Cross-sectional images along eight different meridians (0° to 180°, with a 22.5° step) were obtained. Corneal diameter and the radius of curvatures of both anterior and posterior surfaces along different meridians were analyzed using a custom-developed MATLAB program. By combining the in vivo and ex vivo measurements, geometrical information along eight to 10 different meridians was obtained for each cornea. To develop the geometrical model, curvilinear edges representing anterior and posterior corneal surfaces were sketched along each meridian, and multiple edges that formed the shape of a corneal surface were then blended to obtain a 3D surface for both anterior and posterior cornea. 

Measurement of Depth-Dependent Shear Modulus

The depth-dependent shear modulus of individual corneas was measured via an ex vivo dynamic shear test within 24 hours after inflation testing and OCT measurement. Corneas were stored in Optisol-GS at 37°F until testing was performed. For the shear test, 3-mm-diameter corneal tissue buttons were dissected, fluorescently stained, and mechanically tested in a tissue deformation imaging stage (TDIS) as described in detail in previous studies.^4,38–41 Briefly, a microscope-mounted mechanical device was used, and the stained specimen was held in place between two flat armatures. The stationary armature was connected to a calibrated spring to monitor the external shear force imparted in the specimen, and a moving armature was connected to a closed-loop piezoelectric actuator that allowed control of the applied shear deformation. The TDIS was outfitted with a retainer cup that was filled with PBS to maintain specimen hydration.³⁸ Prior to shear testing, the specimen was also axially compressed until it reached the initial corneal thickness measured with OCT from the geometry measurement. Specimen thickness was equal to the gap between the armatures. The movement of fluorescently labeled stromal keratocytes were tracked under the shear force for each cornea. 

The shear stress and strain, τ₀ and γ₀, were calculated using a particle-tracking algorithm³⁸ and were used to calculate the depth-dependent dynamic shear modulus amplitude as a function of stromal depth according to |G*(z/T)| = τ₀/γ₀. 

Development of the Individualized Corneal Model

The steps in the ex vivo experiment for material characterization and how the measurements were used for individualized corneal modeling are summarized in Figure 4. The geometrical model of each cornea was developed based on the geometrical data obtained from the in vivo topography measurement and the ex vivo OCT measurements. Five material layers with different matrix stiffness (C₁₀) were defined, with values calculated from the ex vivo shear test based on Equation 5:  

\begin{equation}{C_{10\_i,\;\;\;i = a - e}} = {G^*}\left( {z/T} \right)/2\end{equation}

(5)

 

Figure 4.

 

Steps in ex vivo experiment for material characterization and how the measurements were used in developing an individualized corneal model.

Average values of 91.55 kPa and 785.68³⁷ were assumed for material parameters k₁ and k₂. There are no published data regarding the collagen fibril orientation in rhesus monkey corneas, but because rhesus monkeys and human are closely related species the preferential orientations of the collagen fibrils in the monkey cornea model were assumed to be the same as the human cornea model in previous numerical simulations. Contributions from corneal geometry and the known material parameters to the total changes in corneal wavefront aberrations were first quantified though FE modeling. The unknown material parameter κ(ρ, θ) of each cornea was then estimated from the changes in wavefront aberrations of the entire corneal diameter after extrapolation based on Equations 1 to 4 to create a complete individualized model for each ex vivo cornea. The predicted total changes in wavefront aberrations induced by IOP elevation for 6-mm corneal diameter from the model were then compared with the experimental data. 

Results

Numerical Simulation

The estimated distribution of κ(ρ, θ) from investigator B based on the simulated wavefront data for the three corneal models was compared with the initially assumed pattern of κ(ρ, θ) in Figure 5. For all of the three simulated corneas, the magnitudes of κ(ρ, θ) were slightly underestimated. The difference between the assumed and estimated magnitudes of κ(ρ, θ) over the entire cornea had RMS values of 0.07, 0.08, and 0.08 for patterns 1, 2 and 3, respectively, corresponding to 20.7%, 22.9%, and 22.9% error of the overall range of kappa. However, the estimated patterns of fibril dispersion maintained the main characteristics of the assumed patterns. With the estimated κ(ρ, θ) values, the predicted changes in corneal wavefront aberrations induced by IOP elevation from the final FE models for all of the three distribution patterns matched the simulated experimental data well. An example of the comparison between the changes in wavefront aberrations for the three corneal models at an IOP elevation of 40 mm Hg is shown in Figure 6, which shows nearly identical wavefront maps from the simulated experiment and the FE model prediction. In addition, predicted magnitudes of the change in each Zernike coefficient from the FE model were all very close to the experimental data. For estimated fibril dispersion pattern 1, the differences between the model prediction and the experimental data of the changes in all of the Zernike coefficients had RMS values of 0.22, 0.26, and 0.40 µm for IOP elevations of 20, 30, and 40 mm Hg, respectively. The differences were smaller for the corneal model with pattern 2, with RMS values of 0.11, 0.26, and 0.34 µm for the three IOP elevations, respectively. Estimated fibril dispersion pattern 3 gave the most accurate predictions of wavefront aberrations, with RMS values of 0.07, 0.12, and 0.15 µm for the three IOP elevations, respectively. 

Figure 5.

 

Results of numerical simulations. (a) Comparison of the assumed and estimated pattern 1 of the fibril dispersion parameter κ(ρ, θ). (b) Comparison of the assumed and estimated pattern 2 of the fibril dispersion parameter κ(ρ, θ). (c) Comparison of the assumed and estimated pattern 3 of the fibril dispersion parameter κ(ρ, θ).

Figure 6.

 

Comparisons of the predicted changes in corneal wavefront aberrations for the central 6 mm of the FE models with the simulated experimental data at an IOP elevation of 40 mm Hg. (a) FE prediction and experimental data for the cornea with fibril dispersion pattern 1. (b) FE prediction and experimental data for the cornea with fibril dispersion pattern 2. (c) FE prediction and experimental data for the cornea with fibril dispersion pattern 3.

Ex Vivo Experiment

Due to the challenges associated with collecting a complete dataset from three consecutive measurements for each cornea, including wavefront measurement, geometry measurement, and matrix stiffness measurement, a complete FE model with customized corneal geometry, depth-dependent C₁₀, and fibril dispersion κ(ρ, θ) was successfully developed for only three of the six rhesus monkey corneas tested. The measured depth-dependent C₁₀ of the three ex vivo monkey corneas are shown in Figure 7, and Figure 8 shows the estimated distributions of fibril dispersion κ(ρ, θ) based on the changes in corneal wavefront aberrations at the highest IOP elevation, 25 mm Hg. Cornea 1 and cornea 2 are the left and right eyes of the same animal (156 days old) and cornea 3 is from a different animal (10 years old). 

Figure 7.

 

A comparison of the depth-dependent matrix stiffness for the three ex vivo monkey corneas used for finite element model generation. Corneas 1 and 2 were paired corneas from one young animal (156 days old) and cornea 3 was from an older animal (10 years old).

Figure 8.

 

The estimated distributions of collagen fibril dispersion for three ex vivo rhesus monkey corneas.

The predicted IOP-induced changes in corneal wavefront aberrations from the individualized FE model were also compared to the measured wavefront data for the central 6-mm cornea (Fig. 9). The wavefront maps representing the IOP-induced changes predicted by the FE model of cornea 1 and cornea 3 (Figs. 9a, 9c) maintained the overall features of the wavefront map of the experimental data. The individualized model of cornea 1 (Fig. 9a) predicted the overall astigmatism-like shape of the wavefront in the experimental data. The changes in most of the main Zernike coefficients were predicted, including asymmetric aberration astigmatism (Z₃ and Z₅), trefoil (Z₆ and Z₉), coma (Z₇), and quadrafoil (Z₁₀ and Z₁₄). The largest difference was observed in radially symmetric aberration defocus (Z₄), with a value of 0.31 µm. The RMS differences for all of the Zernike coefficients between the model prediction and the experimental data were 0.68, 0.51, and 0.51 µm for IOP elevations of 5, 15, and 25 mm Hg, respectively. Similarly, the individualized model of cornea 3 (Fig. 9c) predicted the overall trefoil-like shape of the wavefront in the experimental data. The magnitudes of IOP-induced changes for cornea 3 are much smaller compared to cornea 1, which may due to stiffening with age. The predicted changes in radially symmetric aberration defocus (Z₄) and asymmetric aberration astigmatism (Z₃) trefoil (Z₆ and Z₉), coma (Z₇), and secondary astigmatism (Z₁₁) were still close to the experimental data. 

Figure 9.

 

Comparisons of the predicted changes in corneal wavefront aberrations for the central 6 mm from the individualized monkey corneal model and the experimental data. (a) FE prediction and experimental data for changes in wavefront aberrations of monkey cornea 1. (b) FE prediction and experimental data for changes in wavefront aberrations of monkey cornea 2. (c) FE prediction and experimental data for changes in wavefront aberrations of monkey cornea 3. The scale in the aberration plots in c is different from a and c due to the relatively small values of the changes in aberrations of cornea 3.

The largest difference was observed in radially symmetric aberration spherical aberration (Z₁₂) and asymmetric aberration secondary astigmatism (Z₁₃), with magnitudes of 0.16 µm and 0.18 µm, respectively. The RMS differences were 0.28, 0.51, and 0.50 µm for IOP elevations of 5, 15, and 25 mm Hg, respectively. For cornea 2 (Fig. 9b), the predicted wavefront map shows a much larger difference, mainly due to the deviation in radially symmetric aberration defocus (Z₄), with a magnitude of 0.96 µm. The smallest and largest differences were observed in astigmatism (Z₅) and defocus (Z₄), with magnitudes of 0.13 and 0.96 µm, respectively. For the three IOP elevations (5, 15, and 25 mm Hg), the corresponding RMS difference for the wavefronts were 0.48, 0.78, and 1.14 µm, respectively. 

Discussion

An individualized biomechanical model of the cornea could improve predictions of optical behaviors; however, characterizing collagen fibril dispersion, particularly in vivo, has been extremely challenging. The present study proposed a novel method of characterizing the spatial distribution of collagen fibril dispersion of individual corneas using IOP-induced changes in corneal aberrations. Preliminary numerical and experimental studies demonstrated the feasibility of the method and established the potential of developing an individualized FE model of the cornea used for predicting optical and mechanical behaviors of the cornea. 

Several techniques, such as x-ray diffraction^15,16 and second harmonic generation imaging microscopy,^42,43 have been used to investigate the collagen fibrils structure of the cornea and their distribution. However, challenges still exist in applying this information to FE modeling of the cornea for biomechanics, where a continuous representation of the collagen fibril distribution over the entire cornea is needed.²² Our preliminary studies suggest that there is an approximately linear correlation between the collagen fibril dispersion and the IOP-induced changes in corneal aberrations, as illustrated in Figure 10. Therefore, we have proposed mathematical equations to correlate the biomechanical properties (geometry and material properties) of the cornea and the IOP-induced changes in its wavefront aberrations. Because IOP-induced changes in wavefront aberrations due to measurable biomechanical parameters such as corneal geometry^35,36 and matrix stiffness,^4,5 can be quantified, the changes due to collagen fibril dispersion can be obtained by subtracting the corresponding contribution from corneal geometry and matrix stiffness from the total changes. A straightforward measurement of corneal wavefront aberration, which is also a continuous representation of the entire corneal surface, may then be used to characterize collagen fibril distribution. 

Figure 10.

 

Examples of the linear correlation between the value of kappa and the changes in corneal optical aberrations induced by an IOP elevation of 15 mm Hg IOP. (a) Change in astigmatism. (b) Change in spherical aberration.

The method of characterizing the collagen fibril dispersion parameter κ(ρ, θ) was first tested through numerical simulations, with the main features estimated well as shown in Figure 5. The observed consistent underestimation of the magnitudes of κ(ρ, θ) may be due to two possible reasons. First, the extrapolation process of the wavefront data from the central 6 mm to a full corneal size may induce errors. Because it is difficult to predict the real wavefront at the periphery, we used a linear extrapolation algorithm and assumed the value of κ at the edge to be 0 based on the known alignment of the fibrils at the limbus.^3,15 Some features at the transition area may not be captured with this method; however, the results were quite reasonable within the 6-mm diameter size. Second, we did not consider the contribution of the interaction effects between important biomechanical parameters (geometry, matrix stiffness, and fibril dispersion) on the IOP-induced changes in corneal wavefront aberrations. Our previous sensitivity study¹⁹ showed that the interaction between the matrix stiffness and fibril dispersion had an impact on the IOP-induced changes in corneal optical aberrations, but the effects of interactions among geometry, matrix stiffness, and fibril dispersion have not yet been clearly quantified. Neglecting these interaction effects may also affect estimation of collagen fibril dispersion. Because collagen fibrils dominate corneal biomechanical behavior at higher pressure,⁴⁴ the prediction error due to the underestimation of κ(ρ, θ) increased with IOP for each pattern. However, the overall predictions of IOP-induced changes in wavefront aberrations for the central 6-mm cornea from the three corneal models were still very close to the simulated experimental data, as shown in the comparison of the wavefront maps in Figure 6. 

For the ex vivo experiments on monkey corneas, the estimated distribution of fibril dispersion (Fig. 8) showed a much larger difference between two different animals (cornea 1 or 2 vs. cornea 3) compared to a paired cornea from the same animal (cornea 1 vs. cornea 2), which may suggest a much larger inter-specimen variability than inter-eye variability. In addition, the matrix stiffness of cornea 3 (10 years old) was approximately two times higher than the values for corneas 1 and 2 (156 days old) (Fig. 7) and showed a smaller change in wavefront aberrations with the same IOP elevation (Fig. 9), thus supporting the finding that corneal stiffness increases with age.^31,45 For the prediction of the IOP-induced changes in wavefront aberrations, larger errors were observed in the FE models developed from ex vivo experiments compared to the theoretical simulations. For the three IOP elevations, the average RMS differences were 0.57, 0.80, and 0.43 µm for corneas 1, 2, and 3, respectively. As shown in Figure 9, the wavefront maps representing the IOP-induced changes predicted from two of the corneal models (corneas 1 and 3 in Figs. 9a and 9c) maintained the overall features of the experimental data, whereas the other one (cornea 2 in Fig. 9b) showed a much larger difference. The characterization of fibril dispersion was less accurate for cornea 2 compared to its paired cornea 1 from the same animal. 

A few factors could lead to errors in the ex vivo experiment. First, it has been known that tissue hydration could significantly affect the characterization of corneal biomechanical properties.^46–48 We tried to carefully control tissue hydration in all stages of the experiment. Corneas were completely submerged in Optisol-GS during the inflation test. After the top chamber was removed for geometry measurement, Optisol-GS was added on top of each cornea every 2 to 3 minutes to keep the anterior surface hydrated. In the last step of the shear test, the thickness of the specimen was compressed to match its initial thickness measured with OCT to avoid the impact of hydration state on the measurement accuracy. However, corneas could still swell to a certain degree due to the long experimental procedure. Although corneas 1 and 2 were tested sequentially in all of the three measurements, cornea 2 was tested 2 to 3 hours later than cornea 1, including the setup time and data collection time required for the experiments, which could have led to a different hydration status and affected the accuracy in its measurement results. For example, wavefront measurement of the cornea is very time sensitive due to reductions in quality of aberration measurements even though the corneas were submerged in Optisol-GS. 

Second, for all three corneas, the main difference was observed in defocus (Z₄), especially for cornea 2. Because defocus is a radially symmetric aberration, this deviation may indicate an offset of matrix stiffness C₁₀. The tested corneal explants were cut into cylinders as required for the shear test, thus effects on integrity due to the cutting may affect the measured material properties.⁴ Corneal tissue must also be highly and evenly cellularized in order to track the cell movement for depth-dependent strain calculation.⁴⁹ Because the shear test was the last step in our experimental sequence, decreased quality of the cornea reduced the number of trackable cells, which may have affected the quality of the data. In addition, external factors such as air bubbles, fluid flow, or specimen slippage during the measurement procedure may also induce certain errors in the measured value of C₁₀. For the same reasons mentioned above in terms of the poor quality of the cornea during the experimental process, we failed to collect a complete dataset for the other three corneas: two corneas with poor transparency for wavefront measurement and one with excessive hydration for shear modulus measurement. 

Characterizing the distribution of collagen fibril dispersion using optical aberrations provides potential for future development of an in vivo individualized biomechanical model of the cornea. Inverse analysis combined with air-puff deformation measurements has been explored to develop an in vivo corneal model.^50,51 However, the derivation of corneal material properties from such iterative procedures based on just apical displacement or one meridional profile of the cornea may not be unique,⁵² which could generate different optical results. Measurements of corneal wavefront aberrations and IOP elevation are the two key factors for characterizing fibril dispersion with the new method proposed in this study, and they are achievable in vivo. For wavefront aberrations of the cornea, there are techniques available for in vivo measurement, such as wavefront sensing^53,54 and corneal topography measurement.^55,56 For in vivo IOP elevation, we have proposed an innovative, safe, non-contact method to temporarily elevate IOP in vivo using an inversion table, and wavefront aberrations may be measured simultaneously.⁵⁷ Collagen fibrils begin to dominate corneal biomechanical behavior at a higher IOP; thus, the maximum IOP elevation was used for more accurate estimation of collagen fibrils in both the numerical simulations (40 mm Hg) and ex vivo experiments (25 mm Hg) in our study. Whether a more modest degree of IOP elevation in vivo is sufficient for accurate estimation of collagen fibril dispersion could depend on the interaction of the matrix stiffness and collagen fibrils of individual corneas and requires further investigation. 

Moreover, several newer technologies, such as optical coherence elastography⁵⁸ and Brillouin optical microscopy,⁵⁹ the results of which might be correlated with the depth-dependent matrix stiffness of the cornea, offer the potential to characterize the depth-dependent matrix stiffness of individual corneas. The development of an in vivo individualized corneal biomechanical model with customized distribution of fibril dispersion using our method might be possible in the future by combining these available techniques, which could be critical for several clinical applications, such as early diagnosis of keratoconus,^60–62 risk evaluation for post-refractive ectasia,^63,64 and the planning of refractive surgery.^65,66 Previous x-ray diffraction studies have reported differences in collagen orientations and distribution between healthy and keratoconus corneas.²⁰ These findings suggest that displacement of the axes of the collagen fibrils and distortion of the orthogonal matrix may be related to the formation of the cone shape of a keratoconus cornea. The ability to characterize the distribution of collagen fibril dispersion could be helpful for early diagnosis and possibly management of keratoconus. For refractive surgery, the loss of the structural integrity of the collagen fibrils associated with excess tissue removal could also be a possible cause for post-refractive ectasia.⁶⁷ In addition, individual variations in the distribution of collagen fibril dispersion could be one of the main biomechanical factors for the large inter-patient variability in the optical outcomes. A better understanding of the distribution of collagen fibril dispersion, especially with further incorporation into an individualized biomechanical model of the cornea, could be helpful for predicting more accurate surgical outcomes, including preventing post-refractive ectasia, in addition to reducing the post-refractive higher order aberrations. 

Several aspects of the current study could be improved in future work. First, we found that small errors in corneal geometry measurement may have a significant impact on the estimation of collagen fibril dispersion. We used the OCT system to measure the posterior corneal surface of ex vivo corneas due to the limitations of our corneal topographer. Due to the intrinsic limitations of OCT, the image signal on the peripheral edge of the corneal surface, especially the posterior surface, is weaker, which may introduce errors into the geometrical model. In vivo topography measurement of both anterior and posterior corneal surfaces could improve the accuracy of geometrical model^68,69 and the following material characterization process. Second, the ex vivo wavefront measurement requires very high corneal quality, so we used rhesus monkey corneas instead of human corneas for better control. It has been known that significant variations exist in the structure of corneal collagen fibrils in different species.^10–12 Although x-ray studies have suggested that structurally the human cornea may mirror that of the adult marmoset more closely than that of any other animal species,⁷⁰ there are still clear differences in their collagen fibril arrangement.^10,11 Rhesus monkey corneas were the only available primate resources we had for the current study, and we assumed that their collagen fibril orientation was the same as the human cornea in our modeling because the exact fibril orientation of the rhesus monkey cornea is still unclear. This is one limitation in the current study that requires further investigation. X-ray diffraction studies on rhesus monkey corneas or further validation of the method on human corneas could be helpful in the future. 

In conclusion, a new method of characterizing the distribution of collagen fibril dispersion of individual corneas was proposed and demonstrated through both numerical simulation and ex vivo experiments. Our findings suggest that spatial variations in the collagen fibril dispersion of the cornea could be estimated from the changes in corneal optical aberrations induced by IOP elevation, and an individualized corneal model with customized distribution of fibril dispersion could improve predicting the optical behaviors of a specific cornea. Although the small sample size does not allow us to reach a definitive conclusion, the results from this exploratory study warrant further investigation. The proposed method could advance development of an in vivo individualized biomechanical model of the cornea, as well as investigations into the impact of corneal biomechanics on the visual performance of patients for clinical applications. 

Acknowledgments

The authors thank Earl Smith III, OD, PhD, and Baskar Arumugam, PhD, University of Houston, for providing the corneas of rhesus monkeys. 

Supported by grants from National Institutes of Health/National Eye Institute (EY014999) and the National Science Foundation (CMMI-1100632), and by Research to Prevent Blindness. 

Disclosure: M. Xu, None; M.A. Ramirez-Garcia, None; H. Narang, None; M.R. Buckley, None; A.L. Lerner, None; G. Yoon, None