Effect of a Gradient Distribution of Cross-Links on the Deformation Behaviors of Corneal Stroma: Theoretical Model and Finite Element Simulation

1 Introduction

Keratoconus is a degenerative and progressive corneal disease characterized by posterior ectasia and corneal thinning. With the development of advanced detection technology, the diagnostic rate of keratoconus has been noticed to keep increasing in recent years (Mas Tur et al., 2017; Buzzonetti et al., 2020; Rahmati et al., 2021; Tian et al., 2021). From a biomechanical point of view, the possible reasons resulting in keratoconus may originate from the abnormal structure and arrangement of collagenous fibers, which weaken the ability of corneas to resist deformation, and make the central or paracentral corneas expand, get thinner and gradually bulge outward into a cone shape. In the early stage, patients with keratoconus could correct the vision problems by wearing glasses or contact lenses. However, these treatments can not fundamentally prevent the development of keratoconus, but are only limited to help improve the vision and refractive degree (Moshirfar et al., 2017; Lin et al., 2018; Vinciguerra et al., 2021). As an alternative treatment, corneal collagen cross-linking (CXL) has been experimentally proven to be able to modify the biomechanical properties of keratoconus by improving the strength and stability of corneal stroma (Wollensak et al., 2003a; Wollensak et al., 2003b). Therefore, it has now been taken as an effective method to clinically treat keratoconus by CXL (Ambekar et al., 2011; Hovakimyan et al., 2012; Raghunathan et al., 2017).

By taking riboflavin as a photosensitizer and irradiation with ultraviolet A (UVA), the stabilization of ectasia and stiffness of corneal stroma in the anterior layer can be increased through the formation of intra- and inter-fibrillar covalent bonds by photosensitized oxidation (Wollensak and Iomdina, 2009; Sorkin and Varssano, 2014; Blackburn et al., 2019; McKay et al., 2019; Rubinfeld et al., 2019; Arrizabalaga and Nollert, 2020). The effectiveness of standard CXL has been clinically and experimentally proved over the last years (Mazzotta et al., 2019; Rubinfeld et al., 2019; Herber et al., 2020), and it is realized that the modification of the UVA irradiation mechanism (Wollensak et al., 2003a) and change of the transmittal mode of riboflavin may help improve the operation efficiency and patient experience (Kymionis et al., 2017; Kobashi and Tsubota, 2020). More importantly, based on clinical observations, the key parameter determining the CXL effect on the corneal stiffening behavior is informed to be the total irradiation energy or say the irradiation dose (Kymionis et al., 2017; Woo et al., 2017; Bao et al., 2018; Marafon et al., 2020; Turhan et al., 2020). Although the variation of the biomechanical responses induced by CXL has been well acknowledged, the existing issue is how to quantitatively characterize the increment of corneal stiffness and strength as a function of the irradiation dose. In order to effectively address this issue and comprehend the underlying deformation mechanisms, increasing interests have been focused on the study of the CXL effect on the biomechanical properties of corneal stroma from the aspects of experimental observations, theoretical analyses and numerical simulations.

In the field of experimental investigation, several different measuring methods, including uniaxial tensile tests, nano-indentation and noncontact Brillouin measurements, have been considered to address the variation of the biomechanical properties of corneal materials after CXL (Kohlhaas et al., 2006; Scarcelli et al., 2012; Scarcelli et al., 2013; Seifert et al., 2014; Dias et al., 2015; Liu et al., 2020b). For instance, Kohlhaas et al. (2006) applied uniaxial tensile tests on two flaps of the porcine and human corneas after the treatment by CXL, and noticed an obvious depth-dependent stiffening behavior that the anterior-treated flap with a thickness of 200 μm was much stiffer than the posterior-treated flap. Similarly, Liu et al. (2020a) compared the tensile experimental data of four overlapped strips after CXL, and found that the elastic modulus decreased non-linearly as a function of the corneal thickness. In addition, by applying an atomic force microscopy-derived indentation method, Seifert et al. (2014) obtained a depth-dependent distribution of the Young’s Modulus for cross-linked porcine corneas, and can be explained by the gradient stiffening effect related to the standard CXL. More interestingly, Scarcelli et al. (2013) indicated the difference of corneal elastic modulus and stiffness for the anterior, middle and posterior stroma through the Brillouin maps, which effectively validated that the noncontact Brillouin microscopy was able to quantify the CXL effect on the corneal mechanical properties. Informed by these experimental works, one may realize that the CXL treatment can effectively facilitate the improvement of corneal stiffness and strength. Moreover, this strengthening behavior tends to decrease with increasing corneal thickness, and ultimately results in an obvious gradient distribution of materials properties for cross-linked corneas.

To further unveil the fundamental deformation mechanisms corresponding to the improvement of the mechanical property of corneal materials after CXL, theoretical analysis within the framework of strain energy density function theory has been mainly focused on in recent years (Scarcelli et al., 2012; Scarcelli et al., 2013; Saez et al., 2014; Tian et al., 2016; Pandolfi et al., 2019; Cornaggia et al., 2020; Holzapfel and Ogden, 2020). A preliminary attempt was considered in the work of Liu et al. (2020b), who applied a collagen fibril crimping model to address the uniaxial tensile data of cross-linked corneal stroma. Different sets of model parameters were calibrated by comparing the theoretical results with the experimental data after UVA irradiation. In order to fundamentally incorporate the CXL effect into the constitutive equations, a 3D strain energy function accounting for the relative orientation and stiffness of cross-links, fibers and their coupling was deduced by Holzapfel and Ogden (2020), which can effectively characterize the increment of materials stiffness with increasing cross-link density. To further quantitatively characterize the stiffening behavior as a function of the UVA irradiation dose, a novel microstructure-based model was recently proposed that the strain energy function was explicitly specified in terms of the density and orientation of cross-links, and was able to predict the uniaxial tensile data at various irradiation doses. Moreover, the increment of corneal strength was observed to be determined by not only the cross-link density but also the cross-link orientation with respect to the collagen fibers (Xiao et al., 2022). However, though great effects have been devoted in the field of theoretical analysis, one of the characteristic features as informed from the experimental measurements, i.e., the cross-links formed in corneal stroma are gradiently distributed along the thickness direction, has not been well addressed in the constitutive laws so far.

