Abstract

The cornea, the transparent tissue in the front of the eye, along with the sclera, plays a vital role in protecting the inner structures of the eyeball. The precise shape and mechanical strength of this tissue are mostly determined by the unique microstructure of its extracellular matrix. A clear picture of the 3D arrangement of collagen fibrils within the corneal extracellular matrix has recently been obtained from the secondary harmonic generation images. However, this important information about the through-thickness distribution of collagen fibrils was seldom taken into account in the constitutive modeling of the corneal behavior. This work creates a generalized structure tensor (GST) model to investigate the mechanical influence of collagen fibril through-thickness distribution. It then uses numerical simulations of the corneal mechanical response in inflation experiments to assess the efficacy of the proposed model. A parametric study is also done to investigate the influence of model parameters on numerical predictions. Finally, a brief comparison between the performance of this new constitutive model and a recent angular integration (AI) model from the literature is given.

1 Introduction

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 [59]. 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,2126]. 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 [2733]. 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.

2 Continuum Mechanical Framework

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 [3840].

2.1 Kinematics.

Let $xR$ represent an arbitrary material point in the fixed reference configuration of the body $BR$. The referential body $BR$ undergoes a motion $x=χ(xR,t)$ to the deformed body $Bt$ with the deformation gradient given by
$F=∇χ, and J=detF>0$
(1)
The right and left Cauchy–Green tensors are given by $C=F⊤F$ and $B=FF⊤$, respectively. The deformation gradient admits the polar decomposition, $F=RU$, where R is the rotation and $U=C$ is the stretch. The distortional part of the deformation gradient is
$Fdis=J−1/3F, and detFdis=1$
(2)
The distortional right and left Cauchy–Green deformation tensors are
$Cdis=Fdis⊤Fdis=J−2/3C andBdis=FdisFdis⊤=J−2/3B$
(3)
respectively. We assume there are two families of collagen fibrils in the corneal stroma with their mean referential directions denoted by unit vectors $a04$ and $a06$, respectively. Additionally, we introduce a unit vector an—normal to the plane spanning by $a04$ and $a06$—to identify the out-of-plane direction. The invariants $I¯1, I¯4, I¯6$ and $I¯n$ are written as
$I¯1=trCdis, I¯i=Cdis:a0i⊗a0i for i=4,6and I¯n=Cdis:an⊗an$
(4)
We use the generalized structure tensor Hi to quantify the dispersion of both families of collagen fibrils [35,41]
$Hi=A1+Ba0i⊗a0i+(1−3A−B)an⊗an for i=4,6$
(5)
with constants A and B written as
$A=2κipκop and B=2κop(1−2κip)$
(6)

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.

2.2 Probability Density Functions for Collagen Fibrils With Dispersion.

The detailed collagen fibril microstructural information could be obtained from the SHG images, which fully characterize their in-plane and out-of-plane angular distributions [13,14,26]. Here, the mean orientation of the collagen fibrils at the reference state is represented by a unit vector N in terms of two Eulerian angles $Θ∈[0,2π]$ and $Φ∈[−π/2,π/2]$. We assume that the base vector $e3$ is the out-of-plane direction (see Fig. 1). We use the bivariate Von-Mises distribution function $ρ(Θ,Φ)=ρip(Θ)ρop(Φ)$ to describe the dispersion of collagen fibrils over the unit sphere [35]
$ρip(Θ)= exp [a cos 2(Θ±ξ)]I0(a) andρop(Φ)=22bπ exp [b(cos 2Φ−1)]erf(2b)$
(7)
where a and b denote the concentration parameters for each distribution function $ρip(Θ)$ and $ρop(Φ)$, ξ denotes the angle between the mean collagen fibril orientation and the base vector $e1$, and $I0(a)$ denotes the modified Bessel function of the first kind of order 0. According to Holzapfel et al. [35], both in-plane and out-of-plane dispersion parameters are defined as
$κip=12π∫02πρip(Θ) sin2ΘdΘ, andκop=14∫−π/2π/2ρop(Φ) cos3ΦdΦ$
(8)
Fig. 1
Fig. 1
Close modal
The closed-form relations between dispersion parameters and concentration parameters are obtained from Eqs. (7) and (8)
$κip=12−I1(a)2I0(a), κop=12−18b+142πb exp(−2b)erf(2b)$
(9)

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 a0 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:

