The inclusion of the epithelium in numerical models of the human cornea

Abstract

We present a patient-specific finite element model of the human cornea that accounts for the presence of the epithelium. The thin anterior layer that protects the cornea from the external actions has a scant relevance from the mechanical point of view, and it has been neglected in most numerical models of the cornea, which assign to the entire cornea the mechanical properties of the stroma. Yet, modern corneal topographers capture the geometry of the epithelium, which can be naturally included into a patient-specific solid model of the cornea, treated as a multi-layer solid. For numerical applications, the presence of a thin layer on the anterior cornea requires a finer discretization and the definition of two constitutive models (including the corresponding properties) for stroma and epithelium. In this study, we want to assess the relevance of the inclusion of the epithelium in the model of the cornea, by analyzing the effects in terms of uncertainties of the mechanical properties, stress distribution across the thickness, and numerical discretization. We conclude that if the epithelium is modeled as stroma, the material properties should be reduced by 10%. While this choice represents a sufficiently good approximation for the simulation of in vivo mechanical tests, it might result into an under-estimation of the postoperative stress in the simulation of refractive surgery.

Appendix A: Constitutive model of the stroma

We model the stroma as a hyperelastic composite, made of an elastic matrix (made of proteoglycans) reinforced with two sets of dispersed collagen fibrils (Pandolfi and Vasta 2012). The uncertainty of the fibril orientation is described with a von Mises distribution (Pandolfi and Vasta 2012) about a main orientation ${\textbf{a}}_M$, with $M=1,2$. The strain energy density is assumed to decompose additively into volumetric, isotropic–isochoric, and anisotropic–isochoric parts in the form

$$\begin{aligned} \Psi = \Psi _\textrm{vol}(J) + \Psi _\textrm{iso}(\overline{I}_1, \overline{I}_2) + \Psi _\textrm{aniso}({\overline{I}^*_{4\, M}}, {\sigma ^2_{I_4\, M}}) \,, \end{aligned}$$

where ${\textbf{F}}= d{\textbf{x}}/d{\textbf{X}}$ is the deformation gradient, ${\textbf{x}}$ are the current coordinates and ${\textbf{X}}$ the reference coordinates, and $J = \det {\textbf{F}}$ is the Jacobian determinant. ${\overline{\textbf{C}}}={\overline{\textbf{F}}}^T {\overline{\textbf{F}}}=J^{-2/3}{\textbf{F}}^T{\textbf{F}}$ is the isochoric Cauchy–Green deformation tensor, and ${\overline{I}_1}$ and ${\overline{I}_2}$ are the first and the second invariants of ${\overline{\textbf{C}}}$,

$$\begin{aligned} {\overline{I}_1}= \textrm{tr} {\overline{\textbf{C}}}\,, \qquad {\overline{I}_2}= 1/2\left[ (\textrm{tr}{\overline{\textbf{C}}})^2 - \textrm{tr}{\overline{\textbf{C}}}^2)\right] \,, \end{aligned}$$

where $\mathrm tr (.)$ denotes the trace operator. The average pseudo-invariant ${\overline{I}^*_{4\, M}}$ is defined as

$$\begin{aligned} {\overline{I}^*_{4\, M}}= \left[ \kappa _M {\textbf{I}}+ (1 - 3\kappa _M) {\textbf{a}}_M \otimes {\textbf{a}}_M \right] : {\overline{\textbf{C}}}\end{aligned}$$

where $\kappa _M$ is a dispersion parameter, $\otimes$ is the tensor product, and ( : ) is the double contraction product. The variance operator ${\sigma ^2_{I_4\, M}}$ is defined as

$$\begin{aligned} {\sigma ^2_{I_4\, M}}={\overline{\textbf{C}}}:{\textbf{a}}_M\otimes {\textbf{a}}_M\otimes {\textbf{a}}_M\otimes {\textbf{a}}_M:{\overline{\textbf{C}}}- {\overline{I}^*_{4\, M}}\end{aligned}$$

The mathematical form of the energies are

$$\begin{aligned} \Psi _\textrm{vol}= & {} \dfrac{1}{4} K \left( J^2 - 1 - 2 \log J \right) \,, \\ \Psi _\textrm{iso}= & {} \dfrac{1}{2} \left[ \mu _1 \left( {\overline{I}_1}-3 \right) + \mu _2 \left( {\overline{I}_2}-3 \right) \right] \,, \quad \mu _1 + \mu _2 = \mu \,, \\ \Psi _\textrm{aniso}= & {} \sum _{M=1}^2\dfrac{k_{1\,M}}{2~k_{2\,M}} \left[ \exp D^*\left( {\overline{I}^*_{4\, M}}\right) - 1 \right] \left( 1 + K^* {\sigma ^2_{I_4\, M}}\right) \,, \end{aligned}$$

where K is the bulk modulus, $\mu = \mu _1 + \mu _2$ is the shear modulus of the soft isotropic matrix, while $k_{1\,M}$ (stiffness-like parameter) and $k_{2\,M}$ (dimensionless rigidity parameters) control the mechanical response of the reinforcing fibers at low and high strains, respectively. The coefficient $D^*({\overline{I}^*_{4\, M}})$ reads

$$\begin{aligned} D^*\left( {\overline{I}^*_{4\, M}}\right) = k_{2\,M}\left( {\overline{I}^*_{4\, M}}- 1 \right) ^2 \,, \end{aligned}$$

and the coefficient $K^*$

$$\begin{aligned} k^* = k_{2\,M} \left[ 1 + 2D^*\left( {\overline{I}^*_{4\, M}}\right) \right] \,. \end{aligned}$$

The model parameters are seven: five with the dimension of a stiffness (shear elastic moduli), i.e., K, $\mu _1$, $\mu _2$, $k_{1\,1}$, $k_{1\,2}$ , and two dimensionless rigidity coefficients $k_{2\,1}$, $k_{2\, 2}$. The interested reader is referred to the original work (Pandolfi and Vasta 2012).