Furthermore, in comparison with the limited application of theoretical models for uniaxial tensile tests, a more rational characterization of the biomechanical behaviors of cross-linked corneas can be achieved through the utilization of finite element simulations (Pandolfi et al., 2019; Cornaggia et al., 2020; Wang and Chester, 2021). For example, Pandolfi et al. (2019) applied a simplified microstructure model that mainly consisted of a trusswork with two main sets of fibers and related cross-links to simulate the expansion deformation of both healthy and keratoconus cornea. Simulated results indicated that the weakening of cross-links can lead to noticeable deformation of corneal stroma, and finally dominate the localized thinning of keratoconus corneas. Wang and Chester (2021) combined a multi-physics model and finite element formulation to analyze the underlying mechanisms induced by cross-links from the aspects of both nano-indentation and inflation tests, which may offer an effective way to help make predictions of a CXL therapy based on the real corneal geometry. Very recently, Cornaggia et al. (2020) simply divided the cross-linked corneas into two regions with/without CXL, and assumed piecewise elastic constants for the two adjacent regions along the corneal thickness direction. By taking the finite element model for the further simulation of inflation tests, it was realized that the increase of the equivalent material strength was mainly ascribed to the localized increment of material stiffness in the anterior layers. For all the finite element simulation works as mentioned above, one may note that the UVA irradiation-induced cross-links are either homogenously or piecewise-linearly distributed in the corneal stroma, which is not consistent with the experimental observations. Therefore, further developments of the finite element simulations should be concentrated on the consideration of the inhomogeneous distribution of cross-links along the corneal thickness direction.

In this work, we intend to combine both theoretical model and finite element simulation to address the influence of gradiently distributed cross-links on the fundamental deformation mechanisms of corneal materials. Within the framework of nonlinear continuum mechanics, a hyperelastic mechanistic model is proposed that takes into account the stiffening contribution of both depth-dependent and irradiation dose-related cross-links. Moreover, the developed constitutive laws are incorporated into the subroutine UMAT of ABAQUS to simulate the inflation tests of porcine cornea treated by CXL. This paper is organized as follows: in Section 2, the constitutive equations for corneal stroma with gradient distribution of cross-links are introduced in details, which contain the expressions for the strain energy density function, Cauchy stress tensor and elasticity tensor of each stiffening component, the setup of the geometrical model for finite element simulation and calibration of model parameters. In Section 3, the developed simulation framework is applied to address the inflation tests of porcine cornea at different irradiation doses. In Section 4, corresponding deformation mechanisms are systematically discussed. Finally, we close with the main conclusions in Section 5.

2 Materials and Methods

2.1 Microstructure Morphology of Corneal Materials After CXL

In order to deduce the constitutive laws addressing the fundamental physical mechanisms, the microstructure information of corneal stroma after CXL is briefly reviewed here. Experimental characterization has indicated that corneas mainly consist of five individual layers, and the stroma layer takes nearly 90% of the corneal thickness. Therefore, the hyperelastic biomechanical behaviors of corneal materials are generally considered to be dominated by the mechanical properties of corneal stroma. At micro-scale, corneal stroma contains 250–400 highly ordered and stacked fibrous lamellae, which are composed of collagen bundles with the same diameter and orientation. Further observation with the wide angle X-ray scattering technique (Aghamohammadzadeh et al., 2004; Meek et al., 2005; Meek and Knupp, 2015) revealed that almost 66% of the collagenous fibers are distributed along the nasal-temporal (N-T) and superior-inferior (S-I) direction, and the rest fibers are randomly distributed in the matrix. Considering the expansion deformation of the cornea under normal intra-ocular pressure (IOP), the angle α between the localized fibers and horizontal reference surface will vary along both the N-T and S-I direction, Figure 1A for an illustration. Therefore, the direction vectors of the two localized fibers can be, respectively, characterized by ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$ with ${\mathbf{a}}_1=\left[-\cos \alpha 0\sin \alpha \right]$ and ${\mathbf{a}}_2\mathrm{=}\left[0\cos \alpha \sin \alpha \right]$. Correspondingly, the projection vectors of ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$ on the projection plane (i.e., the x-o-y plane) give as ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{a}}_1}=\left[-\cos \alpha 00\right]$ and ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{a}}_2}=\left[0\cos \alpha 0\right]$.

FIGURE 1. (A) Side view and (B) anterior view of the corneal stroma with two families of fibers and corresponding cross-links. α is the angle between the localized fibers and horizontal surface. Cross-links distribute symmetrically with respect to the collagen fibers with an angle of β on the projection plane (i.e., the x-o-y plane). ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{a}}_1}$, ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{a}}_2}$, ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{m}}^{\pm }}$ and ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{n}}^{\pm }}$ are the vectors for fibers and cross-links on the projection plane, respectively.