Fig. 2
Fig. 2
Close modal
• 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$.

2.3 Free Energy.

The free energy $ψR$ of corneal stroma per unit reference volume is additively decomposed into (1) isotropic contribution from underlying matrix $ψRm$ and (2) anisotropic contribution from two families of collagen fibrils $ψRfi$
$ψR=ψRm(I¯1,J)+∑i=4,6ψRfi(Cdis,Hi)$
(10)
The matrix domain is treated as a nearly incompressible neo-Hookean material
$ψRm=12G0(I¯1−3)+12K(lnJ)2$
(11)

where G0 denotes the ground state shear modulus and K denotes the bulk modulus.

The mechanical response of collagen fibrils is modeled by the following exponential form [34]:
$ψRfi=k12k2(exp(k2(I¯i∗−1)2)−1) for i=4,6$
(12)
where k1 and k2 denote two the stress-like parameters. The distortional generalized invariant $I¯i∗$ is given by
$I¯i∗=tr(HiCdis)=AI¯1+BI¯i+(1−3A−B)I¯n for i=4,6$
(13)

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).

Based on thermodynamic restrictions, the Cauchy stress is then given through
$T=2J−1F∂ψR∂CF⊤=Tm+∑i=4,6Tfi$
(14)
where
$Tm=J−1[G0(Bdis)0+K(lnJ)1]$
(15a)
$(Bdis)0=Bdis−13tr(Bdis)$
(15b)
and
$Tfi=2J−5/3[k1(I¯i∗−1)exp(k2(I¯i∗−1)2)][FHiF⊤−13tr(HiC)1]$
(16)

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

3 Governing Equations

The balance of linear momentum in the deformed body $Bt$ under the equilibrium condition is given by
$divT=0$
(17)
where T the total Cauchy stress given by Eq. (14). The surface traction on the deformed body surface $∂Bt$ is given by
$t(n)=[[T]]n$
(18)

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.

4 Results and Discussion

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.

4.1 Experimental Measurements.

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.

4.2 Geometry and Boundary Conditions.

The first step in any numerical simulation is to define an accurate geometrical representation of the sample. Since we do not have any information about the exact geometry used in inflation tests [37], a generic but popular form is adopted. In particular, we use a biconic surface equation in a cylindrical coordinate system ${Θ,r,x3}$ for both anterior and posterior surfaces of the cornea [24]
$x3−G+r2E1+1−r2F=0$
(19)
with
$E= cos2(Θ−Θx1)Rx1+ sin2(Θ−Θx1)Rx2$
(20)
and
$F=(Qx1+1) cos2(Θ−Θx1)Rx12+(Qx2+1) sin2(Θ−Θx1)Rx22$
(21)

Here, G is the maximum vertical height at r =0, both $Rx1$ and $Rx2$ are the maximum curvatures of the principal meridians along x1 and x2 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 an (solid line).

Fig. 3
Fig. 3
Close modal
For simplicity, we assume that both families of collagen fibrils share the same in-plane dispersion parameter $κip$. Guided by the previous study on X-ray scattering images [9], the spatial distribution of in-plane dispersion is given by [24]
$κip(Θ)=(κipmin+κipmax2)−(κipmax−κipmin2)cos 4Θ$
(22)
where $κipmax=0.5$ and $κipmin=0.1$ are the maximum and minimum value, respectively. After adding the r dependency, Eq. (22) becomes
$κip(Θ,r)=κipmin+12(κip(Θ)−κipmin)(1−cos 2πrRTZ)$
(23)
where $RTZ=5.5$ mm denotes the radius of the transition zone. Note that we assigned a homogeneous in-plane dispersion $κip=0.5$ in the limbus region. The visualization of Eq. (23) is shown in Fig. 4(a). In the process of assigning the out-of-plane parameter $κop$ across the thickness, we used a local coordinate $s∈[0,1]$ parallel to the out-of-plane unit vectors. The local coordinate s =0 at the anterior surface, while s =1 at the posterior surface (see inset plot in Fig. 5). Guided by the SHG image, where the degree of interweaving between collagen fibrils is found to be varying exponentially across the thickness, we link the out-of-plane dispersion parameter $κop$ to the local coordinate s via the following function:
$κop(s)=κopmin+(κopmax−κopmin)(1−exp(−γds))$
(24)

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)).

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

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.

