Constitutive Modeling of Corneal Tissue: Influence of Three-Dimensional Collagen Fiber Microstructure

The cornea protects the inner contents of the eye against external insults, provides about two-thirds of eye's refractive power, and transmits nearly 90% of the incident light onto the lens [1,2]. The proper optical function of the cornea depends on its ability to maintain its precise shape under physiological loading conditions. The corneal extracellular matrix, stroma, constitutes almost the entire corneal thickness and serves as the key component in providing the mechanical strength necessary to resist external and internal forces.

The microstructure of the stroma resembles a lattice structure, where collagen fibrils are embedded in thin parallel-to-the-surface lamellae [3,4]. The X-ray scattering methods gave detailed information about the preferred collagen fibril orientation in the corneal stroma [5–9]. In particular, it was found that although collagen fibrils are oriented along with two preferred directions, i.e., nasal-temporal (N-T) and inferior–superior (I–S) within the central region of the human cornea, they tend to be aligned circumferentially in the limbus region. Unlike the human cornea, collagen fibrils in bigger mammals such as bovine are found to be aligned mainly in the I–S direction [10,11]. The degree of in-plane dispersion varies in depth, i.e., although collagen fibrils are more aligned along with N–T and I–S directions within the posterior thirds, they are more isotropically oriented within the anterior thirds [12]. The images obtained from the X-ray scattering technique could not fully characterize the 3D dispersion of collagen fibrils through the corneal thickness. Nevertheless, the second harmonic generation (SHG) images, about a decades ago, provided a reconstruction of 3D collagen fibril orientation [13,14]. These images showed that collagen fibrils are highly interwoven in the anterior region but are parallel to each other in the posterior region.

The early works utilized simple linear-elastic or hyperelastic models for representing the corneal constitutive response [15,16]. Later, linear transverse anisotropic models were used to account for the anisotropic response [17,18]. Hyperelastic models considering both isotropic and anisotropic contributions were also used. In these models, the dispersion of collagen fibrils was not considered in early works [19,20], but was added later [11,21–26]. These recent models could be categorized into two groups: the angular integration (AI) models and the generalized structure tensor (GST) models.

The AI-based models have a straightforward formulation, where the free energy corresponding to the continuous collagen fibril distribution is obtained by performing the direct angular integration of an infinitesimal fraction of fibers in a given direction. The statistical description of the collagen fibril distribution could be represented by either distribution probability density function (PDF) or direct extrapolation of the X-ray scattering data. The AI models with different forms of PDFs have been applied to various soft tissues [27–33]. The AI models provide relatively good representations of the mechanical response of biological tissues. However, their main disadvantage is that the numerical implementation of their required direct angular integration scheme is complicated and time-consuming.

On the other hand, the GST models are relatively faster. They use the generalized structure tensor with a dispersion parameter to quantify the collagen fibril dispersion [34]. Once the dispersion parameter is specified, the stretching of collagen fibrils at any given macroscopic deformation is known and the required angular integration can be evaluated. However, this model can be used with the limited number of PDFs for collagen fibril orientation because the derivation of analytical relations between PDFs and dispersion parameters is not trivial [34,35].

The collagen fibril distribution in SHG images suggests that both in-plane and out-of-plane dispersions are essential. In this work, we use a GST model that takes into account both in-plane and out-of-plane collagen fibril distribution throughout the cornea. The in-plane distribution is approximated by fitting the normal distribution function in polar coordinates to the X-ray scattering data [9]. The out-of-plane distribution of collagen fibrils at a given thickness level has been represented by fitting Gaussian curves to the cutoff angle histogram obtained from the SHG images [14]. We numerically implement the proposed GST model in a commercial finite element software abaqus/standard [36] by writing a user-defined material subroutine (UMAT). The model performance is studied by simulating the results of inflation tests [37] using six different collagen fibril distribution of transversely isotropic (T.I.), isotropic (I.), perfect alignment (P.A.), planar dispersion (P.D.), planar isotropic (P.I.), and full-thickness variation (F. V.). A parametric study is also performed to determine the effects of collagen fibril interweaving on stress profiles across the corneal thickness. Finally, it is shown that the proposed model has similar functionality as the available AI model from the literature [26], yet cheaper in terms of the computational expenses.