Besides the statistically distributed collagen fibers as well as the isotropic elastin and proteoglycan matrix, additional covalent cross-links can be formed between adjacent collagen fibers during the process of collagen CXL due to the synergistic effect of ultraviolet irradiation and riboflavin catalysis (Lombardo et al., 2015). The CXL technique can artificially lead to the increase of cross-linking degree, and thus effectively facilitate the strength increment and mechanical stability of corneal stroma (Sorkin and Varssano, 2014; Nohava et al., 2018; Damgaard et al., 2019; Zhang D. et al., 2021; Zhang H. et al., 2021). Given that the distribution of cross-links can not be visualized directly by staining methods or other microscopic techniques so far, it is generally considered that the cross-links tend to distribute symmetrically with respect to the collagen fibers with an angle of β on the projection plane (Holzapfel and Ogden, 2020), as illustrated in Figure 1B. In this way, the direction vectors of the cross-links corresponding to the two main families of fibers can be described by ${\mathbf{m}}^{\pm }$ and ${\mathbf{n}}^{\pm }$ with ${\mathbf{m}}^{\pm }=\left[-\cos \alpha \cos \beta \pm \cos \alpha \sin \beta \sin \alpha \right]$ and ${\mathbf{n}}^{\pm }=\left[\pm \cos \alpha \sin \beta \cos \alpha \cos \beta \sin \alpha \right]$, respectively. Similarly, the projection vectors of ${\mathbf{m}}^{\pm }$ and ${\mathbf{n}}^{\pm }$ on the x-o-y plane yield as ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{m}}^{\pm }}=\left[-\cos \alpha \cos \beta \pm \cos \alpha \sin \beta 0\right]$ and ${\mathbf{P}}_{\text{xoy}}^{{\mathbf{n}}^{\pm }}=\left[\pm \cos \alpha \sin \beta \cos \alpha \cos \beta 0\right]$.

2.2 Microstructure-Based Hyperelastic Strain Energy Density Function

Given the microstructural morphology as mentioned above, corneal stroma after CXL can be taken as a fiber reinforced composite material that has the typical mechanical properties of biological soft materials, e.g., quasi-incompressibility, anisotropy and hyperelasticity (Pandolfi and Manganiello, 2006; Liu et al., 2020a). These biomechanical properties can be well characterized by a unified strain energy density function $\Psi =\Psi \left(\mathbf{C}\right)$, where $\mathbf{C}={\mathbf{F}}^{\text{T}}\mathbf{F}$ is the right Cauchy-Green deformation tensor, and $\mathbf{F}$ is the deformation gradient. In order to address the numerical issue associated with quasi-incompressibility, the strain energy density function for corneal materials is generally decoupled into an isochoric ${\Psi }_{\text{dev}}$ and volumetric ${\Psi }_{\text{vol}}$ part (Gasser et al., 2006), i.e.

where ${\Psi }_{\text{vol}}\left(J\right)=K\left[\left(J^2-1\right)/2-\ln J\right]$ with K the bulk modulus, and $J=\det \left(\mathbf{F}\right)$ is the volume ratio. Here, det is the symbol for the determinant operator. After the CXL treatment, the isochoric part ${\Psi }_{\text{dev}}$ contains three dominant components, i.e., the depth-dependent and irradiation dose-dependent cross-links, statistically distributed collagen fibers and the isotropic matrix, which determine the biomechanical behaviors of cross-linked corneas. Thus, the expression of ${\Psi }_{\text{dev}}$ can be specified as

where ${\Psi }_{\text{C}}$, ${\Psi }_{\text{F}}$ and ${\Psi }_{\text{M}}$ are corresponding to the cross-links, collagen fibers and matrix, respectively. $\overline{\mathbf{C}}={\overline{\mathbf{F}}}^{\text{T}}\overline{\mathbf{F}}$ with $\overline{\mathbf{F}}=J^{-1/3}\mathbf{F}$ is the modified right Cauchy-Green deformation tensor that can be characterized by three modified principal invariants ${\overline{I}}_1$, ${\overline{I}}_2$ and ${\overline{I}}_3$, i.e., ${\overline{I}}_1=\overline{\mathbf{C}}:\mathbf{I}=\text{tr}\overline{\mathbf{C}}$, ${\overline{I}}_2=\frac{1}{2}\left[{\left(\text{tr}\overline{\mathbf{C}}\right)}^2-\text{tr}{\left(\overline{\mathbf{C}}\right)}^2\right]$ and ${\overline{I}}_3=\det \overline{\mathbf{C}}$. Moreover, one can also obtain four additional dimensionless invariants ${\overline{I}}_4$, ${\overline{I}}_6$, ${\overline{I}}_{{\text{m}}^{\pm }}$ and ${\overline{I}}_{{\text{n}}^{\pm }}$, which correspond to the deformation of the two families of orthogonally distributed fibers and adherent cross-links as illustrated in Figure 1B, i.e., ${\overline{I}}_4={\mathbf{a}}_1\cdot \overline{\mathbf{C}}\cdot {\mathbf{a}}_1$, ${\overline{I}}_6={\mathbf{a}}_2\cdot \overline{\mathbf{C}}\cdot {\mathbf{a}}_2$, ${\overline{I}}_{m^{\pm }}={\mathbf{m}}^{\pm }\cdot \overline{\mathbf{C}}\cdot {\mathbf{m}}^{\pm }$ and ${\overline{I}}_{n^{\pm }}={\mathbf{n}}^{\pm }\cdot \overline{\mathbf{C}}\cdot {\mathbf{n}}^{\pm }$.

In the following, let us focus on the deduction of ${\Psi }_{\text{C}}$ by considering three dominant physical mechanisms. 1) First, according to the clinical observation of CXL, it is noticed that the density of cross-links is mainly dominated by the irradiation dose $D$ (Caporossi et al., 2013). In order to quantitatively characterize the variation of the cross-link density as a function of $D$ for corneal strips with a limited thickness, a nonlinear relationship has been proposed for the normalized density $\rho$ of cross-links in our recent work (Xiao et al., 2022) that $\rho \left(\overline{D}\right)={\overline{D}}^m$ with m (m > 0) a dimensionless parameter and $\overline{D}=D/D_{\text{r}}$. Here, $D_{\text{r}}=1{\text{J/cm}}^2$ is a reference irradiation dose. 2) Second, the distribution of cross-links induced by riboflavin and UVA irradiation is inhomogeneous within the corneal stroma. As the cross-linking reaction prefers to occur in the corneal anterior layer, cross-links will be firstly formed at the anterior layer of the stroma. With the increase of corneal thickness h along the irradiation direction, $\rho$ gradually decreases to be zero that separates the stroma layer into two different regions with and without cross-links, respectively (Mazzotta et al., 2007; Spoerl et al., 2007; Mazzotta et al., 2018). In order to further address the evolution of $\rho$ as a function of h, the expression of $\rho$ can be revised as