Before running the simulation, one should pay extra attention to the starting point of the simulation. Given that the in vivo measured dimensions $Xphysio$ of the cornea under a physiological loading of $Pphysio=16$ mmHg, we first obtain the stress-free geometry though a zero-pressure algorithm [43,44]. In the algorithm, the mesh connectivity is kept unchanged while the zero-pressure nodal coordinates $Xk+1$ are iteratively updated through
$Xk+1=Xk+(Xkdef−Xphysio)$
(25)

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)).

Fig. 6
Fig. 6
Close modal

4.3 Comparison.

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.

Fig. 7
Fig. 7
Close modal
Table 1

Material parameters used in the simulation

T.I.P.A.I.P.D.P.I.F.V.
G0 (MPa)0.060.060.060.060.060.06
K (MPa)5.55.55.55.55.55.5
k1 (kPa)20202020200.5
k2400400400400400900
κ[24]01/3N/AN/AN/A
$κip$N/AN/AN/AFigure 5(a) 1/2Figure 5(a)
$κop$N/AN/AN/A1/21/2Figure 5(b)
$γd$N/AN/AN/AN/AN/A2.5
T.I.P.A.I.P.D.P.I.F.V.
G0 (MPa)0.060.060.060.060.060.06
K (MPa)5.55.55.55.55.55.5
k1 (kPa)20202020200.5
k2400400400400400900
κ[24]01/3N/AN/AN/A
$κip$N/AN/AN/AFigure 5(a) 1/2Figure 5(a)
$κop$N/AN/AN/A1/21/2Figure 5(b)
$γd$N/AN/AN/AN/AN/A2.5

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.

Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal

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.

4.4 Parametric Study.

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.

Fig. 10
Fig. 10
Close modal

4.5 Numerical Expensiveness.

In this section, we compare the proposed model against a recent AI model [26], in which the collagen fibril free energy is obtained from the angular integration over the unit sphere Ω
$ψRAI=ψRm+1m∫Ωρ(a0)w(λf)dΩ$
(26)

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.

Fig. 11
Fig. 11
Close modal

5 Concluding Remarks

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].

Funding Data

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

Appendix A: Verification of Finite Element Implementation

To verify the numerical implementations, we compare our simulated results with the analytically tractable solutions. We prescribe a simple shear motion to a matrix cube embedded with a family of collagen fibril with a mean direction of $a04=[ax,ay,0]⊤$ and an out-of-plane direction of $an=[−ay,ax,0]⊤$ (see Fig. 12(a)). To make sure they are unit vectors, the condition of $ax2+ay2=1$ has to be fulfilled. According to Gurtin [47], the simple shear deformation is given by
$[F]=[1γ0010001], [B]=[1+γ2γ0γ10001][C]=[1γ0γγ2+10001]$
(A1)
with $γ=tan θ$ denoting the amount of shear. Referring to Eq. (5), the generalized structure tensor is given by
$[H]=[A+Bax2+Cay2Daxay0DaxayA+Bay2+Cax2000A]$
(A2)
with $C=1−3A−B$ and $D=2B+3A−1$. Next, the generalized invariant in Eq. (13) is now given by
$I4∗=2A+Bax2+Cay2+2Daxayγ+(A+Bay2+Cax2)(γ2+1)$
(A3)
Fig. 12
Fig. 12
Close modal
Since the simple shear deformation is a volume preserved motion (i.e., J =1), the Cauchy stress in Eq. (14) could be written as
$T=G0B+2[k1(I4∗−1)exp(k2(I4∗−1)2)][FHF⊤]+P1$
(A4)

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.

