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 a 1 and a 2 with a 1 = [ − cos ⁡ α 0 ⁡ sin ⁡ α ] and a 2 = [ 0 ⁡ cos ⁡ α ⁡ sin ⁡ α ] . Correspondingly, the projection vectors of a 1 and a 2 on the projection plane (i.e., the x-o-y plane) give as P xoy a 1 = [ − cos ⁡ α 00 ] and P xoy a 2 = [ 0 ⁡ cos α 0 ] .

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 m ± and n ± with m ± = [ − cos ⁡ α ⁡ cos ⁡ β ± cos ⁡ α ⁡ sin ⁡ β ⁡ sin ⁡ α ] and n ± = [ ± cos ⁡ α ⁡ sin ⁡ β ⁡ cos ⁡ α ⁡ cos ⁡ β ⁡ sin ⁡ α ] , respectively. Similarly, the projection vectors of m ± and n ± on the x-o-y plane yield as P xoy m ± = [ − cos ⁡ α ⁡ cos ⁡ β ± cos ⁡ α ⁡ sin ⁡ β 0 ] and P xoy n ± = [ ± cos ⁡ α ⁡ sin ⁡ β ⁡ cos ⁡ α ⁡ cos ⁡ β 0 ] .

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 Ψ = Ψ ( C ) , where C = F T F is the right Cauchy-Green deformation tensor, and 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 Ψ dev and volumetric Ψ vol part (Gasser et al., 2006), i.e. Ψ ( C ) = Ψ dev ( C ¯ , a 1 , a 2 , m ± , n ± ) + Ψ vol ( J ) (1) where Ψ vol ( J ) = K [ ( J 2 − 1 ) / 2 − ln ⁡ J ] with K the bulk modulus, and J = det ( F ) is the volume ratio. Here, det is the symbol for the determinant operator. After the CXL treatment, the isochoric part Ψ 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 Ψ dev can be specified as Ψ dev ( C ¯ , a 1 , a 2 , m ± , n ± ) = Ψ C ( I ¯ 1 , I ¯ m ± , I ¯ n ± ) + Ψ F ( I ¯ 1 , I ¯ 4 , I ¯ 6 ) + Ψ M ( I ¯ 1 ) (2) where Ψ C , Ψ F and Ψ M are corresponding to the cross-links, collagen fibers and matrix, respectively. C ¯ = F ¯ T F ¯ with F ¯ = J − 1 / 3 F is the modified right Cauchy-Green deformation tensor that can be characterized by three modified principal invariants I ¯ 1 , I ¯ 2 and I ¯ 3 , i.e., I ¯ 1 = C ¯ : I = tr C ¯ , I ¯ 2 = 1 2 [ ( tr C ¯ ) 2 − tr ( C ¯ ) 2 ] and I ¯ 3 = det C ¯ . Moreover, one can also obtain four additional dimensionless invariants I ¯ 4 , I ¯ 6 , I ¯ m ± and I ¯ n ± , which correspond to the deformation of the two families of orthogonally distributed fibers and adherent cross-links as illustrated in Figure 1B, i.e., I ¯ 4 = a 1 ⋅ C ¯ ⋅ a 1 , I ¯ 6 = a 2 ⋅ C ¯ ⋅ a 2 , I ¯ m ± = m ± ⋅ C ¯ ⋅ m ± and I ¯ n ± = n ± ⋅ C ¯ ⋅ n ± .