where $N\left(\overline{h}\right)$ is a normalized UVA energy function ranging from 0 to 1, which represents the variation of the absorption degree of UVA energy as a function of the corneal thickness. $\overline{h}=h/h_{\text{r}}$ is the normalized corneal thickness, and $h_{\text{r}}=1\text{mm}$ is a reference corneal thickness. With the assistance of Eq. 3, one can simultaneously address the influence of irradiation dose and depth-dependent cross-link density on the strengthening behavior of corneal stroma after CXL. (3) Third, one may also recall that the cross-links are formed symmetrically with respect to the collagen fibers, and 66% of the fibers distribute along the N-T and S-T directions while the rest are randomly distributed in the matrix. Therefore, two different parts of cross-links should be considered for the expression of ${\Psi }_{\text{C}}$, i.e., the one with random orientations in the matrix and the one with specific orientations (${\mathbf{m}}^{\pm }$ and ${\mathbf{n}}^{\pm }$) along the N-T and S-I directions of collagen fibers. In this way, one can finally have ${\Psi }_{\text{C}}$ as

where the first part of the expression in the square bracket indicates that the homogeneously distributed cross-links in the matrix follow a linearly strengthening law, while the second part represents that the strengthening contribution of cross-links along the ${\mathbf{m}}^{\pm }$ and ${\mathbf{n}}^{\pm }$ directions is characterized in an exponential form. Moreover, $\psi$ is the volume fraction of the cross-links distributed along the ${\mathbf{m}}^{\pm }$ and ${\mathbf{n}}^{\pm }$ directions. L indicates the stiffness modulus of cross-links, and n is a dimensionless parameter.

Comparing with the previous theoretical work (Holzapfel and Ogden, 2020), the advantages of Eq. 4 are obvious: 1) An explicit expression is proposed for the variation of the cross-link density as a function of the corneal thickness and irradiation dose, which can not only facilitate the quantitative characterization of the stiffening and strengthening behavior of cross-linked cornea, but also make it possible to predict the mechanical responses of corneas at various irradiation doses. 2) Consideration of the homogenous and heterogeneous distribution of cross-links in the matrix and along the N-T and S-I directions of fibers, respectively, makes it rational to further address the influence of cross-links on the fundamental deformation mechanisms, i.e., the comparison of Mises stress distribution before and after CXL with finite element simulations.

For corneal stroma without CXL, the biomechanical responses are mainly determined by the anisotropic behaviors of fibers and isotropic properties of the matrix. Considering that parts of the collagen fibers distribute homogeneously in the matrix and the rest along the N-T and S-I directions, ${\Psi }_{\text{F}}$ can then be expressed in a similar form as ${\Psi }_{\text{C}}$, i.e.

where $k_1$ is the stiffness modulus of collagen fibers and $k_2$ is a dimensionless parameter. The first and second term in the square bracket, respectively, represents the stiffening contribution of the randomly distributed and orthogonally orientated fibers. One may note that more sophisticated models have been developed in previous literatures concerning the architecture variation of collagen fibers both in the thickness and planar view (Cheng et al., 2015; Montanino et al., 2018; Gizzi et al., 2021). For instance, in the work of Cheng et al. (2015), the elastic behavior of the stromal collagen network was modelled by considering the probability distribution of three-dimensional collagen orientation for each point in the stroma. Montanino et al. (2018) systematically compared four different models accounting for local variations of the collagen fibril architecture, and numerical results obtained from the simulation of three ideal mechanical tests indicated slight differences between the models in terms of global response and stress distribution. As a comparison with these architecture-based models, the expression of Eq. 5 for collagen fibers is simplified that can not characterize the continuous distribution of materials properties at each point of the cornea, but still can effectively characterize the anisotropic properties of corneal stroma without CXL.

Moreover, the corneal matrix is similar to a polymer that exhibits a quasilinear and hyperelastic isotropic behavior, therefore, the Neo-Hookean model is considered here for ${\Psi }_{\text{M}}$, i.e.

where $C_{10}$ is the shear modulus of the matrix.

By combining Eqs 2–6, the expression of $\Psi \left(\mathbf{C}\right)$ for corneal stroma with gradiently distributed cross-links finally yields as

which involves the contribution of cross-links, collagen fibers and matrix materials, as well as the volumetric part of the strain energy density function. Thereinto, the cross-links and fibers mainly dominate the strength of corneal stroma at a comparatively large deformation, and determine the typical anisotropic biomechanical properties due to the dominant distribution of cross-links and fibers along the N-T and S-I directions. As a comparison, the matrix materials behave homogeneously, and affect the biomechanical behaviors mainly at the initial deformation stage.

2.3 Cauchy Stress Tensor and Elasticity Modulus Tensor

Following the chain rule indicated by (Weiss et al., 1996; Holzapfel et al., 2000; Gasser et al., 2006), it is convenient to further calculate the Cauchy stress tensor $\mathbf{\sigma }$ from Eq. 7 with $\mathbf{\sigma }={\mathbf{F}}^{\text{T}}\mathbf{SF}/J$. Here, $\mathbf{S}=2\partial \Psi /\partial \mathbf{C}$ is the second Piola-Kirchhoff stress tensor. In this way, $\mathbf{\sigma }$ can then be decoupled into the isochoric and volumetric part, as

where ${\mathbf{\sigma }}_{\text{vol}}=K\left(J-1/J\right)\mathbf{{\rm I}}$ and ${\mathbf{\sigma }}_{\text{dev}}=J^{-2/3}\mathrm{\mathbb{P}}:\overline{\mathbf{\sigma }}$ with $\mathrm{\mathbb{P}}=\left(\mathrm{}-\mathbf{I}\otimes \mathbf{I}/3\right)$ the projection tensor. $\mathrm{}$ and $\mathbf{I}$ are, respectively, the fourth-order and second-order unit tensor. Moreover, one can decouple the fictitious Cauchy stress tensor $\overline{\mathbf{\sigma }}$ as

