Computer Methods in Biomechanics and Biomedical

Engineering

ISSN: 1025-5842 (Print) 1476-8259 (Online) Journal homepage: http://www.tandfonline.com/loi/gcmb20

Sensitivity of corneal biomechanical and optical

behavior to material parameters using design of
experiments method

Mengchen Xu, Amy L. Lerner, Paul D. Funkenbusch, Ashutosh Richhariya &

Geunyoung Yoon

To cite this article: Mengchen Xu, Amy L. Lerner, Paul D. Funkenbusch, Ashutosh Richhariya
& Geunyoung Yoon (2018) Sensitivity of corneal biomechanical and optical behavior to material
parameters using design of experiments method, Computer Methods in Biomechanics and
Biomedical Engineering, 21:3, 287-296, DOI: 10.1080/10255842.2018.1447104

To link to this article: https://doi.org/10.1080/10255842.2018.1447104

Published online: 30 Mar 2018.

Submit your article to this journal

Article views: 25

View related articles

View Crossmark data

Full Terms & Conditions of access and use can be found at

http://www.tandfonline.com/action/journalInformation?journalCode=gcmb20
Computer Methods in Biomechanics and Biomedical Engineering, 2018
VOL. 21, NO. 3, 287–296
https://doi.org/10.1080/10255842.2018.1447104

Sensitivity of corneal biomechanical and optical behavior to material parameters

using design of experiments method
Mengchen Xua , Amy L. Lernera,b, Paul D. Funkenbuscha,c, Ashutosh Richhariyad and Geunyoung Yoonb,e
a
Department of Mechanical Engineering, University of Rochester, Rochester, NY, USA; bDepartment of Biomedical Engineering, University of
Rochester, Rochester, NY, USA; cMaterials Science Program, University of Rochester, Rochester, NY, USA; dL V Prasad Eye Institute, Kallam Anji
Reddy Campus, Hyderabad, India; eCenter of Visual Science, Flaum Eye Institute, The Institute of Optics, University of Rochester, Rochester, NY,
USA

ABSTRACT ARTICLE HISTORY

The optical performance of the human cornea under intraocular pressure (IOP) is the result of Received 14 February 2017
complex material properties and their interactions. The measurement of the numerous material Accepted 26 February 2018
parameters that define this material behavior may be key in the refinement of patient-specific KEYWORDS
models. The goal of this study was to investigate the relative contribution of these parameters to the Cornea; finite element
biomechanical and optical responses of human cornea predicted by a widely accepted anisotropic analysis; sensitivity analysis;
hyperelastic finite element model, with regional variations in the alignment of fibers. Design of biomechanics; optics
experiments methods were used to quantify the relative importance of material properties including
matrix stiffness, fiber stiffness, fiber nonlinearity and fiber dispersion under physiological IOP. Our
sensitivity results showed that corneal apical displacement was influenced nearly evenly by matrix
stiffness, fiber stiffness and nonlinearity. However, the variations in corneal optical aberrations
(refractive power and spherical aberration) were primarily dependent on the value of the matrix
stiffness. The optical aberrations predicted by variations in this material parameter were sufficiently
large to predict clinically important changes in retinal image quality. Therefore, well-characterized
individual variations in matrix stiffness could be critical in cornea modeling in order to reliably
predict optical behavior under different IOPs or after corneal surgery.

1. Introduction aberrations including both lower and higher orders (Dupps

and Roberts 2001; Roberts 2002). Under IOP, alterations
The main physiological function of the cornea is to trans-
in the corneal geometry and mechanical structure lead
mit and refract light into the eye to form an image on the
to changes in surface topography and as a consequence,
retina. Mechanically, the distribution of collagen fibrils
changes in optics. Corneal biomechanics and optics are
in the stroma reveals a laminar organization, in which
highly interrelated and the geometrical change associated
interactions between collagen fibrils and ground matrix
with material properties may be the key factor connect-
result in complex material properties such as inhomoge-
ing them. Finite element (FE) biomechanical modeling,
neity, anisotropy, viscoelasticity, and near incompressibil-
coupled with optical analysis, has been used as a powerful
ity. Optically, the anterior and posterior corneal surface
tool to investigate this interaction between corneal biome-
topographies, maintained by IOP, have a significant impact
chanics and optics (Pinsky et al. 2005; Alastrue et al. 2006;
on corneal aberrations, including both lower (defocus and
Pandolfi and Holzapfel 2008). This knowledge makes it
astigmatism) and higher (e.g. spherical aberration) order
possible to predict postoperative surgical outcome and can
aberrations, which degrade optical performance of the
be taken into account in surgical planning to minimize the
eye. Moreover, increased higher-order aberrations remain
biomechanically induced aberrations, further improving
one of the main limitations in laser refractive surgery, such
the patient’s visual performance.
as photorefractive keratectomy (PRK) and laser in situ
Recent advances in non-invasive in vivo corneal imag-
keratomileusis (LASIK). Corneal biomechanical response
ing techniques (Hashemi and Mehravaran 2007; Garcia
to the tissue removal in laser refractive surgery has been
and Rosen 2008) can precisely characterize corneal
proposed as one of the main factors that improve our
geometry and allow for development of patient-specific
understanding of the unexpected induction of optical