In the following, let us focus on the deduction of Ψ 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 ρ of cross-links in our recent work (Xiao et al., 2022) that ρ ( D ¯ ) = D ¯ m with m (m > 0) a dimensionless parameter and D ¯ = D / D r . Here, D r = 1 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, ρ 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 ρ as a function of h, the expression of ρ can be revised as ρ ( D ¯ , h ¯ ) = [ D ¯ N ( h ¯ ) ] m (3) where N ( h ¯ ) 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. h ¯ = h / h r is the normalized corneal thickness, and h r = 1 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 Ψ C , i.e., the one with random orientations in the matrix and the one with specific orientations ( m ± and n ± ) along the N-T and S-I directions of collagen fibers. In this way, one can finally have Ψ C as Ψ C ( I ¯ 1 , I ¯ m ± , I ¯ n ± ) = ρ ( D ¯ , h ¯ ) L 2 n [ 2 ( 1 − ψ ) ( I ¯ 1 − 3 ) + ψ ∑ i = m ± , n ± ( e n ( I ¯ i − 1 ) 2 − 1 ) ] (4) 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 m ± and n ± directions is characterized in an exponential form. Moreover, ψ is the volume fraction of the cross-links distributed along the m ± and n ± 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, Ψ F can then be expressed in a similar form as Ψ C , i.e. Ψ F ( I ¯ 1 , I ¯ 4 , I ¯ 6 ) = k 1 2 k 2 [ 2 ( 1 − ψ ) ( I ¯ 1 − 3 ) + ψ ∑ i = 4,6 ( e k 2 ( I ¯ i − 1 ) 2 − 1 ) ] (5) 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 Ψ M , i.e. Ψ M ( I ¯ 1 ) = C 10 ( I ¯ 1 − 3 ) (6) where C 10 is the shear modulus of the matrix.