The remainder of this paper is organized as follows. In Sec. 2, we review the continuum mechanical framework and present the main constitutive equations. The governing equation is briefly summarized in Sec. 3. The numerical results are shown in Sec. 4. Finally, we finish in Sec. 5 with some concluding remarks. In Appendix  A, we present the details of our code verification.

This section covers the large deformation kinematics required for describing the hyperelastic anisotropic behavior of the corneal stroma. A similar framework has been previously applied to soft materials [38–40].

Note that $κip$ and $κop$ in the above expression are in-plane and out-of-plane dispersion parameters whose characteristics will be discussed in the following.

where $κip∈[0,1]$ and $κop∈[0,1/2]$ are dispersion parameters, and $I1(a)$ is the modified Bessel function of the first kind of order 1. In Fig. 2, we project the total PDF in Eq. (7) onto the surface of a unit sphere with different combinations of in-plane and out-of-plane dispersion parameters; here one family of fibers with orientation a₀ that is aligned with the unit vector $N=[1,0,0]⊤$ is considered. The out-of-plane normal is set to be $an=[0,0,1]⊤$⁠. As $a→0$ and $b→0$⁠, the collagen fibrils are evenly distributed. Conversely, as $a→∞$ and $b→∞$⁠, the collagen fibrils are perfectly aligned long with the mean orientation. The collagen fibrils are isotropically distributed within $x1−x2$ plane as $a→0$ and $b→∞$⁠, and are isotropically distributed within $x1−x3$ plane as $a→∞$ and $b→0$⁠. Accordingly, the generalized structure tensor H for one family of fibers could be simplified into five special cases:

• •

  Perfect alignment—$H=a0⊗a0$⁠;

• •

  Isotropic dispersion—$H=(1/3)1$⁠;

• •

  Transversely isotropic—$H=κ1+(1−3κ)a0⊗a0$ when $κ=1−2κop$⁠;

• •

  Planar dispersion—$H=kI+(1−2κ)a0⊗a0$ when I is the 2D identity and k is the dispersion parameter in the plane;

• •

  Planar isotropic—$H=(1/2)I$⁠.

where G₀ denotes the ground state shear modulus and K denotes the bulk modulus.

It is worth noting that the collagen fibril free energy does not have any volumetric contribution to the total free energy. Furthermore, collagen fibrils are not able to withstand any compression, so if $I4¯≤1$ and $I6¯≤1$⁠, the free energy $ψRfi$ is completely omitted in Eq. (10).

are the stress contributions from the underlying matrix and collagen fibrils, respectively.

where n is the out-normal to $∂Bt$⁠, and $[[•]]$ is the jump operator, defined as the difference between the quantity inside and outside the domain.

The proposed model is implemented numerically in abaqus/standard [36] by writing a UMAT, and its verification is found in Appendix  A. In this section, we investigate the capabilities of the proposed model by simulating the standard inflation test.

We use the previous inflation experimental results of Anderson et al. [37]. In these experiments, porcine corneal samples, with the narrow ring of surrounding scleral tissue, were mounted such that the portion that connects the limbus and the sclera was fully fixed. An internal pressure with a maximum value of 100 mmHg was gradually (quasi-static condition) applied to the samples' posterior surface. Meanwhile, the apical displacement was continuously monitored by a CCD laser displacement sensor and plotted against the pressure.

Here, G is the maximum vertical height at r = 0, both $Rx1$ and $Rx2$ are the maximum curvatures of the principal meridians along x₁ and x₂ directions, respectively. $Θx1$ is the direction of the steepest principal meridian, both $Qx1$ and $Qx2$ are the asphericity parameters in the directions $Θx1$ and $Θx1+π/2$⁠, respectively.

