This was a simulation study using finite element modeling. The method for the finite element model generation is briefly described here from another recent study.
11 Three-dimensional (3-D) geometry of a patient cornea was created from Pentacam (OCULUS Optikgerate Gmbh). The Pentacam provided Cartesian coordinates of the anterior and posterior corneal surface, which were used for creating a 3-D volume. Epithelium thickness of the cornea was measured with RTVue (Optovue, Inc., Fremont, CA, USA). Finite element mesh was created with 8-noded linear hexahedral elements (TrueGrid; XYZ Scientific Applications, Inc., Livermore, CA, USA). A total of 3312 hexahedral elements were used to represent the corneal volume.
Figures 1A and
1B show a cross section of the preoperative and postoperative mesh (with flap/cap) of the central cornea, respectively. PRK finite element mesh did not have any flap/cap zone. An anisotropic, hyperelastic, fiber-dependent material model with material incompressibility was chosen.
1 The material model accounted for the orthogonal arrangement of fibers in the central cornea, depth dependency of angular direction of the interweaving fibers, and reorientation of the in-plane peripheral collagen fibers to circumferential direction.
11 The hyperelastic material model was represented by free energy density (
ψ):
\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}\tag{1}\psi \left( C \right) = {\psi _m}\left( {{I_1},{I_3}} \right) + {\psi _{f \mbox{-} plane}}\left( C \right) + {\psi _{f \mbox{-} cross}}\left( C \right)\end{equation}
where
I1 = C : 1 and
I3 = det[
C] were the 1st and 3rd invariants of the deformation tensor.
C was the right Cauchy-Green deformation tensor and det[
C] represented determinant of the tensor. The isotropic energy density of the matrix (
ψm) was described by:
\begin{equation}\tag{2}{\psi _m}\left( {{I_1},{I_3}} \right) = {C_1}\left( {{{\overline I }_1} - 3} \right) + {C_2}{\left( {{{\overline I }_1} - 3} \right)^2} + {{{{\left( {{I_3} - 1} \right)}^2}} \over {D_1}}\end{equation}
where
Display Formula\({\bar I_1} = I_3^{ - 1/3}{I_1}\) was the distortional part of
I1,
I3 was the determinant of deformation gradient tensor, and
C1s were the material constants.
D1 was the bulk modulus to enforce incompressibility and was assumed to be 10
−6. The fiber energy density (
ψf-plane) was described by:
\begin{equation}\tag{3}{\psi _{f \mbox{-} plane}}\left( C \right) = \mathop \smallint \limits_{ - \pi }^\pi {W_f}\left( {{\lambda _{f \mbox{-} plane}}} \right){D_{plane}}\left( \theta \right)d\theta \end{equation}
\begin{equation}\tag{4}{W_{f \mbox{-} plane}}\left( {{\lambda _{f \mbox{-}plane}}} \right) = {{k_1} \over {2{k_2}}}\exp \left( {{k_2}{{\left[ {{\lambda _{f \mbox{-}plane}} - 1} \right]}^2}} \right) - {{k_1} \over {2{k_2}}}\end{equation}
represented the energy density of in-plane lamellar collagen fibers with stretch,
λf-plane equal to
Display Formula\(\sqrt {A.CA} \), where
A = [cos
θ, sin
θ, 0]
T was the local direction vector of the fibers.
k1 and
k2 were the material constants.
Dplane(
θ) represented a weighted average of the energy density of the fiber families at each integration point of the element. It also represented the change in the preferred direction of the fibers from orthogonal in the central cornea to circumferential near the limbus.
\begin{equation}\tag{5}{\psi _{f \mbox{-} cross}}\left( C \right) = \mathop \smallint \limits_{ - \pi }^\pi {W_f}\left( {{\lambda _{f \mbox{-} cross}}} \right){D_{cross}}\left( \theta \right)d\theta \end{equation}
\begin{equation}\tag{6}{W_{f \mbox{-} cross}}\left( {{\lambda _{f \mbox{-} cross}}} \right) = {{k_1} \over {2{k_2}}}\exp \left( {{k_2}{{\left[ {{\lambda _{f \mbox{-} cross}} - 1} \right]}^2}} \right) - {{k_1} \over {2{k_2}}}\end{equation}
represented the energy density of crosslink fibers between the lamellae with stretch,
λf-cross equal to
Display Formula\(\sqrt {B.CB} \), where
B was the direction vector of the crosslink fibers.
B was determined by taking a cross-product of
A with the surface normal and then rotating it out-of-plane around
A by an angle
ξ. The angle
ξ was assumed to be a function of depth and was modeled as follows
11:
\begin{equation}\tag{7}\xi = {28.6^\circ }{{\exp \left( {3.19\left[ {1 - s} \right]} \right) - 1} \over {\exp \left( {3.19} \right) - 1}}\end{equation}
was the nondimensional local thickness. The angle
ξ was evaluated at each element centroid.
Dplane(
θ) was kept equal to
Dcross(
θ). From the above equations, the Cauchy stress was determined by:
\begin{equation}\tag{8}\sigma = {1 \over {\sqrt {I_3} }}F{{\partial \psi } \over {\partial C}}{F^T}\end{equation}
where
F was the deformation gradient tensor. The epithelium was modeled as an isotropic, hyperelastic, incompressible material [
Display Formula\(\psi \left( C \right) = {\psi _m}\left( {{I_1},{I_3}} \right)\) only] with
c1 = 5 kPa (Pascal) and
c2 = 0.0 kPa.