with

where $\overline{\mathbf{B}}=\overline{\mathbf{F}}{\overline{\mathbf{F}}}^{\text{T}}$ is the modified left Cauchy-Green deformation tensor. ${\overline{\mathbf{A}}}_4={\mathbf{a}}_1\cdot \overline{\mathbf{F}}$, ${\overline{\mathbf{A}}}_6={\mathbf{a}}_2\cdot \overline{\mathbf{F}}$, ${\overline{\mathbf{X}}}_{{\mathbf{m}}^{\pm }}={\mathbf{m}}^{\pm }\cdot \overline{\mathbf{F}}$ and ${\overline{\mathbf{X}}}_{n^{\pm }}={\mathbf{n}}^{\pm }\cdot \overline{\mathbf{F}}$. ${\overline{\mathbf{\sigma }}}_{\text{C}}$, ${\overline{\mathbf{\sigma }}}_{\text{F}}$ and ${\overline{\mathbf{\sigma }}}_{\text{M}}$ are, respectively, the Cauchy tensor component corresponding to the cross-links, fibers and matrix. Combining Eqs 8–10, the Cauchy stress tensor for corneal stroma with CXL can then be specified as

For further finite element simulation, the derivation of the fourth-order elasticity tensor is also required. According to its fundamental definition, one can firstly have $ℂ=4{\partial }^2\Psi /\partial \mathbf{C}\partial \mathbf{C}$, and by using the push-forward operator, the elasticity tensor is then transformed as ${ℂ}^J=\mathbf{FF}ℂ{\mathbf{F}}^{\text{T}}{\mathbf{F}}^{\text{T}}/J$. Under the deformation framework with large strain, the Zaremba-Jaumann rate is considered for the stress change rate (Huang et al., 2016; Fehervary et al., 2020). Then, the elasticity tensor related to the Jaumann rate gives as

where $\delta$ is the Kronecker delta and

2.4 Finite Element Simulations for the Inflation Tests of Cornea With CXL

In this work, the proposed hyperelastic constitutive equations for corneal stroma with gradiently distributed cross-links are implemented into the subroutine UMAT of ABAQUS to simulate the inflation tests of porcine cornea before and after cross-linking. We firstly introduce the details about the establishment of the finite element model, which contains the basic parameters for the geometrical model, boundary conditions and element mesh, etc. Then, the setup of the gradiently distributed cross-links in the model is introduced by referring to corresponding experimental data. Moreover, the dominant steps for the calibration of model parameters are specified by comparing the simulated results with the experimental data under different UVA irradiation doses.

2.4.1 Geometrical Model and Boundary Conditions

The 3D geometrical parameters, as listed in Table 1, are informed from the experimental measurements of porcine cornea (Cornaggia et al., 2020), which include the anterior diameter $\left(D_{\text{NT}}^{\text{A}}\right)$, posterior diameter $\left(D_{\text{NT}}^{\text{P}}\right)$, anterior curvature $\left(R_{\text{NT}}^{\text{A}}\right)$ and posterior curvature $\left(R_{\text{NT}}^{\text{P}}\right)$ along the N-T meridian; the anterior diameter $\left(D_{\text{SI}}^{\text{A}}\right)$, posterior diameter $\left(D_{\text{SI}}^{\text{P}}\right)$, anterior curvature $\left(R_{\text{SI}}^{\text{A}}\right)$ and posterior curvature $\left(R_{\text{SI}}^{\text{P}}\right)$ along the S-I meridian; the edge thickness $t_0$, central thickness $t_{\text{c}}$ and height of the anterior chamber $H$. The center of the anterior curvature along the N-T meridian locates at the origin of coordinates, and the thickness direction is consistent with the z axis. Correspondingly, the S-I and N-T direction for 66% of the fiber distribution are parallel to the x and y axis, respectively. During the inflation tests, corneal deformation is limited by the adjacent sclera, therefore, the edge of porcine cornea is fixed while the rest surfaces are free of constraints. IOP is applied on the posterior surface of the cornea, which gradually increases up to 30 mmHg as the same considered in the experimental work (Cornaggia et al., 2020). The 3D mesh for the cornea is illustrated in Figure 2, where the geometrical model is meshed with 32,400 hexahedral elements of type C3D8 and 36,571 nodes. Numerical convergence is ensured when the size of the mesh is set as 0.1 mm along the thickness direction and 0.3 mm for the other directions.

TABLE 1. Geometrical parameters (unit in mm) for the finite element simulation of porcine cornea (Cornaggia et al., 2020).

FIGURE 2. 3D schematic (N-T meridional section) of the geometrical model for the inflation test of porcine cornea.

2.4.2 Numerical Implementation of the Gradiently Distributed Cross-Links

In order to numerically implement the gradiently distributed cross-links in the finite element model, one should firstly determine the normalized UVA energy function $N\left(\overline{h}\right)$ as expressed in Eq. 3. According to the experimental work of (Liu et al., 2020b), the UVA energy is noticed to get saturated near the anterior surface of the cornea, and then gradually decreases with increasing corneal thickness until a critical depth is reached that the UVA energy becomes zero. By normalizing the distribution data of the UVA energy with the maximum UVA energy, the experimental version of $N\left(\overline{h}\right)$ is then obtained, as illustrated by the dots in Figure 3A. By fitting these normalized experimental data in the layers with/without cross-links, the numerical expression of $N\left(\overline{h}\right)$ yields as

FIGURE 3. (A) Normalized UVA energy compared between the experimental data (Liu et al., 2020b) and fitting function. (B) Numerical distribution of the normalized UVA energy function in ABAQUS.