We use referential unit vectors $a04$ and $a06$ to represent the two mean orientations of collagen fibrils (see Fig. 3(a)). In the central region, two families of the collagen fibrils are running from N–T (dashed-line) and I–S (dotted-line) directions in a 3D curved fashion. In the limbus region, one family of collagen fibrils (dashed-line) is running circumferentially, and another (dotted-line) is pointing outward from the center. Additionally, the out-of-plane direction is denoted as the unit vector a_n (solid line).

where $κopmin=1/3$ and $κopmax=1/2$ are minimum and maximum value, respectively, and the constant γ_d controls the nonlinearity of the function (see Fig. 4(b)).

The geometry is discretized into U3D8 elements with six elements spanning the thickness, and only a quarter of the entire mesh is presented for clarity (see Fig. 3(b)). For boundary conditions, we fully fix the surface linking the limbus and sclera and applied an internal pressure of P = 100 mmHg to the posterior surface.

where $Xk$ and $Xkdef$ denote the zero-pressure and deformed coordinates at kth step. Meanwhile, the mean collagen fibril orientations are consistently mapped back to the zero-pressure configuration. Here, the iteration is terminated based on the global error $e=||Xkdef−XPhysio||∞$⁠. The parameters of the biconic Eq. (19) used for the anterior surface under P = 16 mmHg are obtained from previous studies [24], i.e., $Rx1=7.71$ mm, $Rx2=7.87$ mm, $Θx1=0.51π, Qx1=Qx2=−0.41$ and G =2.52 mm. The parameters used for the posterior surface under P = 16 mmHg are $Rx1=6.36$ mm, $Rx2=6.69$ mm, $Θx1=0.51π,Qx1=Qx2=−0.52$ and G = 1.89 mm. We plot the physiological coordinates as stars, while the deformed and zero-pressure coordinates at each iteration as triangles and squares (see Fig. 6(a)). It is observed that the global error is minimized quickly within about ten iterations (see Fig. 6(b)).

We consider six different collagen fibril distributions, i.e., T.I., P.A., I., P.D., P.I., and F.V., in our simulation, which are based on material parameters given in Table 1. The parameters are selected, such that the apical rise–pressure curves of both F.V. and T.I. fall onto the experimental data as close as possible (see Fig. 7(a)). The other four cases are simulated using their respective dispersion parameters while keeping mechanical parameters unchanged.

We compare the contours of Von-Mises stress among all cases under internal pressure of P = 100 mmHg (see Fig. 8). For the cases of T.I., P.A., and P.D., the contours share a similar pattern—a cross mark in the central region—indicating that the collagen fibrils along N-T and I-S directions are under tension. In the case of F.V., the Von-Mises stress is much lower at the anterior surface. It is because collagen fibrils near the posterior surface are almost perfectly aligned, making them exhibit an earlier stretch-locking than the collagen fibrils at the anterior surface. In the cases of I. and P.I., no significant stretching of the collagen fibrils is observed.

Echoing the main focus of the current study—modeling structural variation in collagen fibrils across the corneal thickness—we plot in Fig. 9 the side view of the same Von-Mises stress contour. The F.V. case could predict a reasonable stress profile across the thickness that is in line with the observed collagen fibril distribution in SHG images. On the other hand, the stress is found to be more concentrated at the anterior surface in the cases of T.I., P.A., and P.D. There is no apparent stress gradient across the thickness in the cases of I. and P.I.

In Fig. 7(a), the simulated apical rise–pressure curves are plotted as lines, while the experimental data obtained from Anderson et al. [37] are plotted as circles. The cases of T.I. and F.V. can capture the experimental results quite well. The case of P.A. exhibits the earliest stretch-locking behavior, while the case of I. shows no fiber engagement under the same boundary conditions. Interestingly, for planar cases of P.D. and P.I., the collagen fibrils exhibit a relatively earlier stretch-locking behavior caused by the narrower dispersion space.