By combining Eqs 2–6, the expression of Ψ ( C ) for corneal stroma with gradiently distributed cross-links finally yields as Ψ ( C ) = ρ ( D ¯ , h ¯ ) L 2 n [ 2 ( 1 − ψ ) ( I ¯ 1 − 3 ) + ψ ∑ i = m ± , n ± ( e n ( I ¯ i − 1 ) 2 − 1 ) ] + k 1 2 k 2 [ 2 ( 1 − ψ ) ( I ¯ 1 − 3 ) + ψ ∑ i = 4,6 ( e k 2 ( I ¯ i − 1 ) 2 − 1 ) ] + C 10 ( I ¯ 1 − 3 ) + K ( J 2 − 1 2 − ln ⁡ J ) (7) 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.

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 σ from Eq. 7 with σ = F T SF / J . Here, S = 2 ∂ Ψ / ∂ C is the second Piola-Kirchhoff stress tensor. In this way, σ can then be decoupled into the isochoric and volumetric part, as σ = σ dev + σ vol (8) where σ vol = K ( J − 1 / J ) Ι and σ dev = J − 2 / 3 ℙ : σ ¯ with ℙ = (  − I ⊗ I / 3 ) the projection tensor.  and I are, respectively, the fourth-order and second-order unit tensor. Moreover, one can decouple the fictitious Cauchy stress tensor σ ¯ as σ ¯ = 2 J F ∂ Ψ ∂ C ¯ F T = 2 J F ( ∂ Ψ C ∂ C ¯ + ∂ Ψ F ∂ C ¯ + ∂ Ψ M ∂ C ¯ ) F T = σ ¯ C + σ ¯ F + σ ¯ M (9) with { σ ¯ C = 2 ρ ( D ¯ , h ¯ ) L ( 1 − ψ ) n B ¯ + ψ ρ ( D ¯ , h ¯ ) L ∑ i = m ± , n ± ( I ¯ i − 1 ) e n ( I ¯ i − 1 ) 2 X ¯ i ⊗ X ¯ i σ ¯ F = 2 k 1 ( 1 − ψ ) k 2 B ¯ + ψ k 1 ∑ i = 4,6 ( I ¯ i − 1 ) e k 2 ( I ¯ i − 1 ) 2 A ¯ i ⊗ A ¯ i σ ¯ M = 2 C 10 B ¯ (10) where B ¯ = F ¯ F ¯ T is the modified left Cauchy-Green deformation tensor. A ¯ 4 = a 1 ⋅ F ¯ , A ¯ 6 = a 2 ⋅ F ¯ , X ¯ m ± = m ± ⋅ F ¯ and X ¯ n ± = n ± ⋅ F ¯ . σ ¯ C , σ ¯ F and σ ¯ 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 σ = 2 [ C 10 + k 1 ( 1 − ψ ) k 2 + ρ ( D ¯ , h ¯ ) L ( 1 − ψ ) n ] ( B ¯ − 1 3 I ¯ 1 I ) + ψ k 1 ∑ i = 4,6 ( I ¯ i − 1 ) e k 2 ( I ¯ i − 1 ) 2 ( A ¯ i ⊗ A ¯ i − 1 3 I ¯ i I ) + ψ ρ ( D ¯ , h ¯ ) L ∑ i = m ± , n ± ( I ¯ i − 1 ) e n ( I ¯ i − 1 ) 2 ( X ¯ i ⊗ X ¯ i − 1 3 I ¯ i I ) + K ( J − 1 J ) I (11)

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 ∂ 2 Ψ / ∂ C ∂ C , and by using the push-forward operator, the elasticity tensor is then transformed as ℂ J = FF ℂ F T F 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 ℂ i j k l σ J = ℂ i j k l J + 1 2 ( δ i k σ j l + δ j k σ i l + δ i l σ k j + δ j l σ i k ) = 1 2 s T ( δ i k B ¯ j l + δ j k B ¯ i l + δ i l B ¯ k j + δ j l B ¯ i k ) − s T ( 2 3 δ k l B ¯ i j − 2 3 δ i j B ¯ k l + 2 9 I ¯ 1 δ i j δ k l ) + ∑ z = 4,6 s F ( I ¯ z ) ( A ¯ z i A ¯ z j − 1 3 I ¯ z δ i j ) ( A ¯ z k A ¯ z l − 1 3 I ¯ z δ k l ) + ∑ z = 4,6 1 2 c F ( I ¯ z ) ( δ i k A ¯ z j A ¯ z l + δ j k A ¯ z i A ¯ z l + δ i l A ¯ z k A ¯ z j + δ j l A ¯ z i A ¯ z k ) − ∑ z = 4,6 c F ( I ¯ z ) ( 2 3 δ k l A ¯ z i A ¯ z j − 2 3 δ i j A ¯ z k A ¯ z l + 2 9 I ¯ z δ i j δ k l ) + ∑ z = m ± , n ± s C ( I ¯ z ) ( X ¯ z i X ¯ z j − 1 3 I ¯ z δ i j ) ( X ¯ z k X ¯ z l − 1 3 I ¯ z δ k l ) + ∑ z = m ± , n ± 1 2 c C ( I ¯ z ) ( δ i k X ¯ z j X ¯ z l + δ j k X ¯ z i X ¯ z l + δ i l X ¯ z k X ¯ z j + δ j l X ¯ z i X ¯ z k ) − ∑ z = m ± , n ± c C ( I ¯ z ) ( 2 3 δ k l X ¯ z i X ¯ z j − 2 3 δ i j X ¯ z k X ¯ z l + 2 9 I ¯ z δ i j δ k l ) + 2 J K δ i j δ k l (12) where δ is the Kronecker delta and { s T = 2 J [ C 10 + k 1 ( 1 − ψ ) k 2 + ρ ( D ¯ , h ¯ ) L ( 1 − ψ ) n ] s F ( I ¯ z ) = ∂ Ψ F ∂ I ¯ z = 2 k 1 J Ψ k 1 ( I ¯ z − 1 ) e k 2 ( I ¯ z − 1 ) 2 s C ( I ¯ z ) = ∂ Ψ C ∂ I ¯ z = 2 J ψ ρ ( D ¯ , h ¯ ) L ( I ¯ i − 1 ) e n ( I ¯ i − 1 ) 2 c F ( I ¯ z ) = ∂ 2 Ψ F ∂ I ¯ z ∂ I ¯ z = 2 k 1 J ψ k 1 [ 1 + 2 k 2 ( I ¯ z − 1 ) 2 ] e k 2 ( I ¯ z − 1 ) 2 c C ( I ¯ z ) = ∂ 2 Ψ C ∂ I ¯ z ∂ I ¯ z = 2 k 1 J ψ ρ ( D ¯ , h ¯ ) L [ 1 + 2 n ( I ¯ z − 1 ) 2 ] e n ( I ¯ z − 1 ) 2 (13)

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.

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 J / cm 2 to 10.8 J / 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  J / cm 2 and 2 . 7  J / 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.

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  J / cm 2 and 10 . 8  J / 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 . 7  J / cm 2 to 10 . 8  J / 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.

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 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 0 J / cm 2 to 10 . 8  J / cm 2 , and the IOP equals 30 mmHg. When D = 0 J / 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 . 7  J / cm 2 to 10 . 8  J / 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.

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 . 7  J / cm 2 to 10 . 8  J / 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.

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 0 J / cm 2 to 10 . 8  J / 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.