Here, one should note that Eq. 14 is a specific fitting law corresponding to the data of (Liu et al., 2020b), and it can be changed depending on the research objects and irradiation conditions. Moreover, it is a simplified treatment to take the same gradient distribution form for the whole cornea. In practice, this initial condition can hardly be guaranteed due to the various diffusion abilities of riboflavin and irradiation conditions. However, it is comparatively difficult to obtain the actual cross-link distribution profile for the whole cornea from microstructural observations so far.

Next, according to the coordinates of the Gauss points of each element, the normalized cross-link energy for each element is obtained by submitting $h=R-\sqrt{x^2+y^2+z^2}$ into Eq. 14. Here, x, y and z are the coordinates of each Gauss point. In this way, the inhomogeneous distribution of the normalized cross-link energy can be effectively incorporated into the subroutine UMAT, Figure 3B for an illustration. By further considering the specific irradiation dose in Eq. 3, the gradiently distributed cross-links are finally obtained in ABAQUS.

2.4.3 Parameter Calibration of the Constitutive Equations

The constitutive equations presented in Eq. 7 contain nine model parameters. Thereinto, $C_{10}$ for the matrix materials, $k_1$ and $k_2$ for the collagen fibers, and $m$, $n$ and $L$ for the cross-links need to be determined by comparing the numerical results with corresponding experimental data. The remaining three can be obtained from previous literature. For instance, $\psi =66\%$ represents the volume ratio of the fibers distributed along the N-T and S-I directions, which is informed by the X-ray diffraction data (Aghamohammadzadeh et al., 2004; Meek et al., 2005; Meek and Knupp, 2015). $K=100\text{ MPa}$ is considered to ensure the numerical convergence (Huang et al., 2016). $\beta$ is the angle between the cross-links and fibers, and is assumed to be 45° (Holzapfel and Ogden, 2020). In the following, we describe a detailed parameter calibration method for the determination of $C_{10}$, $k_1$, $k_2$, $m$, $n$ and $L$.

1) Perform inflation tests for the cornea without CXL to obtain the experimental IOP-apex displacement (P-r) relationships. Compare the simulated results with the experimental data to obtain $C_{10}$, $k_1$ and $k_2$.

2) Perform inflation tests at different UVA irradiation doses to get the experimental P-r relationships treated by CXL. Then, transform the P-r relationships with/without CXL into the equivalent stress-strain $\left(\sigma -\epsilon \right)$ relations following the conversion formula proposed by (Anderson et al., 2004), i.e.

where $t=\left(t_{\text{o}}+t_c\right)/2$ is the average corneal thickness, and $\upsilon =0.50$ is the Poisson’s ratio for the cornea. Here, one should note that the P-r and $\sigma -\epsilon$ relations are actually equivalent but only presented in different forms to help calibrate the model parameters affected by CXL.

3) Consider the transformed ${\text{σ}}^{D_i}-\epsilon$ relations at three irradiation doses $D_i$ ($i=\mathrm{1,2,3}$ and $D_1<D_2<D_3$). Then, obtain the Cauchy stress component related to the cross-links ${\text{σ}}_{\text{C}}^{D_i}$ through ${\text{σ}}_{\text{C}}^{D_i}={\text{σ}}^{D_i}-{\text{σ}}^{D_0}$ with $D_0=0$ representing the condition without CXL.

4) Plot the ${\text{σ}}_{\text{C}}^{D_3}/{\text{σ}}_{\text{C}}^{D_1}-\epsilon$ and ${\text{σ}}_{\text{C}}^{D_2}/{\text{σ}}_{\text{C}}^{D_1}-\epsilon$ relations, and a straight line is expected according to Eq. 3. Then, one can have ${\text{σ}}_{\text{C}}^{D_2}/{\text{σ}}_{\text{C}}^{D_1}={\left(D_2/D_1\right)}^{m_1}$ and ${\text{σ}}_{\text{C}}^{D_3}/{\text{σ}}_{\text{C}}^{D_1}={\left(D_3/D_1\right)}^{m_2}$, which yield the values of $m_1$ and $m_2$. Take the average value of $m_1$ and $m_2$ to finally obtain the value of parameter $m$.

5) Perform finite element simulation for the inflation tests at one given irradiation dose, e.g., $D_1$. Compare the simulated results with related experimental data to obtain the fitting values of $n$ and $L$.

3 Results

3.1 IOP-Apex Displacement Relations of Cornea at Different Irradiation Doses

We firstly apply the proposed constitutive laws and finite element model as presented in Section 2 to simulate the inflation test of porcine corneas. As informed by the experimental work of (Boschetti et al., 2021), 90 excised and de-epithelized porcine corneas are considered in the test, and the irradiation dose ranges from 0 $\text{J}/{\text{cm}}^2$ to 10.8 $\text{J}/{\text{cm}}^2$ by using a constant irradiance of 9 mW/cm². Following the model parameter calibration process as mentioned in Section 2.4.3, we can not only obtain the transformed experimental data as illustrated in Figure 4, but also get the model parameters with and without the influence of CXL, as listed in Table 2. Thereinto, the values of $C_{10}$ and $k_1$ are close to the effective range for porcine cornea (Cornaggia et al., 2020). With the calibrated constitutive equations, we are then able to compare the P-r relations at different irradiation doses. In Figure 5, the evolution of IOP with respect to the apex displacement is compared between the simulated results (lines) and experimental data (dots) when $D=0\text{ J}/{\text{cm}}^2$ and $2\text{.}7\text{ J}/{\text{cm}}^2$. It can be seen that: 1) Both the simulated results with and without CXL effect match well with corresponding experimental data when IOP increases up to 30 mmHg. 2) At the same value of IOP, the apex displacement of corneal materials after cross-linking is obviously much smaller than that without cross-linking. Corresponding mechanisms are mainly ascribed to the formation of covalent bonds between adjacent proteoglycan and collagen molecules, which make the cross-linked cornea be more efficient in withstanding the external load. Therefore, cornea after CXL can better resist the geometrical deformation when stimulated by the increase of IOP.