CONTACT Geunyoung Yoon gyoon@ur.rochester.edu

© 2018 Informa UK Limited, trading as Taylor & Francis Group
288  M. XU ET AL.

geometrical models (Simonini and Pandolfi 2015). criterion to exclude the contribution of collagen fibers
Although a patient-specific material model is an impor- under compression and remove a blunt error that leads
tant goal, one of the challenges has been to quantify indi- to a stress discontinuity. The modified strain energy func-
vidual material properties of the cornea. Inverse analysis tion that describes the anisotropic hyperelastic material
combined with experimental measurements has been behavior of the cornea was shown in Equations (1)–(4):
explored to estimate material parameters (Roy and Dupps ( ( el )2 )
2011; Roy et al. 2015). However, the solution of such iter- ) 1 J −1
𝜓 = C10 Ī1 − 3 + − lnJ el
(
ative procedures based on just apical displacement or one D 2
meridional profile of the cornea may not be unique or N (1)
k1 ∑ { ⟨ 2 ⟩
stable (Kabanikhin 2008), which could potentially gener- exp k2 Ē 𝛼 − 1
}
+
ate different optical results. Recent technologies such as 2k2 𝛼=1
Brillouin Optical Microscopy (Scarcelli et al. 2012) and
shear wave speed imaging using ultrasound (Nguyen et def (
Ē 𝛼 = 𝜅 Ī1 − 3 + (1 − 3𝜅) Ī4(𝛼𝛼) − 1 , (2)
) ( )
al. 2012) have demonstrated the potential to quantify the
instantaneous modulus of the cornea at a certain IOP in
vivo. However, these techniques cannot fully characterize 1𝜋
𝜅= ∫ 𝜌(𝜃)sin3 𝜃d𝜃, (3)
the biomechanical behavior of corneal stroma, such as 40
the stress stiffening effect or non-linearity observed ex
vivo. Meanwhile, each method involves a potential risk or 1 1 𝜋
𝜌(𝜃) = exp (bcos2𝜃), I = ∫ exp (bcos2𝜃)d𝜃, (4)
financial burden, so clinical studies may need to prioritize 2𝜋I 𝜋0
measurement of the most critical model inputs.
where C10 and k1 are stress-like parameters referring to
For the purpose of reducing the complexity in quan-
matrix stiffness and fiber stiffness, respectively; k2 is a
tifying individual material properties, it is essential to
dimensionless parameter referring to fiber nonlinearity
understand the relative significance of numerous mate-
and D represents the inverse of the volumetric modulus,
rial parameters in both corneal biomechanics and optics.
which was set to be negligibly small (1 × 10 kPa ~ kPa−1)
−6 −1

Therefore, the aim of this work is to investigate the sen-

for all the analyses considering corneal nearly incom-
sitivity of corneal biomechanical and optical behavior to
pressibility; Jel is the elastic volume ratio; Ī1 is the first
various material parameters using a design of experiments
invariant of the modified right Cauchy-green tensor;
approach. This was done based on a hyperelastic aniso-
Ē 𝛼 characterizes the deformation of the family of fibers
tropic FE model that incorporates the spatial variation in
with mean direction and Ī4𝛼𝛼 represents the square of the
fibril dispersion. The quantitative biomechanical outcome
stretch along the direction of the 𝛼-th family of fibers.
of apical displacement, optical outcomes such as spheri-
Two fiber families are considered (N = 2) and the degree
cal aberration and refractive power were analyzed over a
of their dispersion about preferential orientation is rep-
range of possible material parameters.
resented by the parameter 𝜅. It has been shown by X-ray
diffraction that fibril dispersion exhibits a spatial varia-
2. Materials and methods tion over the corneal surface (Boote et al. 2005; Pandolfi
and Holzapfel 2008) with fibers strongly aligned along the
2.1. Finite element model of cornea
nasal-temporal (NT) and superior-inferior (SI) meridians
A three-dimensional (3-D) corneal model was gener- and more dispersed in the four lobes (the four quadrant
ated using FEM software (ABAQUS, Dassault Systemes regions of corneal surface). This spatial variation follows
Simulia Corp.). To focus on the relative influence of mate- a transversely isotropic and π-periodic, normalized, von
rial parameters, we used a simplified axially symmetric Mises distribution (Equation 4), with the concentration
average corneal geometry (0.572 mm central corneal parameter b(𝜌, 𝜃) as a function of the distance from the
thickness, 0.779 mm peripheral corneal thickness and optical axis, 𝜌, and the angle of any meridian with the
7.97 mm radius of curvature of anterior surface) from positive horizontal axis, 𝜃. The distribution of b(𝜌, 𝜃)
previous measurements (Elsheikh et al. 2007). To describe was developed based on previous models (Pandolfi and
the microstructure of the cornea in a more computation- Holzapfel 2008). High values of b indicate highly oriented
ally economical way, a modified Gasser-Holzapfel-Ogden’s sets of fibers while low values of b define randomly ori-
strain energy function from the original GOH model ented fibers. To simplify the continuous variable b, we
(Gasser et al. 2006) with pre-integrated fiber distribution divided the cornea into different regions including five
was adopted in the commercial FE code. The modified sections in each lobe, NT and SI sections, a limbus section
GOH model used a different tension-compression switch and a transition-limbus section (Figure 1). Each section
 COMPUTER METHODS IN BIOMECHANICS AND BIOMEDICAL ENGINEERING  289

has an average value of b through the depth, established 2.2. Stress-free corneal geometry
based on 𝜌 and 𝜃 of all the nodes included in the respec-
Corneal topographies obtained from in vivo imaging
tive section. Fibers are mostly aligned in the NT, SI and
techniques are always under physiological IOP. Early
limbus sections. The degree of fibril dispersion increased
contributions of FE studies highlighted the importance
in each of the four lobe regions and reached maximum in
of achieving the pre-existing stress state or stress-free con-
the central section in each lobe.
figuration of the cornea to accurate numerical modeling
A convergence analysis of the finite element mesh was
of corneal biomechanics (Pinsky et al. 2005; Pandolfi and
performed based on the simulation of a general inflation
Manganiello 2006). In our study, the stress-free configu-
test with physiological IOP 15 mmHg. The number of
ration of the average corneal geometry under physiolog-
elements through the corneal thickness varied from 5 to
ical IOP of 15 mmHg was identified through the iterative
17 elements (5, 7, 11 and 17), corresponding to a total
algorithm developed based on Elsheikh’s work (Elsheikh
number of nodes ranging from 33,264 to 293,496. The
et al. 2013) for its advantages in implementation with
mesh was considered to be converged when the relative
our FE software and the high efficiency in iterations for
changes in all three measures, apical displacement, SA
ophthalmic applications. The algorithm incorporated
and refractive power of the cornea were smaller than
the consistent mapping of collagen fibril orientation and
assigned thresholds. A threshold of 0.005 mm for api-
dispersion with an initial set of averaged material coef-
cal displacement was selected based on reported exper-
ficients (C10 = 24.83 kPa, k1 = 91.55 kPa, k2 = 785.68,
imental data of apical displacement (Elsheikh et al.
b = 0.875–2.9) (Kok et al. 2014). Both corneal geometry
2007) and thresholds of 0.1 μm and 0.25D were used
and material properties have a significant impact on its
for SA and refractive power considering clinical signif-
biomechanical and optical behaviors, however, the geom-
icance, respectively. Based on the convergence analysis,
etry of the cornea can be measured precisely with several
the FE model was generated with 17 elements through
existing imaging techniques (Hashemi and Mehravaran
the thickness, resulting in a total of 293,496 nodes. To
2007; Garcia and Rosen 2008). To focus on the sensitivity
represent the much greater stiffness of the limbus, the
analysis of corneal material parameters, the same average
peripheral edge of the model were fixed with a simplified
stress-free corneal geometry was used as the starting point
zero displacement boundary condition (Alastrué 2005).
for the geometry in all subsequent simulations. As confir-
For the loading condition, uniform IOPs were applied
mation, apical displacement and optical aberrations were
to the posterior surface of the cornea, ranging from 5
analyzed for the stress-free model with average material
to 30 mm Hg.
properties to compare with previous inflation experiments
and normal optical aberrations values (Wang et al. 2003;
Elsheikh et al. 2007).

2.3. Optical analysis

The 3-D anterior and posterior corneal surface geome-
try data were exported from each finite element simu-
lation at three levels of IOP (10, 15 and 20 mmHg) to
compute refractive power and SA for the central 6 mm
corneal diameter. A custom-developed Matlab (2013, The
MathWorks, United States) program was used to charac-
terize the corneal surface map from simulation directly
using Zernike circle polynomials. An optical ray-tracing
program (Zemax, Radiant Zemax, LLC) was used to cal-
culate individual Zernike coefficients representing differ-
ent types of total wavefront aberrations induced by both
anterior and posterior corneal surfaces.
Figure 1. Fiber organization in the finite element model of cornea
with five different sections corresponding to different variations
in the concentration parameter b: 1- nasal-temporal (NT) and
2.4. Sensitivity study
superior-inferior (SI) section; 2- limbus section; 3- transition-
limbus section; 4- five sections in each lobe. Four material parameters were investigated, including
Notes: High values of b define highly oriented set of fibers. Low values of
b define dispersed set of fibers. From section 1 to section 4, the values of b matrix stiffness (C10), fiber stiffness (k1), fiber nonlinear-
decreased gradually. ity (k2) and the concentration parameter (b). A two-level
290  M. XU ET AL.

16 treatment condition full factorial array, which can 3. Results

estimate all factors and interaction effects, was used to
3.1. Stress-free corneal geometry
test the importance of these factors and their interactions
on the biomechanical and optical behavior of the cornea. To evaluate the estimation of the stress-free geometry, the
Variability in parameter values is incorporated into the stress-free state was pressurized again to compare with
sensitivity design by using two ‘levels’ (−1 and +1) for the target corneal surface topography under 15 mmHg
each parameter and assigning specific parameter values IOP. The difference between the target and pressurized
to each. The influence of each parameter is evaluated by corneal surface topography was negligibly small, with a
running the model for different treatment conditions, maximum RMS value of 3.0E-04 mm and an average RMS
each reflecting a specific combination of parameter levels. value of 2.9E-05 mm, demonstrating the reliability of the
The high (+1) and low (−1) levels used for C10, k1 and k2 iterative procedure.
were selected from previous experimental data and inverse
analyses for human corneas (Elsheikh et al. 2007; Kok et
3.2. Validation of corneal behavior with average
al. 2014). For parameter b, since there was considerable
material properties
variation between different corneas in the proportion of
collagen fibrils oriented in the main orientations (Boote et Both the magnitudes and the nonlinearity of the pre-
al. 2005), the assumption was made that in the distribution dicted apical displacement with average material prop-
map (Figure 1): b (−1) = 1.55–2.9 and b (+1) = 0.2–2.9. erties compared well with inflation test data of human
By setting b in this way, the fibrils in NT, SI and limbus corneas (Figure 2) (Elsheikh et al. 2007). In addition, the
sections remained highly aligned between the two levels, values of SA and refractive power from the simulation
however the dispersion degree in the other sections varied, were 0.06, 0.33 and 0.57 μm and 47.8, 48.1 and 47.6 D
especially in the lobed region. The levels selected for each for 10, 15 and 20 mmHg IOP, respectively. These are all
factor are summarized in Table 1. reasonable aberration values for normal human corneas
The sensitivity analysis was performed based on (Wang et al. 2003).
changes in biomechanical and optical responses at
three physiological IOPs (10, 15 and 20 mmHg). For the
3.3. Sensitivities of biomechanical and optical
model, the stress-free geometry is independent of the
behavior
mechanical properties and, therefore, provides a con-
venient reference against which the response may be The ANOM results showed the magnitude of the effects
measured. The apical displacement in mm of the cornea (the difference in the response attributable to the change
is defined as the difference between the apical position in the factor between its two levels) for all of the factors.
with the IOP applied and the stress-free, reference, state. An example of ANOM results for 15 mmHg is shown
Similarly, the change in refractive power in diopters (D) in Figure 3. For apical displacement, the difference in
and the SA in μm are determined relative to the stress- response observed for k1 between the +1 and −1 level
free reference state because the stress-free cornea is not was the largest, with a magnitude of 0.11 mm, while the
optically perfect, i.e. has some aberrations. Analysis of
means (ANOM) was used to quantify the magnitude of
the effect and analysis of variance (ANOVA) was used to
determine the variance (sum of squares or SS) attributa-
ble to each factor and interaction. The percentage of total
sum of squares (%TSS) for each factor can then be used
to compare relative contributions to the total variance
of the response.

Table 1. Values of the two levels for material parameters investi-

gated in sensitivity analysis.
Level
Parameters −1 +1
Matrix stiffness (C10: kPa) 2.75 46.91 Figure 2. Comparison of apical displacements result from model
Fiber stiffness (k1: kPa) 21.7 161.4 with average material properties and inflation tests data of
Fiber nonlinearity (k2) 291.28 1280.07 different corneas with age range 50–95 years old from Elsheikh’s
Concentration parameter (b) 1.55–2.9 0.2–2.9
experiments (Elsheikh et al. 2007).
 COMPUTER METHODS IN BIOMECHANICS AND BIOMEDICAL ENGINEERING  291

Figure 3. (a) Graphical presentation of ANOM results for variations in apical displacement at 15 mmHg IOP. Each line reflects the variation
from the + 1 to −1 level. Vertical axis depicts the corresponding average response for each factor and interaction with +1 and −1 level.
The slope of each line represents how much difference in apical displacement is caused by varying factors from +1 to −1 level. (b)
Graphical presentation of ANOM results for variations in spherical aberration change at 15 mmHg IOP. Each line reflects the variation
from the +1 to −1 level. Vertical axis depicts the corresponding average response for each factor and interaction with +1 and −1 level.
The slope of each line represents how much difference in spherical aberration change is caused by varying factors from +1 to −1 level.
(c) Graphical presentation of ANOM results for variations in refractive power change at 15 mmHg IOP. Each line reflects the variation from
the +1 to −1 level. Vertical axis depicts the corresponding average response for each factor and interaction with +1 and −1 level. The
slope of each line represents how much difference in refractive power change is caused by varying factors from +1 to −1 level.

differences observed for C10 and k2 were slightly smaller. for three pressures, respectively. The effects of k1 and k2
In contrast, the largest differences in optical responses increased with IOP while the effect of C10 decreased. For
were observed with changes in C10, with magnitudes of the optical behaviors, the matrix stiffness C10 always had
1.43 μm and 5.19 D for change in SA and refractive power, the predominant impact, 70%TSS for change in SA and
respectively. These aberration values are sufficiently large 79%TSS for change in refractive power in average for three
to cause a clinically noticeable degradation in retinal pressures. In addition to C10, the concentration parameter
image quality. b (15% TSS) and interaction term C10 × b (11%TSS) also
The ANOVA for the biomechanical and optical had considerable impact on change in SA. However, for
results judged by %TSS for the three IOPs are shown in refractive power, all of the other factors made only minor
Figure 4(a)–(c). Some of the interaction terms had minor contributions. The significance of C10 decreased with IOP
effects (accounted for less than 1%TSS). For example, for changes in SA and, conversely, increased for changes
the four factor interaction (C10 × k1 × k2 × b) was con- in refractive power.
sistently negligible and only one three factor interaction
(C10 × k1 × k2) was not negligible. The two factor interac-
4. Discussion
tions that were important were generally consistent with
the most important factors. For apical displacement, the Uncertainty about individual variations in corneal mate-
results showed that C10, k1 and k2 all contributed to the rial properties may limit our ability to predict corneal
variance, accounting for an average 22, 31 and 22%TSS biomechanical and optical response to IOPs and tissue
292  M. XU ET AL.

Figure 4. (a) Percentage of total sum of squares of identified important factors for different pressures (15–25 mmHg) to the variation
in apical displacement. (b) Percentage of total sum of squares of identified important factors for different pressures (15–25 mmHg) to
the variation in spherical aberration change. (c) Percentage of total sum of squares of identified important factors for different pressures
(15–25 mmHg) to the variation in refractive power change.

removal in refractive surgery. The present study examined (Holzapfel et al. 2004; Holzapfel and Ogden 2010). From
the relative importance of material parameters on corneal the numerical point of view, the compressed fibers need
apical displacement and optical aberrations. The matrix to be excluded in the strain energy function to describe
stiffness was identified as the most significant material the anisotropic hyperelastic behavior of the cornea. One
parameter contributing to both biomechanical and optical of the controversial issues of the GOH model is how
behaviors of the cornea under different IOPs. Interactions to account for the fiber distribution stiffness when the
between the material parameters were found to have a rel- stretch in the main fiber direction is less (or equal) than
atively small impact. This suggests that a less comprehen- 1, especially in the case of dispersed fibers (Cortes et
sive approach, concentrating only on the factors identified al. 2010; Pandolfi and Vasta 2012). This singularity has
as important, might be adequate in some future theoretical been handled through a modified tension-compres-
and experimental studies. sion switch criterion in the current FE code in Abaqus
The anisotropic hyperelastic corneal model with (Dassault Systemes Simulia Corp) to exclude the con-
a modified Gasser-Holzapfel-Ogden’s strain energy tribution of collagen fibers under compression. This
function (Gasser et al. 2006; Pandolfi and Holzapfel approach removes a blunt error that leads to a stress
2008) implemented in commercial FE code Abaqus discontinuity in the original GOH model (Gasser et al.
was selected for this study based on its advantages in 2006). Although the treatment of the switch for dis-
lower computational cost compared to the Angular persed fibers in Abaqus could still be improved using
Integration models (Sacks 2003; Driessen et al. 2005; more sophisticated approaches, such as a new pre-in-
Asejczyk-Widlicka et al. 2009). In soft fibrous material, tegrated proposal presented in Latorre’s work (Latorre
such as cornea, collagen fibers are the principal compo- and Montans 2016), we think that the exclusion of the
nent that provides the mechanical stiffness and strength. contribution of collagen fibers under compression is
Because of the wavy (crimped) structure, it is assumed still an acceptable approximation for our current goal
that collagen fibers are not able to support compres- of analyzing the overall trend of corneal responses to
sion since the fibers would buckle under compression variations in physiological IOP loading.
 COMPUTER METHODS IN BIOMECHANICS AND BIOMEDICAL ENGINEERING  293

4.1. Contribution of matrix and fibrils consistent with previous findings (Pandolfi and Holzapfel
2008). The current results further quantified its impact on
Three factors (matrix stiffness C10, fiber stiffness k1and
one of the most common corneal high order aberrations,
fiber nonlinearity k2) all contributed to corneal apical dis-
SA. Because of the increased degree of fibril dispersion,
placement, but their contributions changed with IOP over
more surface deformation occurs in the four lobe sections
the range examined. Both k1 and k2 became more impor-
of the cornea compared to the NT/ SI sections. This ele-
tant than C10 with increasing pressure (Figure 4(a)), which
vation difference directly contributed to the resulting SA
supports the understanding of the nonlinear stress-strain
value, as seen by comparison of the 0° and 45° meridian
behavior of the cornea (Ruberti et al. 2011). Under lower
profiles (Figure 5). After surface fitting, the softest cornea
pressure, variations in the matrix contributed to apical dis-
(Figure 5(a)) shows the largest elevation difference and
placement more, resulting in an initial matrix dominated
was more curved towards the edge, which resulted in the
linear behavior. As the pressure increased, fibers in tension
largest positive SA value of 2.32 μm among all conditions.
gradually produced more resistance, corresponding to a
In contrast, for the stiffest cornea (Figure 5(b)), the eleva-
fibril-dominated highly nonlinear behavior.
tion difference was smallest. The overall geometry of the
Optical behaviors of the cornea were most affected by
peripheral cornea was flatter which resulted in the largest
C10. Relative contributions of other factors and interac-
negative SA of 0.13 μm. However, the results were based
tions were much lower. Among these, the concentration
on one sample assumption of the variation in b, more
parameter b was a factor that could not be ignored for
realistic analysis based on X-ray data would be necessary
analyzing SA. It has been known from X-ray diffraction
in future studies.
studies that from central to peripheral cornea, collagen
fibers covert from a strongly aligned state along NT, SI
meridians to a more dispersed state in the transition zone 4.2. Indications for individualized corneal modeling
(the four lobes section) and again a strongly aligned state
in the limbus area (Boote et al. 2005; Abahussin et al. Our sensitivity results suggest prioritizing the assessment
2009). However, the dispersion degree of collagen fibers of C10 and b for individualized material models to bet-
in the transition zone may vary among individual corneas ter understand corneal optical behavior, thus potentially
and it is important to include this in a sensitivity analy- reducing the complexity of both the inverse analysis and
sis. It has been shown that there is a clearly considerable direct measurements. Several newer technologies provide
variation in the proportion of fibrils oriented within 45° the potential to characterize the shear or elastic modulus
sectors of the NT, SI meridians between different healthy of individual cornea, such as shear wave speed imaging
human cornea specimens (Boote et al. 2005). In addition, using ultrasound (Nguyen et al. 2014), optical coher-
the spatial distribution of fibril dispersion might also be ence elastography (Ford et al. 2011) and Brillouin optical
different between healthy and diseased corneas, such as microscopy (Scarcelli et al. 2012). However, their correla-
keratoconus cornea (Meek and Quantock 2001). To inves- tion with the real matrix stiffness still requires additional
tigate the potential impact of the variation in fibril disper- investigation. Further quantification of the fibril concen-
sion degree, we maintained the higher limit of b in areas tration parameter b is also essential to understand optical
where fibers are known to be highly aligned, while incor- changes. X-ray diffraction has been used to investigate
porating the variability in the alignment of other regions collagen fibril structure and its relation to the mechanical
by adjusting the lower limit on b from 0.2 (nearly isotropic properties of the cornea (Boote et al. 2005). However, the
dispersion) to an assumed value of 1.55 (partial of the variation in fibril dispersion among different corneas has
fibers are strongly oriented and partial of the fibers are dis- not been fully characterized.
persed). In this fashion, the analysis can assess the mod-
el’s sensitivity to the variability in the dispersion degree 4.3. Comparison with previous studies and future
(Section 4), and also retaining an appropriate pattern with improvements
respect to the variation of dispersion with location on the
cornea. Our pilot study showed that fibril dispersion was Previous sensitivity study (Ariza-Gracia et al. 2016) using
the dominant parameter in the model with uniform fibril the same material model identified fiber stiffness k1 as the
dispersion (Alastrué 2005; Ariza-Gracia et al. 2016). By most influential material parameter on corneal apical dis-
incorporating the anisotropic fibril dispersion revealed by placement during air applanation. In contrast, the current
X-ray studies (Boote et al. 2005), the significance of differ- study found that a combination of C10, k1 and k2 contrib-
ent material parameters was affected. Therefore, the spatial uted to the apical displacement in inflation test. The differ-
variation of the collagen fibril dispersion is important to ence in material sensitivities may be a result of differences
consider in investigating corneal behaviors numerically, in model descriptions, including the loading condition
294  M. XU ET AL.

Figure 5. (a) The cross-sectional surface profile and displacement map (15 mmHg IOP) along 0° and 45° meridians with respect to the
radial distance from the center of cornea for extreme combination of material parameters: C10 (−1), k1 (−1), k2 (−1), b(+1), corresponding
to the softest material set. Fibrils are most dispersed in the lobe areas, which lead to the largest elevation difference the lobes and NT
band. (b) The cross-sectional surface profile and displacement map (15 mmHg IOP) along 0° and 45° meridians with respect to the radial
distance from the center of cornea for extreme combination of material parameters: C10 (+1), k1 (+1), k2 (+1), b(−1), corresponding to
the stiffest material set. Dispersion degree of fibrils in the lobe areas was lowest and the elevation difference was smallest.

(air applanation vs. inflation), the fibril dispersion pattern current study of the cornea under physiological IOP, an
(uniform vs. regionally varied) and the selected magni- improved structure tensor (Cortes et al. 2010; Pandolfi
tudes of k1 (130 MPa vs. 91.55 kPa on average). Future and Vasta 2012) would not change the sensitivity results,
studies will need to consider the potential impact of these but would be needed for predicting accurate internal
differences. First, the relative role of different material stress of the cornea. The regional variation in fibril dis-
parameters may strongly depend on the specific loading persion appeared to be an acceptable approximation for
scenario. An independent sensitivity study may be neces- the current case, as it captured the overall trend in topog-
sary depending on the analytical purpose. Second, further raphy changes and resulting effects on optical behavior.
experimental characterization of fibril dispersion patterns However, for applications focusing on localized material
among human corneas is important to better justify model and optical behavior, a continuous variation would be nec-
choices. Finally, accurate quantification of the parame- essary. The current study also used a cornea-only model
ter range in a broad population might be necessary to with a simplified limbus boundary condition instead of
improve the sensitivity analysis results. a whole-eye model considering the cost of the iterative
Several other improvements are also suggested for stress-free algorithm. The cornea-only model can predict
future work. The current model used uniform material sufficiently accurate refractive power results (Alastrué
properties throughout the corneal depth. However, more 2005) and a similar trend of apical displacement (Roy
detailed description of depth variation in matrix stiffness and Dupps 2009), however, several other studies (Pandolfi
and fiber dispersion (Scarcelli et al. 2012; Abass et al. 2015) and Holzapfel 2008; Elsheikh et al. 2010) has pointed out
may be important, especially for refractive surgery stud- the impact of compliant boundary conditions of the cor-
ies. The material model selected in the current study has nea-only model on simulating the biomechanical response
drawbacks in stress estimation in several loading condi- of the cornea. To avoid the concern, the correlation of the
tions such as uniaxial, biaxial tension and shear. For the biomechanical and optical responses predicted by the FE
 COMPUTER METHODS IN BIOMECHANICS AND BIOMEDICAL ENGINEERING  295

models with fixed and compliant boundary conditions in Alastrue V, Calvo B, Pena E, Doblare M. 2006. Biomechanical
our sensitivity analysis was investigated, with a R2 value modeling of refractive corneal surgery. J Biomech Eng-T
of 0.929, 0.997 and 0.849 for the apical displacements, SA Asme. 128:150–160.
Ariza-Gracia MA, Zurita J, Pinero DP, Calvo B, Rodriguez-
and refractive power, respectively. The sensitivity trends Matas JF. 2016 May. Automatized patient-specific
should not be affected but, for predicting accurate val- methodology for numerical determination of biomechanical
ues of biomechanical and optical behaviors, a whole-eye corneal response. Ann Biomed Eng. 44:1753–1772.
model may still be needed. In addition, in modeling the Asejczyk-Widlicka M, Srodka W, Krzyzanowska-Berkowska P.
general response of the cornea to physiological IOP, we did 2009. The biomechanical modelling of the refractive surgery.
Optik. 120:923–933.
not account for its multiphasic characteristics. However,
Boote C, Dennis S, Huang YF, Quantock AJ, Meek KM. 2005
in cases where viscoelastic material behavior might be Jan. Lamellar orientation in human cornea in relation to
critical, our approach for sensitivity analysis may be useful mechanical properties. J Struct Biol. 149:1–6.
in developing an understanding of the relevant parameters Cheng X, Petsche SJ, Pinsky PM. 2015 Feb. A structural model
in those material models (Cheng et al. 2015). for the in vivo human cornea including collagen-swelling
interaction. J R Soc Interface R Soc. 6:12.
Cortes DH, Lake SP, Kadlowec JA, Soslowsky LJ, Elliott DM.
4.4. Conclusion 2010 Oct. Characterizing the mechanical contribution of
fiber angular distribution in connective tissue: comparison
The sensitivity analysis indicated that matrix stiffness was of two modeling approaches. Biomech Model Mechanobiol.
the most significant biomechanical parameter affecting 9:651–658.
corneal optical behavior. For understanding and predict- Driessen NJB, Bouten CVC, Baaijens FPT. 2005 Jun. A
structural constitutive model for collagenous cardiovascular
ing higher-order aberrations such as spherical aberration,
tissues incorporating the angular fiber distribution. J
fibril dispersion was also an important parameter, and Biomech Eng-T Asme. 127:494–503.
should be prioritized for further characterization. These Dupps WJ, Roberts C. 2001 Nov–Dec. Effect of acute
findings will be beneficial in developing individualized biomechanical changes on corneal curvature after
corneal models from either inverse analysis or experimen- photokeratectomy. J Refract Surg. 17:658–669.
tal measurements, and will advance the ability to predict Elsheikh A, Wang D, Pye D. 2007 Oct. Determination of the
modulus of elasticity of the human cornea. J Refract Surg.
surgery outcomes for individual patients. 23:808–818.
Elsheikh A, Geraghty B, Rama P, Campanelli M, Meek KM.
2010 Oct 6. Characterization of age- related variation in
Disclosure statement corneal biomechanical properties. J Roy Soc Interface.
No potential conflict of interest was reported by the authors. 7:1475–1485.
Elsheikh A, Whitford C, Hamarashid R, Kassem W, Joda A,
Buchler P. 2013 Feb. Stress free configuration of the human
Funding eye. Med Eng Phys. 35:211–216.
Ford MR, Dupps WJ Jr, Rollins AM, Roy AS, Hu Z. 2011
This work was supported by the Research to Prevent Blindness Jan–Feb. Method for optical coherence elastography of the
and the National Science Foundation [grant number CMMI- cornea. J Biomed Opt. 16:016005.
1100632]. Garcia JPS, Rosen RB. 2008 Nov–Dec. Anterior segment
imaging: optical coherence tomography versus ultrasound
biomicroscopy. Ophthal Surg Las Im. 39:476–484.
ORCID Gasser TC, Ogden RW, Holzapfel GA. 2006 Feb 22. Hyperelastic
Mengchen Xu http://orcid.org/0000-0003-4443-1809 modelling of arterial layers with distributed collagen fibre
orientations. J Roy Soc Interface. 3:15–35.
Hashemi H, Mehravaran S. 2007 Oct. Central corneal thickness
References measurement with Pentacam, Orbscan II, and ultrasound
devices before and after laser refractive surgery for myopia.
Abahussin M, Hayes S, Cartwright NEK, Kamma-Lorger J Cataract Refract Surg. 33:1701–1707.
CS, Khan Y, Marshall J, Meek KM. 2009 Nov. 3D collagen Holzapfel GA, Ogden RW. 2010 Jun 8. Constitutive modelling
orientation study of the human cornea using x-ray diffraction of arteries. P Roy Soc a-Math Phy. 466:1551–1596.
and femtosecond laser technology. Invest Ophthalmol Vis Holzapfel GA, Gasser TC, Ogden RW. 2004 Apr. Comparison
Sci. 50:5159–5164. of a multi-layer structural model for arterial walls with a
Abass A, Hayes S, White N, Sorensen T, Meek KM. 2015 Mar fung-type model, and issues of material stability. J Biomech
6. Transverse depth-dependent changes in corneal collagen Eng-T Asme. 126:264–275.
lamellar orientation and distribution. J Roy Soc Interface. Kabanikhin SI. 2008 Jul. Definitions and examples of inverse
12:20140717. and ill-posed problems. J Inverse Ill-Pose P. 16:317–357.
Alastrué V. 2005. Biomechanical modeling of refractive corneal Kok S, Botha N, Inglis HM. 2014 Dec. Calibrating corneal
surgery. J Biomech Eng. 128:150. material model parameters using only inflation data: an ill-
296  M. XU ET AL.

posed problem. Int J Num Methods Biomed Eng. 30:1460– Roberts C. 2002 Sep–Oct. Biomechanics of the cornea and
1475. wavefront-guided laser refractive surgery. J Refract Surg.
Latorre M, Montans FJ. 2016 Apr. On the tension-compression 18:S589–592.
switch of the gasser-ogden-holzapfel model: analysis and a Roy AS, Dupps WJ. 2009 Oct. Effects of altered corneal stiffness
new pre-integrated proposal. J Mech Behav Biomed Mater. on native and postoperative LASIK corneal biomechanical
57:175–189. behavior: a whole-eye finite element analysis. J Refract Surg.
Meek KM, Quantock AJ. 2001 Jan. The use of x-ray scattering 25:875–887.
techniques to determine corneal ultrastructure. Prog Retinal Roy AS, Dupps WJ. 2011 Jan. Patient-specific modeling of
Eye Res. 20:95–137. corneal refractive surgery outcomes and inverse estimation
Nguyen TM, Aubry JF, Touboul D, Fink M, Gennisson JL, of elastic property changes. J Biomech Eng-T Asme.
Bercoff J, Tanter M. 2012 Aug. Monitoring of cornea elastic 133:011002.
properties changes during UV-A/riboflavin-induced corneal Roy AS, Kurian M, Matalia H, Shetty R. 2015 Aug. Air-puff
collagen cross-linking using supersonic shear wave imaging: associated quantification of non-linear biomechanical
a pilot study. Invest Ophthalmol Vis Sci. 53:5948–5954. properties of the human cornea in vivo. J Mech Behav
Nguyen TM, Aubry JF, Fink M, Bercoff J, Tanter M. 2014 Biomed Mater. 48:173–182.
Nov. In vivo evidence of porcine cornea anisotropy using Ruberti JW, Roy AS, Roberts CJ. 2011. Corneal biomechanics
supersonic shear wave imaging. Invest Ophthalmol Vis Sci. and biomaterials. Ann Rev Biomed Eng. 13(13):269–295.
55:7545–7552. Sacks MS. 2003 Apr. Incorporation of experimentally-derived
Pandolfi A, Holzapfel GA. 2008 Dec. Three-dimensional fiber orientation into a structural constitutive model
modeling and computational analysis of the human cornea for planar collagenous tissues. J Biomech Eng-T Asme.
considering distributed collagen fibril orientations. J 125:280–287.
Biomech Eng-T Asme. 130:061006. Scarcelli G, Pineda R, Yun SH. 2012 Jan. Brillouin optical
Pandolfi A, Manganiello F. 2006 Nov. A model for the human microscopy for corneal biomechanics. Invest Ophthalmol
cornea: constitutive formulation and numerical analysis. Vis Sci. 53:185–190.
Biomech Model Mechanobiol. 5:237–246. Simonini I, Pandolfi A. 2015 Jun 22. Customized finite element
Pandolfi A, Vasta M. 2012 Jan. Fiber distributed hyperelastic modelling of the human cornea. PLoS One. 10:e0130426.
modeling of biological tissues. Mech Mater. 44:151–162. Wang L, Dai E, Koch DD, Nathoo A. 2003 Aug. Optical
Pinsky PM, van der Heide D, Chernyak D. 2005 Jan. aberrations of the human anterior cornea. J Cataract Refract
Computational modeling of mechanical anisotropy in the Surg. 29:1514–1521.
cornea and sclera. J Cataract Refract Surg. 31:136–145.