Figure 12(b) compares the analytical solutions against the numerical solutions for the shear stress given by
$T12=G0γ+2[k1(I4∗−1)exp(k2(I4∗−1)2)][α+βγ]$
(A5)

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 G0, 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.

References

1.
Hatami-Marbini
,
H.
, and
Etebu
,
E.
,
2013
, “
Hydration Dependent Biomechanical Properties of the Corneal Stroma
,”
Exp. Eye Res.
,
116
, pp.
47
54
.10.1016/j.exer.2013.07.016
2.
Meek
,
K. M.
, and
Knupp
,
C.
,
2015
, “
Corneal Structure and Transparency
,”
Prog. Retinal Eye Res.
,
49
, pp.
1
16
.10.1016/j.preteyeres.2015.07.001
3.
Cogan
,
D. G.
,
1951
, “
Applied Anatomy and Physiology of the Cornea
,”
,
55
, pp.
329
359
.https://pubmed.ncbi.nlm.nih.gov/14835711/
4.
Maurice
,
D. M.
,
1957
, “
The Structure and Transparency of the Cornea
,”
J. Physiology
,
136
(
2
), pp.
263
286
.10.1113/jphysiol.1957.sp005758
5.
Meek
,
K. M.
,
Blamires
,
T.
,
Elliott
,
G. F.
,
Gyi
,
T. J.
, and
Nave
,
C.
,
1987
, “
The Organisation of Collagen Fibrils in the Human Corneal Stroma: A Synchrotron x-Ray Diffraction Study
,”
Curr. Eye Res.
,
6
(
7
), pp.
841
846
.10.3109/02713688709034853
6.
Newton
,
R. H.
, and
Meek
,
K. M.
,
1998
, “
Circumcorneal Annulus of Collagen Fibrils in the Human Limbus
,”
Investigative Ophthalmol. Visual Sci.
,
39
(
7
), pp.
1125
1134
.https://iovs.arvojournals.org/article.aspx?articleid=2161745
7.
Newton
,
R. H.
, and
Meek
,
K. M.
,
1998
, “
The Integration of the Corneal and Limbal Fibrils in the Human Eye
,”
Biophys. J.
,
75
(
5
), pp.
2508
2512
.10.1016/S0006-3495(98)77695-7
8.
,
H.
,
Newton
,
R. H.
, and
Meek
,
K. M.
,
2004
, “
X-Ray Scattering Used to Map the Preferred Collagen Orientation in the Human Cornea and Limbus
,”
Structure
,
12
(
2
), pp.
249
256
.10.1016/j.str.2004.01.002
9.
Boote
,
C.
,
Dennis
,
S.
, and
Meek
,
K.
,
2004
, “
Spatial Mapping of Collagen Fibril Organisation in Primate Cornea-an x-Ray Diffraction Investigation
,”
J. Struct. Biol.
,
146
(
3
), pp.
359
367
.10.1016/j.jsb.2003.12.009
10.
Hayes
,
S.
,
Boote
,
C.
,
Lewis
,
J.
,
Sheppard
,
J.
,
Abahussin
,
M.
,
Quantock
,
A. J.
,
Purslow
,
C.
,
Votruba
,
M.
, and
Meek
,
K. M.
,
2007
, “
Comparative Study of Fibrillar Collagen Arrangement in the Corneas of Primates and Other Mammals
,”
Anatom. Rec.
,
290
(
12
), pp.
1542
1550
.10.1002/ar.20613
11.
Nguyen
,
T.
, and
Boyce
,
B.
,
2011
, “
An Inverse Finite Element Method for Determining the Anisotropic Properties of the Cornea
,”
Biomech. Model. Mechanobiol.
,
10
(
3
), pp.
323
337
.10.1007/s10237-010-0237-3
12.
Abahussin
,
M.
,
Hayes
,
S.
,
Cartwright
,
N. E. K.
,
Kamma-Lorger
,
C. S.
,
Khan
,
Y.
,
Marshall
,
J.
, and
Meek
,
K. M.
,
2009
, “
3d Collagen Orientation Study of the Human Cornea Using X-Ray Diffraction and Femtosecond Laser Technology
,”
Invest. Ophthalmol. Visual Sci.
,
50
(
11
), pp.
5159
5164
.10.1167/iovs.09-3669
13.
Morishige
,
N.
,
Wahlert
,
A. J.
,
Kenney
,
M. C.
,
Brown
,
D. J.
,
Kawamoto
,
K.
,
Chikama
,
T.-I.
,
Nishida
,
T.
, and
Jester
,
J. V.
,
2007
, “
Second-Harmonic Imaging Microscopy of Normal Human and Keratoconus Cornea
,”
Invest. Ophthalmol. Visual Sci.
,
48
(
3
), pp.
1087
1094
.10.1167/iovs.06-1177
14.
Winkler
,
M.
,
Chai
,
D.
,
Kriling
,
S.
,
Nien
,
C. J.
,
Brown
,
D. J.
,
Jester
,
B.
,
Juhasz
,
T.
, and
Jester
,
J. V.
,
2011
, “
Nonlinear Optical Macroscopic Assessment of 3-d Corneal Collagen Organization and Axial Biomechanics
,”
Invest. Ophthalmol. Visual Sci.
,
52
(
12
), pp.
8818
8827
.10.1167/iovs.11-8070
15.
Bryant
,
M. R.
,
Velinsky
,
S. A.
,
Plesha
,
M. E.
, and
Clarke
,
G. P.
,
1987
, “
,”
Eye Contact Lens
,
13
(
4
), pp.
238
242
.https://pubmed.ncbi.nlm.nih.gov/3453772/
16.
Hanna
,
K. D.
,
Jouve
,
F. E.
, and
Waring
,
G. O.
,
1989
, “
Preliminary Computer Simulation of the Effects of Radial Keratotomy
,”
Arch. Ophthalmol.
,
107
(
6
), pp.
911
918
.10.1001/archopht.1989.01070010933044
17.
Pinsky
,
P. M.
, and
Datye
,
D. V.
,
1991
, “
A Microstructurally-Based Finite Element Model of the Incised Human Cornea
,”
J. Biomech.
,
24
(
10
), pp.
907
922
.10.1016/0021-9290(91)90169-N
18.
Bryant
,
M. R.
, and
McDonnell
,
P. J.
,
1996
, “
Constitutive Laws for Biomechanical Modeling of Refractive Surgery
,”
ASME J. Biomech. Eng.
,
118
(
4
), pp.
473
481
.10.1115/1.2796033
19.
Alastrué
,
V.
,
Calvo
,
B.
,
Peña
,
E.
, and
Doblaré
,
M.
,
2006
, “
Biomechanical Modeling of Refractive Corneal Surgery
,”
ASME J. Biomech. Eng.
,
128
(
1
), pp.
150
160
.10.1115/1.2132368
20.
Pandolfi
,
A.
, and
Manganiello
,
F.
,
2006
, “
A Model for the Human Cornea: Constitutive Formulation and Numerical Analysis
,”
Biomech. Model. Mechanobiol.
,
5
(
4
), pp.
237
246
.10.1007/s10237-005-0014-x
21.
Lanir
,
Y.
,
1983
, “
Constitutive Equations for Fibrous Connective Tissues
,”
J. Biomech.
,
16
(
1
), pp.
1
12
.10.1016/0021-9290(83)90041-6
22.
Pinsky
,
P. M.
,
van der Heide
,
D.
, and
Chernyak
,
D.
,
2005
, “
Computational Modeling of Mechanical Anisotropy in the Cornea and Sclera
,”
J. Cataract Refractive Surg.
,
31
(
1
), pp.
136
145
.10.1016/j.jcrs.2004.10.048
23.
Boyce
,
B.
,
Jones
,
R.
,
Nguyen
,
T.
, and
Grazier
,
J.
,
2007
, “
Stress-Controlled Viscoelastic Tensile Response of Bovine Cornea
,”
J. Biomech.
,
40
(
11
), pp.
2367
2376
.10.1016/j.jbiomech.2006.12.001
24.
Pandolfi
,
A.
, and
Holzapfel
,
G. A.
,
2008
, “
Three-Dimensional Modeling and Computational Analysis of the Human Cornea Considering Distributed Collagen Fibril Orientations
,”
ASME J. Biomech. Eng.
,
130
(
6
), p.
061006
.10.1115/1.2982251
25.
Pandolfi
,
A.
,
Fotia
,
G.
, and
Manganiello
,
F.
,
2009
, “
Finite Element Simulations of Laser Refractive Corneal Surgery
,”
Eng. Comput.
,
25
(
1
), pp.
15
24
.10.1007/s00366-008-0102-5
26.
Petsche
,
S. J.
, and
Pinsky
,
P. M.
,
2013
, “
The Role of 3-d Collagen Organization in Stromal Elasticity: A Model Based on X-Ray Diffraction Data and Second Harmonic-Generated Images
,”
Biomech. Model. Mechanobiol.
,
12
(
6
), pp.
1101
1113
.10.1007/s10237-012-0466-8
27.
Sacks
,
M. S.
,
2003
, “
Incorporation of Experimentally-Derived Fiber Orientation Into a Structural Constitutive Model for Planar Collagenous Tissues
,”
ASME J. Biomech. Eng.
,
125
(
2
), pp.
280
287
.10.1115/1.1544508
28.
Driessen
,
N. J.
,
Bouten
,
C. V.
, and
Baaijens
,
F. P.
,
2005
, “
A Structural Constitutive Model for Collagenous Cardiovascular Tissues Incorporating the Angular Fiber Distribution
,”
ASME J. Biomech. Eng.
,
127
(
3
), pp.
494
503
.10.1115/1.1894373
29.
Alastrué
,
V.
,
Martinez
,
M.
,
Menzel
,
A.
, and
Doblaré
,
M.
,
2009
, “
On the Use of Non-Linear Transformations for the Evaluation of Anisotropic Rotationally Symmetric Directional Integrals. Application to the Stress Analysis in Fibred Soft Tissues
,”
Int. J. Numer. Methods Eng.
,
79
(
4
), pp.
474
504
.10.1002/nme.2577
30.
Raghupathy
,
R.
, and
Barocas
,
V. H.
,
2009
, “
A Closed-Form Structural Model of Planar Fibrous Tissue Mechanics
,”
J. Biomech.
,
42
(
10
), pp.
1424
1428
.10.1016/j.jbiomech.2009.04.005
31.
Federico
,
S.
, and
Gasser
,
T. C.
,
2010
, “
Nonlinear Elasticity of Biological Tissues With Statistical Fibre Orientation
,”
J. R. Soc. Interface
,
7
(
47
), pp.
955
966
.10.1098/rsif.2009.0502
32.
Ateshian
,
G. A.
,
Rajan
,
V.
,
Chahine
,
N. O.
,
Canal
,
C. E.
, and
Hung
,
C. T.
,
2009
, “
Modeling the Matrix of Articular Cartilage Using a Continuous Fiber Angular Distribution Predicts Many Observed Phenomena
,”
ASME J. Biomech. Eng.
,
131
(
6
), p.
061003
.10.1115/1.3118773
33.
Gasser
,
T. C.
,
Gallinetti
,
S.
,
Xing
,
X.
,
Forsell
,
C.
,
Swedenborg
,
J.
, and
Roy
,
J.
,
2012
, “
Spatial Orientation of Collagen Fibers in the Abdominal Aortic Aneurysm's Wall and Its Relation to Wall Mechanics
,”
Acta Biomater.
,
8
(
8
), pp.
3091
3103
.10.1016/j.actbio.2012.04.044
34.
Gasser
,
T. C.
,
Ogden
,
R. W.
, and
Holzapfel
,
G. A.
,
2006
, “
Hyperelastic Modelling of Arterial Layers With Distributed Collagen Fibre Orientations
,”
J. R. Soc. Interface
,
3
(
6
), pp.
15
35
.10.1098/rsif.2005.0073
35.
Holzapfel
,
G. A.
,
Niestrawska
,
J. A.
,
Ogden
,
R. W.
,
Reinisch
,
A. J.
, and
Schriefl
,
A. J.
,
2015
, “
Modelling Non-Symmetric Collagen Fibre Dispersion in Arterial Walls
,”
J. R. Soc. Interface
,
12
(
106
), p.
20150188
.10.1098/rsif.2015.0188
36.
Abaqus/Standard
,
2019
,
Abaqus Reference Manuals
,
Dassault Systemes Simulia
,
Providence, RI
.
37.
Anderson
,
K.
,
El-Sheikh
,
A.
, and
Newson
,
T.
,
2004
, “
Application of Structural Analysis to the Mechanical Behaviour of the Cornea
,”
J. R. Soc. Interface
,
1
(
1
), pp.
3
15
.10.1098/rsif.2004.0002
38.
Bosnjak
,
N. S.
,
Wang
,
S.
, and
Chester
,
S. A.
,
2017
, “
Modeling Deformation-Diffusion in Polymeric Gels
,”
Poromechanics VI
,
Paris, France
, July 9–13, pp.
141
148
.10.1061/9780784480779.017
39.
Bosnjak
,
N.
,
Wang
,
S.
,
Han
,
D.
,
Lee
,
H.
, and
Chester
,
S. A.
,
2019
, “
Modeling of Fiber-Reinforced Polymeric Gels
,”
Mech. Res. Commun.
,
96
, pp.
7
18
.10.1016/j.mechrescom.2019.02.002
40.
Holzapfel
,
G. A.
,
2000
,
Nonlinear Solid Mechanics: A Continuum Approach for Engineering
,
Wiley
,
New York
.
41.
Niestrawska
,
J. A.
,
Viertler
,
C.
,
Regitnig
,
P.
,
Cohnert
,
T. U.
,
Sommer
,
G.
, and
Holzapfel
,
G. A.
,
2016
, “
Microstructure and Mechanics of Healthy and Aneurysmatic Abdominal Aortas: Experimental Analysis and Modelling
,”
J. R. Soc. Interface
,
13
(
124
), p.
20160620
.10.1098/rsif.2016.0620
42.
Meek
,
K. M.
, and
Newton
,
R. H.
,
1999
, “
Organization of Collagen Fibrils in the Corneal Stroma in Relation to Mechanical Properties and Surgical Practice
,”
J. Refractive Surg.
,
15
(
6
), pp.
695
699
.10.3928/1081-597X-19991101-18
43.
Riveros
,
F.
,
Chandra
,
S.
,
Finol
,
E. A.
,
Gasser
,
T. C.
, and
Rodriguez
,
J. F.
,
2013
, “
A Pull-Back Algorithm to Determine the Unloaded Vascular Geometry in Anisotropic Hyperelastic Aaa Passive Mechanics
,”
Ann. Biomed. Eng.
,
41
(
4
), pp.
694
708
.10.1007/s10439-012-0712-3
44.
Ariza-Gracia
,
M. Á.
,
Zurita
,
J.
,
Piñero
,
D. P.
,
Calvo
,
B.
, and
Rodríguez-Matas
,
J. F.
,
2016
, “
Automatized Patient-Specific Methodology for Numerical Determination of Biomechanical Corneal Response
,”
Ann. Biomed. Eng.
,
44
(
5
), pp.
1753
1772
.10.1007/s10439-015-1426-0
45.
Hatami-Marbini
,
H.
,
2014
, “
Hydration Dependent Viscoelastic Tensile Behavior of Cornea
,”
Ann. Biomed. Eng.
,
42
(
8
), pp.
1740
1748
.10.1007/s10439-014-0996-6
46.
Hatami-Marbini
,
H.
, and
Maulik
,
R.
,
2016
, “
A Biphasic Transversely Isotropic Poroviscoelastic Model for the Unconfined Compression of Hydrated Soft Tissue
,”
ASME J. Biomech. Eng.
,
138
(
3
), p.
031003
.10.1115/1.4032059
47.
Gurtin
,
M. E.
,
Fried
,
E.
, and
Anand
,
L.
,
2010
,
The Mechanics and Thermodynamics of Continua
,
Cambridge University Press
,
Cambridge, UK
.