FIGURE 4. (A) Experimental data for the IOP-apex displacement (P-r) relationships of porcine cornea at different irradiation doses (Boschetti et al., 2021). (B) Equivalent stress-strain ($\sigma -\epsilon$) relations at various irradiation doses following the conversion formula of (Anderson et al., 2004). (C) Evolution of ${\text{σ}}_{\text{C}}^{D_i}$ ($i=\mathrm{1,2,3}$) with respect to the strain $\epsilon$ at three irradiation doses, i.e., $D_1=2\text{.}7\text{ J}/{\text{cm}}^2$, $D_2=8\text{.}10\text{ J}/{\text{cm}}^2$ and $D_3=10\text{.}8\text{ J}/{\text{cm}}^2$. (D) Variation of ${\text{σ}}_{\text{C}}^{D_3}/{\text{σ}}_{\text{C}}^{D_1}$ and ${\text{σ}}_{\text{C}}^{D_2}/{\text{σ}}_{\text{C}}^{D_1}$ as a function of $\epsilon$.

TABLE 2. Calibrated model parameters for porcine cornea with CXL.

FIGURE 5. Comparison of the P-r relations of porcine cornea between the experimental data (Boschetti et al., 2021) and simulated results when $D=0\text{ J}/{\text{cm}}^2$ and $2\text{.}7\text{ J}/{\text{cm}}^2$.

In order to further verify the rationality and accuracy of the proposed constitutive equations and calibrated parameters, we also try to compare the predicted P-r results with corresponding experimental data (Cornaggia et al., 2020)by submitting $D=8.10\text{ J}/{\text{cm}}^2$ and $10\text{.}8\text{ J}/{\text{cm}}^2$ into the calibrated model. As shown in Figure 6, a rational agreement is still achieved for both the fitted and predicted results when $D$ increases from $2\text{.}7\text{ J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$, which obviously indicates that the hyperelastic constitutive equations developed in this work can well characterize the biomechanical behaviors of corneal stroma with gradiently distributed cross-links for a wide range of irradiation dose. In addition, for the same value of IOP, the apex displacement becomes further smaller with increasing irradiation dose, and it is determined by the inhomogeneous increment of cross-link density from the anterior layer of the corneal stroma.

FIGURE 6. Comparison of the P-r relations of porcine cornea between the experimental data (Boschetti et al., 2021), predicted and simulated results when $D=2\text{.}7\text{ J}/{\text{cm}}^2$, $8\text{.}10\text{ J}/{\text{cm}}^2$ and $10\text{.}8\text{ J}/{\text{cm}}^2$.

4 Discussion

4.1 Distribution of the Mises Stress at Various Irradiation Doses

After the analysis of macroscopic materials deformation, we can move forward to discuss the underlying deformation mechanisms with the assistance of finite element simulations. As informed from previous literature, the distribution of cross-links in the corneal stroma is inhomogeneous. In order to get a sophisticated comprehension of the variation of corneal strength at different corneal thicknesses, the Mises stress (unit in MPa) distribution on the anterior, middle and posterior surface of corneal stroma is compared at various irradiation doses when the IOP equals 30 mmHg, Figure 7 for an illustration. It is noticed that: 1) For the cornea without CXL (Figures 7A–C), the region with the maximum stress locates on the intersecting line between the side and posterior surface. This is consistent with the heterogeneous distribution of fibers and geometrical feature of the corneal stroma, which makes the x axis direction the main load bearing direction during the inflation test. As a comparison, the region with the maximum stress becomes the anterior surface for the cornea after CXL (Figures 7D,G,J). This feature is ascribed to the special microstructure distribution that cross-links with a maximum value of density locate on the anterior region, and makes it become the dominant load bearing tissue. 2) On the middle surface, the distribution of Mises stress varies from an elliptical pattern (Figures 7B,E) to a square pattern (Figures 7H,K) with the increase of irradiation dose. The primary cause originates from the fact that the angle between the cross-links and fibers along the N-T and S-I directions is set as 45° for corneal stroma. Therefore, the dominant load bearing directions gradually vary from the two main fiber distribution directions to the two main cross-link distribution directions when the irradiation dose increases. 3) As a comparison, the Mises stress on the posterior surface tends to decrease with the increase of irradiation dose (Figures 7C,F,I,L), which is mainly determined by the inhomogeneous distribution of cross-links along the thickness direction, and no cross-links exist on the posterior surface of cornea stroma.