Here, we investigate the influence of decay rate constant γ_d on the mechanical response of corneal stroma under inflation. Figure 7(b) compares the apical rise–pressure curves under different values of γ_d. It is seen that the “stretch locking” behavior occurs earlier as γ_d increases. It is because a large fraction of collagen fibrils with the narrower out-of-plane dispersion is engaged in the deformation. The case with $γd=2.5$ has a good agreement with the experimental results. More interestingly, as Fig. 10 shows, the spots of high-stress concentration move from the posterior to the anterior side as γ_d increases. This result could be useful in terms of predicting the location of the highest stress across the thickness with the known collagen fibril architecture.

where $dΩ=sin ΦdΦdΘ$ denotes the differential unit sphere, $a0$ denotes the referential collagen fibril orientation, $ψRm$ denotes the matrix free energy, and w denotes the collagen fibril strain energy, which is a function of fiber stretch $λf$⁠. The integral in Eq. (26) is normalized by $m=∫Ωρ(a0)dΩ$⁠. We implement the AI model numerically in abaqus/standard [36] by writing UMAT. The simulated apical rise–pressure curve using the AI model is shown in Fig. 7(a) as the dotted black line, which shows a good agreement with the experimental data. The CPU elapsed time of simulations under the different number of unconstrained degrees-of-freedom $ndof$ is recorded. It is found that the CPU elapsed time ratio between the two models $tAI/tGST$ is proportional to the number of unconstrained degrees-of-freedom $ndof$ (see Fig. 11). Since the double angular integration in Eq. (26) is evaluated by the Gaussian quadrature scheme, the number of Gauss points has a significant impact on the numerical expensiveness, where the time ratio is nearly 100 at $ndof=140$ by using 64 Gauss points. Overall, the proposed model has almost the same feature as the AI one, but it is cheaper in terms of the simulation cost.

This work develops a continuum mechanics model considering collagen fibril out-of-plane dispersion across the corneal thickness. In particular, the function that links the out-of-plane dispersion parameter to the corneal thickness serves as one of the important contributions of the current work. The proposed model is numerically implemented, and its capabilities are demonstrated by performing numerical simulations of inflation experiments using six different collagen fibril orientations, i.e., transversely isotropic, isotropic, perfect alignment, planar dispersion, planar isotropic, and full-thickness variation.

The results show that the proposed model can replicate the experimental pressure displacement curves very well. It also predicts a reasonable stress profile across the corneal thickness. A parametric study on the decay rate constant indicates that the stress profile across the thickness is sensitive to the collagen fibril out-of-plane dispersion. From the perspective of computational expenses, compared to a recent AI model [26], the performance of the proposed model stands out because while it requires much less computational power, it has almost the same functionality.

Looking toward the future, more work yet remains. For example, the model could be strengthened by incorporating more detailed collagen fibril structural information—the interactions between collagen fibril layers. The model could also be extended by adding dissipative mechanisms such as viscoelasticity and considering the corneal gel-like behavior—poroelasticity and fluid migration [45,46].

• National Science Foundation, the division of Civil, Mechanical, and Manufacturing Innovation (Grant 1635290; Funder ID: 10.13039/100000001).

where P is a constitutively indeterminate pressure that is introduced to satisfy the incompressibility constraint. On the numerical side, a single element (U3D8) in abaqus/standard [36] is prescribed with the same deformation. We also take $K=103G0$ to secure a nearly incompressible condition in the simulations.

where two constants are $α=(2B+3A−1)axay$ and $β=A+Bay2+(1−3A−B)ax2$⁠, respectively. The stress is normalized by the matrix shear modulus G₀, and three different combinations of dispersion parameters are considered. The excellent agreement between the numerical and analytical solutions suggests that our numerical implementation is fully verified.