FIGURE 7. Comparison of the Mises stress (unit in MPa) distribution for the anterior surface (A,D,G,J), middle surface (B,E,H,K) and posterior surface (C,F,I,L) of porcine cornea when the IOP equals 30 mmHg. The irradiation dose increases from $\text{0 J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$.

Figure 8 illustrates the meridional N-T sectional view of the Mises stress (unit in MPa) distribution for corneal stroma when the irradiation dose increases from $\text{0 J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$, and the IOP equals 30 mmHg. When $D=\text{0 J}/{\text{cm}}^2$, a comparatively homogeneous distribution of the Mises stress (Figure 8A) is noticed in the central region of corneal stroma. As a comparison, when D increases from $2\text{.}7\text{ J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$, an obviously gradient stress distribution is observed in Figures 8B–D that the Mises stress decreases from the anterior layer surface along the thickness direction until it becomes zero. This thickness-dependent distribution of Mises stress turns to be increasingly obvious with the increase of irradiation dose. Recall Eq. 3 for the depth-dependent density of cross-links, one may note that the increment rate of the cross-link density with the irradiation dose is the fastest in the anterior layer when $m>0$. Therefore, the gradient change of materials strength becomes increasingly obvious for the cross-linked corneas with a comparatively high value of irradiation dose.

FIGURE 8. Meridional N-T sectional view (A–D) of the Mises stress (unit in MPa) distribution for porcine cornea when the IOP equals 30 mmHg. The irradiation dose increases from $\text{0 J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$.

4.2 Comparison of Different Strengthening Mechanisms Affected by CXL

In order to get a sophisticated comprehension of the contribution of different microstructures to the materials strength of cross-linked corneas, the Mises stress distribution (unit in MPa) on the middle surface of corneal stroma is compared between the matrix, fiber and cross-link component. The irradiation dose increases from $2\text{.}7\text{ J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$, and the IOP equals 30 mmHg (corresponding to comparatively large deformations). It can be seen in Figure 9 that for corneas with CXL, the distribution pattern intrinsically varies for the three different microstructures, namely, the Mises stress of the matrix, fiber and cross-link component exhibits a circular, zonal and square pattern, respectively. This is fundamentally determined by the various organizations of microstructures at the micro-scale, i.e., the matrix materials are homogenously distributed in corneal stroma, while the fibers are heterogeneously distributed in the N-T and S-I directions, and the cross-links are preferentially distributed along the two main families of fibers with a specific angle. Corresponding to these microstructural distribution features, both the N-T and S-I direction, as well as the 45° direction with respect to the fibers become the main load bearing directions. Moreover, with the increase of irradiation dose, the square pattern of cross-links becomes increasingly obvious, while the circle pattern of the matrix changes to be square, and the zonal pattern of the fibers gets weakened and finally disappears, which is ascribed to the fact that the density of cross-links keeps increasing that makes the materials strength of corneas after CXL be dominated by the riboflavin/UVA induced cross-links.

FIGURE 9. Comparison of the Mises stress (unit in MPa) distribution on the middle surface of porcine cornea between the matrix (A,D,G), fiber (B,E,H) and cross-link (C,F,I) component. The IOP equals 30 mmHg, and the irradiation dose increases from $2\text{.}7\text{ J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$.

4.3 Displacement Distribution at Various Irradiation Doses

Last but not the least, the effect of CXL on the deformation behavior of corneal materials can also be well characterized through the analysis of the distribution pattern of macroscopic displacement. As illustrated in Figure 10, the meridional N-T sectional view of the displacement (unit in mm) distribution for porcine cornea is compared before and after the expansion deformation. The IOP equals 30 mmHg, and the irradiation dose varies from $\text{0 J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$. It is known that: 1) For the cornea with/without CXL, materials deformation mainly concentrates at the center part due to the constraint conditions applied at the corneal edges, Figure 2 for an illustration. 2) With the increase of irradiation dose, the apex displacement at 30 mmHg decreases from 0.257 to 0.106 mm, which clearly indicates that the ability of corneal stroma to resist the inflation deformation is effectively elevated after cross-linking. This deformation feature is fundamentally consistent with the Mises stress distribution patterns as presented in Figures 7–9 that the increase of cross-link density at a comparatively high value of irradiation dose can help improve the materials strength. From the clinic point of view, the ability to control corneal deformation and increase the corneal stability through the application of CXL offers a promising method to cure the congenital corneal diseases like keratoconus.

FIGURE 10. Meridional N-T sectional view (A–D) of the displacement (unit in mm) distribution for porcine cornea when the IOP equals 30 mmHg. The irradiation dose increases from $\text{0 J}/{\text{cm}}^2$ to $10\text{.}8\text{ J}/{\text{cm}}^2$.

Conclusion

The method of CXL has long been realized as an effective technique to medically improve the stiffness and strength of collagen-based biomaterials, e.g., corneal stroma. However, corresponding deformation mechanisms from a biomechanical point of view have not been well comprehended so far, especially for the influence of gradiently distributed cross-links on the thickness-dependent strengthening behavior of corneal stroma after CXL. In this work, both theoretical model and finite element simulation are combined to help address the fundamental deformation mechanisms of cross-linked corneal stroma under inflation tests. The dominant innovations and conclusions are as follows:

1) Within the framework of strain energy density function theory, a hyperelastic model is proposed for corneal stroma that contains the stiffening contributions of gradiently distributed cross-links, collagen fibers and matrix materials. Thereinto, the strengthening mechanism induced by CXL is characterized by a mechanism-based strain energy function, which not only explicitly depends on the density of cross-links as a function of the corneal thickness and UVA irradiation dose, but also relies on the discrepant distribution of cross-links in the materials matrix and along the orthogonally distributed fibers.

2) The developed constitutive laws are incorporated into the subroutine UMAT of ABAQUS to simulate the deformation behaviors of corneal inflation tests at different UVA irradiation doses. Related methods for the setup of the heterogeneous distribution of cross-links along the thickness direction of corneal stroma and calibration of model parameters based on experimental data are both specified for the numerical simulation. It is realized that the simulated results not only can fit well with corresponding experimental data with/without CXL, but also are able to predict the macroscopic IOP-apex displacement relations as a function of the UVA irradiation dose. Further attentions, from the aspects of clinical practice, will be focused on the analysis of the effect of CXL on the mechanical behaviors of diseased corneas, e.g. keratoconus and corneas after refractive surgery.

3) The established finite element simulation framework offers a valid platform to address the underlying stiffening mechanisms induced by different microstructures at different UVA irradiation conditions. For the case without CXL, the orthogonally distributed fibers along the N-T and S-I directions dominate the heterogeneous deformation behaviors. As a comparison, for the condition with CXL, the stiffening behavior induced by cross-links dramatically varies along the thickness direction of corneal stroma. The unveiled mechanisms may theoretically help design and guide the operation of CXL in clinics.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

Manuscript writing and review: CX, XX, and YY. Model development and finite element simulation: CX. Design of this study and revision of the manuscript: XX and YY. All authors read and approved the final manuscript.

Funding

This work was supported by the National Nature Science Foundation of China (NSFC) under Contract No’s 11802344 and 12172384. XX thanks the initial funding supported by Central South University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References