journal of the mechanical behavior of biomedical materials 48 (2015) 173–182

Available online at www.sciencedirect.com

www.elsevier.com/locate/jmbbm

Research Paper

Air-puff associated quantification of non-linear biomechanical properties of the human cornea in vivo Abhijit Sinha Roya,n, Mathew Kurianb, Himanshu Mataliac, Rohit Shettyc,d,1 a

Imaging, Biomechanics and Mathematical Modeling Solutions, Narayana Nethralaya, Bangalore, India Cornea and Cataract Services, Narayana Nethralaya, Bangalore, India c Corneal and Refractive Services, Narayana Nethralaya, Bangalore, India d Narayana Nethralaya, Bangalore, India b

art i cle i nfo

ab st rac t

Article history:

With the advent of newer techniques to correct refraction such as flapless laser procedure

Received 22 January 2015

and collagen crosslinking, in vivo estimation of corneal biomechanical properties has

Received in revised form

gained importance. In this study, a new 3-D patient specific inverse finite element method

5 April 2015

of estimating corneal biomechanical properties from air-puff applanation was developed.

Accepted 9 April 2015

The highlight of the model was inclusion of patient-specific corneal tomography, fiber

Available online 20 April 2015

dependent hyperelastic model, cross links between collagen lamellae and epithelium layer.

Keywords:

A lumped mass, spring and dashpot model was included to model the resistance to motion

Cornea

and deformation of the eye globe caused by air-puff applanation. 10 normal eyes of 10

Biomechanics

human subjects were used for the study. 3-D finite element models were constructed and

Ansiotropic

custom routines were scripted for performing the inverse calculations. The model for each

Collagen fiber

eye was perturbed to estimate the effect of measured intraocular pressure on the

Hyperelastic

estimated biomechanical variables. The study demonstrated that the inverse method was effective in quantification of material properties and was sensitive to intraocular pressure alterations. Specifically, in vivo fiber dependent hyperelastic biomechanical properties of human corneas were estimated for the first time. & 2015 Elsevier Ltd. All rights reserved.

1.

Introduction

There are several refractive procedures where biomechanics of the cornea plays an important role in outcomes of the treatment (Ruberti et al., 2011; Terai et al., 2012). For example, laser-assisted in situ keratomileusis (LASIK) reshapes the anterior stroma and causes a biomechanical change n

Corresponding author. Tel.: þ91 9740566833. E-mail address: [email protected] (A. Sinha Roy). 1 Vice-Chairman.

http://dx.doi.org/10.1016/j.jmbbm.2015.04.010 1751-6161/& 2015 Elsevier Ltd. All rights reserved.

(Roberts, 2002; Pedersen et al., 2014). Collagen crosslinking in keratoconus cornea strengthens the stroma biomechanically to slow the progression of the disease (Terai et al., 2012). While biomechanics plays a key role in several other corneal treatments, quantification of classical non-linear viscoelastic properties of the cornea became possible only recently (Ruberti et al., 2011; Piñero and Alcón, 2014). A recent review

174

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

summarized several studies that quantified the biomechanical properties ex vivo using either human or animal corneal tissue (Ruberti et al., 2011). Tensile modulus of a cornea can be determined ex vivo by either inflation testing of corneal buttons or by tensile stretching of rectangular cut corneal specimens (Wollensak et al., 2003; Ruberti et al., 2011). However, these methods cannot be replicated in vivo in routine clinical practice (Ruberti et al., 2011). Newer technologies may provide the in vivo quantification. Shear wave speed imaging using ultrasound can measure shear and Young's modulus of the corneal tissue in vivo (Nguyen et al., 2012; Torricelli et al., 2014; Touboul et al., 2014). The technique can resolve depth and radial variation in the modulus of the cornea (Nguyen et al., 2012; Torricelli et al., 2014; Touboul et al., 2014). However, it approximates tissue behavior as a linear, elastic material and requires contact with the tissue during measurement (Nguyen et al., 2012; Torricelli et al., 2014; Touboul et al., 2014). Optical coherence elastography (OCE) is another technique that can resolve depth dependent strain variation in corneal tissue subject to contact stress ex vivo (Ford et al., 2011, 2014a, 2014b). However, correlation between OCE mechanical strain and linear modulus of the cornea has not been reported yet (Ford et al., 2011, 2014a, 2014b). A modified non-contact OCE method can measure shear wave speed in the cornea but this method assumes the cornea as a linear, isotropic material (Ford et al., 2014b). Brillouin scattering uses an optical method to quantify the corneal biomechanical properties (Scarcelli et al., 2012a, 2012b, 2013). By measuring the spectral shift in the wavelength of the scattered light, Brillouin modulus can be measured in normal and crosslinked corneal tissue ex vivo (Scarcelli et al., 2012a, 2012b) and recently in vivo (Scarcelli et al., 2013). The Brillouin modulus may correlate with Young's modulus, though the correlation with in vivo properties requires further study. However, none of these investigational techniques described the stress stiffening effect or non-linear biomechanical properties of corneal stroma observed ex vivo (Ruberti et al., 2011). The US Food and Drug Administration (FDA) approved the ocular response analyzer (ORA) as the first tonometer to quantify corneal biomechanical properties in vivo (Hallahan et al., 2014a, 2014b; Piñero and Alcón, 2014; Sedaghat et al., 2010; Terai et al., 2012; Ventura et al., 2013). The device applies a high pressure air-puff over 20 ms and measures the intensity of infrared beam of light reflected off from the anterior corneal surface during deformation. The ORA measures changes in the biomechanical properties of the cornea through temporal and geometrical parameters that may be related to the biomechanical properties (Terai et al., 2012; Ventura et al., 2013; Hallahan et al., 2014a, 2014b). However after collagen crosslinking, these surrogate indices may not detect the change in corneal elastic and/or viscous moduli in vivo (Terai et al., 2012; Hallahan et al., 2014b). ORA indicates biomechanical weakening after refractive surgery (Terai et al., 2012) but no relationship has been established yet between the change in ORA variables and the corresponding change in modulus. Recently another air-puff applanation device, known as the Corvis-ST (CoST), was introduced (Bak-Nielsen et al., 2014; Faria-Correia et al., 2013; Shen et al., 2014; Tejwani et al.,

2014). Similar to the ORA, CoST applies an air-puff over 30 ms on the cornea. The air-puff causes deformation of the cornea and several images of a cross-section of the deforming cornea are acquired by a high speed Scheimpflug camera in CoST. Since CoST measures corneal deformation, mechanically relevant properties of the cornea may be derived through advanced mathematical algorithms and continuum biomechanical theories. In recent studies, such an algorithm estimated the change in the stiffness of the cornea after collagen crosslinking in keratoconus corneas in vivo (Sinha Roy et al., 2012, 2013) and also estimated linear corneal biomechanical properties in human eyes in vivo using a 2-D inverse finite element model (Kling et al., 2014). Therefore, this study aimed to develop an inverse algorithm and 3-D finite element (FE) representation of the whole eye for determination of nonlinear, fiber dependent biomechanical properties of the cornea. The biomechanical estimation included the contribution of the cellular matrix, in-plane collagen fibers and cross-links between lamellae that provide the depth-dependent shear resistance. In this study, known load (air-puff pressure) and deformation (corneal images acquired during deformation) determined the material properties. This is also known as an inverse method.

2.

Methods and materials

The magnitude of deformation of the cornea depends on the biomechanical properties of the cornea and globe (sclera, aqueous and vitreous humor), resistance offered by fatty tissues and muscles, magnitude of the transient airpressure incident on the cornea, cross-sectional area of the air-puff and air–cornea interaction (fluid–structure interaction) at the interface, thickness of tissue layers and intraocular pressure. Thus, estimation of the biomechanical properties of the cornea needs a multi-physics approach.

2.1.

Measurement of corneal deformation

The protocol was approved by the ethics committee of Narayana Nethralaya Multi-specialty hospital and written informed consent was obtained from the study participants. 10 eyes of 10 subjects with normal corneas were measured. The mean Goldmann applanation IOP was 15.270.7 mmHg. Mean flat and steep axis keratometry was 43.171.3 D and 43.771.4 D, respectively. Mean central corneal thickness was 527.2715.3 μm. CoST captured images of a cross-section of the deforming cornea during applanation. By simple edge detection routines, the anterior corneal edge was segmented and converted to Cartesian co-coordinates. A total of 140 images during a single air-puff test were captured. By subtracting the instantaneous coordinates of the edge with coordinates at the beginning of the air-puff applanation, the displacement of the anterior edge of the cornea was calculated. These calculations were performed real-time by the device and extracted for further analyses. For each eye, 3 replicate measurements were obtained and used for the inverse estimation of corneal biomechanical properties.

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

2.2.

CFD simulations

Normal pressure force and tangential shear stress caused deformation of the cornea during the applanation test. Both the deformation and interface forces can be modeled simultaneously by a two-way coupled fluid–structure interaction simulation. The relative magnitudes of the pressure force and shear stress were first estimated with a 2-D transient axisymmetric computational fluid dynamics (CFD) model (Fig. 1). The 2-D profile of the cornea was derived from the Gullstrand eye model (Pandolfi and Manganiello, 2006). The geometry was meshed with quadrilateral elements. Fig. 1 shows the dimen
 sions of the CFD model. No-slip boundary condition U ¼ 0 was applied at the cornea and sclera surface. Transient pressure [p(t)] boundary condition was applied at the nozzle exit, where t was the time (Fig. 1). A built-in pressure sensor measured the transient pressure profile, which was applied as a load using 2-D (function of space and time) cubic spline interpolation. Other boundary conditions applied to the model are shown in Fig. 1. The k–ω shear stress transport with low Reynolds number correction turbulence model (Fluent v14, Ansys Inc. USA) was used to solve the continuity and momentum equations. Air was simulated as an ideal gas due to expected low density variation in the flow domain. Fully coupled solver with second order upwinding was used to solve the continuity and momentum equations. Automatic time stepping was used with a maximum time step size of 0.1 ms. The CFD simulations were parameterized to vary the distance (L) between the exit nozzle and corneal apex: 10, 11 and 12 mm (Fig. 1). The results of the CFD simulation showed that the shear stresses were lower than normal pressure force by a factor of 100 at peak air-puff pressure at all magnitudes of L. Thus the shear stress was negligible compared to the normal pressure force and only the normal pressure force distribution obtained from CFD simulations was applied as a load on the anterior corneal surface in subsequent inverse FE simulations. Thus, inverse estimation

Fig. 1 – 2-D axisymmetric computational fluid dynamics model to compute pressure distribution on the anterior corneal surface. The air-puff originated from the nozzle. Uimplied that velocity vector along the boundary was zero. Un implied that the velocity normal to the edge was zero. p(t) was the temporal pressure boundary condition applied at the nozzle exit.

175

of corneal biomechanical properties did not require coupled fluid–structure interaction (cornea–air interface) simulations at the air–cornea interface

2.3.

3-D finite element model for inverse simulations

The method for 3-D corneal geometry generation was described previously in several studies (Sinha Roy and Dupps, 2011; Sinha Roy et al., 2012, 2013). The shape of the anterior and posterior corneal surfaces was derived from Scheimpflug imaging (Pentacam, Oculus Optikegrate Gmbh, Germany). The x, y and z Cartesian coordinates were used to generate 3-D anterior and posterior surfaces of the cornea in Solid Works 2010 (Dassault Systèmes SolidWorks Corporation, Waltham, USA). These 3-D surfaces were exported to TruGrid for 3-D volume meshing (XYZ Scientific Applications, Inc., Livermore, USA). Only the cornea was used for the simulations. An epithelium was modeled in each cornea with uniform thickness of 53 μm. During applanation of the cornea, the motion and deformation of the globe contributed to measured deformation (Kling and Marcos, 2013; Kling et al., 2014). However, the entire globe, fatty tissues and muscles cannot be imaged with current ocular imaging methods. Therefore, a simple spring–dashpot model coupled with a lumped mass was used to simulate the globe motion and deformation. The lumped mass accounted for the inertia of the globe since air-puff caused a dynamic (time dependent) deformation. The eye globe was approximated by a parallel network of linear spring (Kz) along the z axis and dashpot (μ) connected to the lumped mass (m) (Fig. 2A,B). The linear spring stiffness along x (Kx) and y (Ky) direction was assumed to be infinite (Fig. 2B) assuming that the primary deformation during applanation was along the z axis with minimal deformation along the x and y axis at the limbus (Kling et al., 2014). It was assumed that the optimized value of m would be similar to the mass of the eye globe, and that the resistance provided by the fatty tissue and muscles intrinsically contributed to the optimized magnitudes of the parallel network spring–dashpot model. The parallel network was connected to the limbal edge of the cornea (Fig. 2B). One parallel network was connected between each node along the limbus and ground (indicates zero displacement). Henceforth, all displacements mentioned in the text and in figures will be along the z-axis. The cornea was modeled as a hyperelastic, fiber dependent elastic material. Viscoelasticity of the cornea was not considered as the duration of the air-puff test was assumed to be too short to achieve thermodynamic equilibrium (Petsche and Pinsky, 2013). The distribution of air-pressure on the anterior corneal surface derived from the CFD simulations was applied as a boundary condition normal to the anterior epithelium surface (Fig. 2B). Since the distribution of pressure varied both spatially and temporally, cubic spline interpolation was used to interpolate between time and space, i.e., pressure as a function of time and radius measured from the geometric center of the cornea. The corn eal geometry was meshed with 8-noded hexahedral linear elements for an incompressible material. Four layers of elements were modeled through the corneal thickness. The finite element model was solved using Abaqus (Dassault Systèmes, Maryland, USA).

176

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

Fig. 2 – (A) Top-view of the cornea in the 3-D inverse finite element model of human eye. The shaded circular region shows a schematic of the area where the air-puff was incident on the cornea. Springs Kx, Ky and Kz were connected to the limbal surface; (B) cross-section view of the cornea with lumped mass (m) and parallel network spring–dashpot model to provide the resistance to globe deformation and motion caused by the air-puff pressure.

2.4.

Fiber-dependent material model of the cornea

Several fiber dependent hyperelastic material models have been proposed for the cornea (Pandolfi and Manganiello, 2006; Pinsky et al., 2005; Studer et al., 2010; Nguyen and Boyce, 2011). In this study, a combination of the free energy density was used to best replicate the in vivo fiber structure of the cornea. The hyperelastic material model is represented by free energy density (ψ) ψ ðCÞ ¼ ψ m ðI1 ; I3 Þ þ ψ f plane ðCÞ þ ψ f cross ðCÞ

ð1Þ

where I1 ¼ C:1 and I3 ¼ det[C] are the 1st and 3rd invariants of the deformation tensor. C was the right Cauchy-Green deformation tensor and det[C] represents determinant of the tensor. The isotropic energy density of the matrix (ψm) is described by
 ðI3  1Þ ψ m ðI1 ; I3 Þ ¼ c1 I1 3 þ c2 ðI1  3Þ2 þ D1

2

ð2Þ

 1=3 I1 ¼ I3 I1 is

where the distortional part of I1, I3 is the determinant of deformation gradient tensor and ci's are the material constants. D1 is the bulk modulus to enforce incompressibility and is assumed to be 10  6. The fiber energy density (ψf-plane) is described by (Petsche and Pinsky, 2013; Nguyen and Boyce, 2011): Z π Wf ðλf plane ÞDplane ðθÞdθ ð3Þ ψ f plane ðCÞ ¼ π

Wf plane ðλf plane Þ ¼

  2  k1 k1  exp k2 λf plane  1 2k2 2k2

ð4Þ

ψf-plane represented the energy density of in-plane lamellar pffiffiffiffiffiffiffiffiffiffiffiffiffiffi collagen fibers with stretch, λf-plane equal to A UCA, where T A¼ [cos θ, sin θ, 0] is the local direction vector of the fibers. k1 and k2 are the material constants. Dplane(θ) represented a weighted average of the energy density of the fiber families at each integration point of the element and also represented the change in the preferred direction of the fibers from orthogonal in the central cornea to circumferential near the

limbus (Nguyen and Boyce, 2011). Z π ψ f cross ðCÞ ¼ Wf ðλf cross ÞDcross ðθÞdθ π

Wf cross ðλf cross Þ ¼

  2  k1 k1  exp k2 λf cross 1 2k2 2k2

ð5Þ ð6Þ

ψf-cross represented the energy density of crosslink fibers pffiffiffiffiffiffiffiffiffiffiffiffiffi between the lamellae with stretch, λf-cross equal to BUCB, where B is 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 ξ (Studer et al., 2010; Nguyen and Boyce, 2011). Since in vivo imaging of this angle was not possible in human subjects, the angle ξ was assumed to be a function of depth and is modeled as follows (Petsche and Pinsky, 2013; Winkler et al., 2013): ξ ¼ 28:61

expð3:19½1 sÞ  1 expð3:19Þ 1

ð7Þ

where s is the non-dimensional local thickness (Petsche and Pinsky, 2013; Winkler et al., 2013). The angle ξ was evaluated at each element centroid. Dplane(θ) was kept equal to Dcross(θ) (Winkler et al., 2013). From the above equations, the Cauchy stress is determined by: 1 ∂ψ σ ¼ pffiffiffiffi F FT I3 ∂C

ð8Þ

where F is the deformation gradient tensor. The epithelium was modeled as an isotropic, hyperelastic, incompressible
 material ψ ðCÞ ¼ ψ m ðI1 ; I3 Þ only with c1 ¼ 5 kPa and c2 ¼ 0.0 kPa (Sinha Roy et al., 2014).

2.5.

Inverse algorithm for material property estimation

The biomechanical property variables that were optimized by the inverse algorithm were c1 (kPa), c2 (MPa), k1 (kPa), k2, Kz (N/ m), μ (Pa s) and m (g), where Pa is Pascal. c1 and c2 described the continuum behavior of the extracellular matrix. k1 and k2 described the continuum behavior of the collagen fibers. Kz, μ and m were a simplified aggregate representation of the

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

inertia of the globe due its mass, assumed linear viscoelastic resistance offered by the fatty tissues and orbital muscles. The inverse algorithm was based on the recent method developed for estimation of corneal stiffness change after collagen crosslinking in keratoconus corneas (Sinha Roy et al., 2012, 2013). The optimization function that is minimized by the inverse algorithm is as follows: Minimize χ 2 ¼

r X n X

ðdCoST  dFE Þ2

177

were then repeated by adding Gaussian noise (Mathworks, MA, USA) with a signal to noise ratio of 10.0 to the virtual corneal deformation data. Even with significant noise in the deformation data, estimated χ 2 was still of the order of 10  4 indicating once again that the inverse model was well-posed. In all the virtual tests of the inverse model, the difference between the optimized and assumed magnitudes of the material constants was less than 0.1%.

ð9Þ

2.6.

i¼1j¼1

where dCost and dFE are the displacement of the corneal edge measured by CoST and computed by the 3-D FE model, respectively. r and n are the total number of replicate measurements and the number of data points evaluated on the corneal edge, respectively. In this study, r was 3 for all the eyes. Function χ2 was minimized using the Levenberg–Marquardt algorithm. The gradient of the function χ2 with respect to the biomechanical property variables was computed numerically. Since minimization by Levenberg–Marquardt algorithm may yield a local minimum, the simulations were repeated with different initial estimates of the parameters to enhance the probability of finding the global minimum. To ensure uniqueness of the optimized solution, this method of optimization was recently adopted in another inverse modeling study (Pandolfi and Manganiello, 2006). When the relative change in the magnitude of all the variables was less than 1% between two successive iterations of the inverse method, the simulation was terminated. For each eye, the computations were repeated 3 times to evaluate the sensitivity of the estimated biomechanical property variables to measured IOP (IOP¼13, 15 and 17 mmHg). This would also account for variability in IOP measurements which are effected by corneal biomechanical properties (Ruberti et al., 2011). All biomechanical variables were reported as mean7standard deviation for each IOP group. An ill-posed inverse model can result in multiple optimized biomechanical properties for the same measured corneal deformation (Kok et al., 2014). To ensure that the inverse model was well defined, “virtual” corneal deformation data was computed by assuming arbitrarily chosen magnitudes of the material constants. Then, the virtual deformation data was used an input in the inverse model with different initial estimates of material constants to determine whether the inverse model was able to converge to the assumed magnitudes of the material constants. Table 1 summarizes the results of the virtual tests of the inverse model. Since χ 2 was of the order of 10  4 (rows 1–3 in Table 1), it indicated that the inverse model was well posed to converge to the assumed magnitudes of material constants. The virtual tests described in rows 1–3 of Table 1

Comparison of corneal biomechanical properties

The air-puff model computed corneal deformation combined with the globe deformation. Therefore, it was necessary to assess the biomechanics of each cornea independent of globe mechanics. Further, analyzing all the property variables among eyes may get complex in a clinical setting. Therefore, direct comparison of magnitude of deformation of the corneas based on the estimated property variables independent of the globe property variables may be a better method. Thus, virtual corneal inflation simulations were performed similar to ex-vivo corneal button inflations tests (Ruberti et al., 2011). The corneal apex rise was computed as a function of uniform pressure increase applied to the posterior surface of the cornea in a forward FE model by restricting the displacement of the limbus completely. Thus for each eye, three different corneal apex rise vs. increase in uniform pressure applied to the posterior surface were obtained. This also simplified the comparison of the stiffness of the different corneas to a single curve as a function of just IOP and 3-D geometry for each cornea instead of comparing the individual property variables. Thus if two corneas had the same thickness but different estimated biomechanical property variables, the more compliant the cornea, the more the corneal apical rise in the virtual inflation tests.

3.

Results

Fig. 3 shows the axial velocity contour plot at time t¼ 0.015 s. The peak magnitude of the axial velocity was 195 m/s. Fig. 4A shows the temporal distribution of air pressure at the nozzle exit and deformation amplitude (displacement of the corneal apex) measured by CoST in one eye. There was negligible phase difference between the peak air pressure and peak apical displacement. This indicated that corneal viscous properties did not contribute significantly to the measured displacement. Fig. 4B shows the spatial distribution of pressure at different times during the applanation. For comparison, the pressure values at each time were normalized with respect to maximum pressure (at distance ¼0.0 m) at the

Table 1 – Results of virtual testing of the inverse model. The magnitude of the optimization function, χ 2 , with and without added Gaussian noise is also shown below. Assumed material constants c1 (kPa), c2 (MPa), k1 (kPa), k2, Kz (N/m), μ (Pa s) and m (g)

Initial value c1 (kPa), c2 (MPa), k1 (kPa), k2, Kz (N/m), μ (Pa s) and m (g)

χ 2 (mm2)

χ 2 with Gaussian noise (mm2)

60, 3, 20, 300, 1.4, 4.0, 0.04 70, 5, 35, 200, 1.5, 5.5, 0.06 80, 6, 45, 350, 1.7, 4.7, 0.05

120, 6, 40, 100, 2.0, 2.0, 0.01 140, 10, 25, 150, 2.0, 2.0, 0.01 160, 12, 70, 100, 2.0, 2.0, 0.01

1  10  4 1.25  10  4 1.1  10  4

1.9  10  4 2.25  10  4 1.67  10  4

178

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

Fig. 3 – Contour plot of axial velocity distribution at peak airpuff pressure.

Fig. 4 – (A) Measured deformation amplitude (displacement of corneal apex) from 1 eye. The air-puff pressure at the nozzle exit was superimposed; (B) normalized air pressure incident on the anterior corneal surface vs. distance from geometric center of the cornea. The actual pressure at a radius was divided by the pressure magnitude at corneal geometric center (0.0 on the x-axis) at the same time to calculate the normalized pressure.

same time in Fig. 4B. There was no perceptible difference between the pressure distributions at different time points. Fig. 5A shows the measured (CoST) vs. estimated deformation amplitude for the same eye shown in Fig. 4A. The solid lines were the estimated (from inverse simulations) deformation amplitudes at different IOP's¼ 13, 15 and 17 mmHg. The estimated deformation amplitudes for each IOP was the same as that for the other two IOP's as the inverse model simply perturbed the material constants to optimize the estimated deformation amplitude to the measured deformation amplitude. Fig. 5B shows the measured and estimated displacement along the anterior edge of the cornea at different time points. In both figures, the estimated displacements confirmred well with the measured displacements by CoST. Fig. 6A–C shows the estimated corneal apical displacement as a function of uniform pressure increase from 0 to 60 mmHg for all the eyes plotted together. In these figures, the plots were obtained by restraining the limbus as explained in the methods section under “comparison of

Fig. 5 – (A) Estimated deformation amplitude from 3-D inverse finite element model of the eye shown in Fig. 4A. The red, green and blue lines show the estimated deformation amplitude for different intraocular pressures (13, 15 and 17 mmHg); (B) comparison of actual and estimated (model) displacement of the anterior corneal edge at different times for the same eye reported in Fig. 4A and in inset A. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

179

corneal biomechanical properties” and by using each of the three estimated sets of biomechanical property variables one set at a time. The non-linear response of the cornea was evident from these results. Non-linear regression was performed and plotted in Fig. 7 shows inverse simulations estimated a stiffer cornea, when the IOP was lower and vice versa for the same deformation amplitude curve. For

Fig. 7 – Comparison of simulated intraocular pressure on the posterior corneal surface vs. displacement of the corneal apex using the non-linear regression computed in Fig. 6. For the same deformation amplitude curve, higher the intraocular pressure, the more compliant the cornea. example, the displacement of the corneal apex was 0.1 mm when a pressure of 15 mmHg was applied on the posterior corneal surface, if the estimated property variables were derived with a measured in vivo IOP of 15 mmHg (Fig. 7). However, the same was  0.09 mm when a pressure of 15 mmHg was applied on the posterior corneal surface, if the estimated property variables were derived with a measured in vivo IOP of 13 mmHg (Fig. 7). The average values of the biomechanical variables are reported in Table 2.

4.

Fig. 6 – (A), (B) and (C) show the simulated intraocular pressure on the posterior corneal surface vs. displacement of the corneal apex similar to corneal button inflations test. (A), (B) and (C) were plotted using estimated biomechanical property variables at different intraocular pressures of 13, 15 and 17 mmHg, respectively. Non-linear regression shown by solid lines was performed.

Discussion

These are exciting times as there are several techniques in development for quantification of corneal biomechanical properties (Torricelli et al., 2014; Touboul et al., 2014; Ford et al., 2011, 2014a, 2014b; Scarcelli et al., 2012a, 2012b, 2013). However, in vivo quantification of non-linear, fiber dependent properties of the cornea is needed. In this study, a method was established and validated for quantification of the corneal biomechanical properties. A recent study used a 2-D inverse finite element model to estimate corneal linear elastic and viscous properties from CoST air-puff deformation (Kling et al., 2014). Significant improvements that were achieved in this study over the 2-D method (Kling et al., 2014) are as follows: (a) 3-D patient specific tomography was implemented; (b) fiber (in-plane and between the lamellae) distribution and non-linear elastic properties were implemented; (c) integration of the globe, fatty tissues and orbital muscles into a lumped parameter, parallel network, viscoelastic material formulation was achieved; (d) simultaneous optimization of the biomechanical parameters instead of sequential approach (Kling et al., 2014) was achieved. Unlike the 2-D study where corneal viscous properties were considered, the results from this study (Fig. 4A) support the hypothesis that duration of the air-puff deformation may not be long enough to sense phase lag between pressure and deformation (Petsche and Pinsky, 2013). At peak deformation (Fig. 4A at time¼ 0.015 s), the globe deformation had little contribution to the overall measured deformation (around 1/10th) (Kling et al., 2014). However, the globe deformation continued even after the airpuff pressure reduced to zero indicating some viscous effects

180

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

Table 2 – Summary of estimated biomechanical variables for all the eyes. Mean7standard deviation of the variables are reported. IOP (mmHg)

c1 (kPa)

c2 (Mpa)

k1 (kPa)

k2 (kPa)

Kz (N/m)

μ (Pa s)

m (g)

χ2 (mm2)

13 15 17

76723.8 73722.8 69722

5.3371.18 4.8179.97 4.8271.14

0.4170.04 0.3970.05 0.3970.04

276744 308776 293748

1.6070.29 1.6170.33 1.6370.36

6.3673.72 5.9973.72 5.9773.59

0.05770.009 0.05870.011 0.05970.014

0.05670.007 0.05670.005 0.06170.009

after corneal deformation ceased, e.g., even at time t¼ 0.03 s, the deformation amplitude was greater than 0.0 mm (Fig. 4A). Future applications in degenerate and surgically treated corneas may be possible with this 3-D inverse finite element model because such corneas have more asymmetry shape, which may impact the uniqueness of the optimized material parameters. The mean lumped mass (m) was estimated to be 0.058 g (Table 1). Since the limbal surface had 264 nodes, the total estimated lumped mass was 15.31 g ( ¼264  0.058). This was remarkably similar to globe volume (mean of 14.83 g in females and 15.99 in males) measured in adult humans (Forbes et al., 1985). This study also evaluated the impact of measured intraocular pressure on the estimated properties. For the same deformation amplitude curve, higher the intraocular pressure, the more compliant the cornea (Table 1). Thus, accurate knowledge of intraocular pressure may be a confounding factor in assessment of post-refractive surgery and crosslinking induced biomechanical changes. Measured IOP may vary depending on the device used as well. To understand the confounding effect of IOP measurement, all the eyes were simulated at different IOP's to define an upper and lower bound of the biomechanical parameters within a 72 mmHg variation in the measured IOP (Table 2). With continuous advancements in IOP measurement, these bounds will require re-estimation. Further layer specific biomechanical information, e.g., Bowman's layer, may be included in the 3-D geometry of the cornea as data becomes available in future studies on the relative stiffness of individual layers of the cornea. A limitation of this study was that the estimated biomechanical property variables could not be validated for the same corneas using conventional ex vivo methods due to the in vivo setting. Hence, optimized values of the property variables were only referred to as “estimates”. Another limitation of the study was that CoST reported deformation only along one crosssection, which was along horizontal meridian (or 01 meridian) of the cornea at the start of the corneal deformation measurement. This imaging plane remained fixed in the device during the measurement. However, the deformation due to the air-puff may have occured in 3 dimensions depending on the level of asymmetry of the 3-D cornea and this necessitated the use of 3D finite element model. In the inverse model, the axial deformation (along z axis in Fig. 2B) measured by CoST in the fixed imaging plane was compared with the corresponding deformation of the edge on the anterior surface of the cornea along the horizontal plane. Only the anterior surface was used for the optimization as the measurement of posterior surface may be suspect in surgically treated corneas, which would limit the applicability of this inverse model post-treatment (Ruberti et al., 2011). With further advancement of the imaging method, it may

be possible to include more imaging planes in the inverse simulations and more realistic estimates of the biomechanical property variables may be obtained. Future advancements in air-puff deformation measurements may also allow adjustment to the duration of the air-puff and peak pressure applied to anterior surface of the cornea in order to quantify viscous properties. Future advancements in in vivo imaging of fiber distribution (Fukuda et al., 2013) may also result in improved estimates of biomechanical property variables. Another limitation of the present study was that the spatial distribution of pressure profile may vary during deformation of the cornea and cannot be accounted by only CFD models. However, the present study provided evidence to support that spatial variation of the pressure profile may have minimal impact on the results from the inverse model. Firstly, most of the kinetic energy of the airpuff was concentrated in the central 2.5 mm zone even at the corneal plane. This is shown in Fig. 3 and also in another recent study (Ariza-Gracia et al., 2015). Secondly since significant shear stresses will be concentrated in regions of steep velocity gradients, the maximal shear stress was localized at the central 2.5 mm zone at the cornea. Thirdly, even at peak applanation pressure and corresponding maximal air-puff velocity, the shear stress was lower by an order of magnitude of 100 indicating that shear stress may have little role in influencing the temporal and spatial change in corneal deformation. To summarize, this study established an inverse modeling approach to estimate the non-linear, fiber dependent, biomechanical properties of the cornea in an in vivo setting. The application and utility of the estimated property variables needs to be tested in larger cohorts of normal and disease subjects similar to studies conducted with the ORA (Piñero and Alcón, 2014).

Financial disclosures ASR has intellectual property on biomechanical modeling of the human eye through Cleveland Clinic Innovations, USA. ASR has received research funding from Carl Zeiss Inc., USA, Avedro Inc., USA and Topcon Inc., USA in past three years on biomechanical modeling of cornea. MK, HM and RS have no financial interests to declare.

Acknowledgments The authors wish to thank Oculus Optikegerate Gmbh, Germany for providing access to the research software that performed the corneal edge detection in CoST and reported the edge shape in Cartesian co-ordinates.

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

r e f e r e n c e s

´ ., Zurita, J.F., Pin˜ero, D.P., Rodriguez-Matas, J.F., Ariza-Gracia, M.A Calvo, B., 2015. Coupled biomechanical response of the cornea assessed by non-contact tonometry. A simulation study. PLoS One 10, e0121486. Bak-Nielsen, S., Pedersen, I.B., Ivarsen, A., Hjortdal, J., 2014. Dynamic scheimpflug-based assessment of keratoconus and the effects of corneal cross-linking. J. Refract. Surg. 30 (6), 408–414. Faria-Correia, F., Ramos, I., Valbon, B., Luz, A., Roberts, C.J., Ambro´sio Jr., R., 2013. Scheimpflug-based tomography and biomechanical assessment in pressure-induced stromal keratopathy. J. Refract. Surg. 29 (5), 356–358. Ford, M.R., Dupps Jr, W.J., Rollins, A.M., Sinha Roy, A., Hu, Z., 2011. Method for optical coherence elastography of the cornea. J. Biomed. Opt. 16 (1), 016005. Ford, M.R., Rollins, A.M., Dupps, W.J. Jr., 2014a. Quantitative in vivo corneal elastography by Doppler shear wave imaging. ARVO Meeting Abstracts, vol. 55, pp. 3724. Ford, M.R., Sinha Roy, A., Rollins, A.M., Dupps Jr., W.J., 2014b. Serial biomechanical comparison of edematous, normal, and collagen crosslinked human donor corneas using optical coherence elastography. J. Cataract Refract. Surg. 40 (6), 1041–1047. Forbes, G., Gehring, D.G., Gorman, C.A., Brennan, M.D., Jackson, I.T., 1985. Volume measurements of normal orbital structures by computed tomographic analysis. Am. J. Roentgenol. 145 (1), 149–154. Fukuda, S., Yamanari, M., Lim, Y., Hoshi, S., Beheregaray, S., Oshika, T., Yasuno, Y., 2013. Keratoconus diagnosis using anterior segment polarization-sensitive optical coherence tomography. Invest. Ophthalmol. Vis. Sci. 54 (2), 1384–1391. Hallahan, K.M., Sinha Roy, A., Ambrosio Jr., R., Salomao, M., Dupps Jr., W.J., 2014a. Discriminant value of custom ocular response analyzer waveform derivatives in keratoconus. Ophthalmology 121 (2), 459–468. Hallahan, K.M., Rocha, K., Sinha Roy, A., Randleman, J.B., Stulting, R.D., Dupps, W.J. Jr., 2014b. Effects of corneal collagen crosslinking on ocular response analyzer waveform-derived variables in keratoconus and post-LASIK ectasia. ARVO Meeting Abstracts, vol. 55, pp. 3722. Kling, S., Bekesi, N., Dorronsoro, C., Pascual, D., Marcos, S., 2014. Corneal viscoelastic properties from finite-element analysis of in vivo air-puff deformation. PLoS One 9 (8), e104904. Kling, S., Marcos, S., 2013. Contributing factors to corneal deformation in air puff measurements. Invest. Ophthalmol. Vis. Sci. 54 (7), 5078–5085. Kok, S., Botha, N., Inglis, H.M., 2014. Calibrating corneal material model parameters using only inflation data: an ill-posed problem. Int. J. Numer. Method Biomed. Eng. 30, 1460–1475. Nguyen, T.M., Aubry, J.F., Touboul, D., Fink, M., Gennisson, J.L., Bercoff, J., Tanter, M., 2012. Monitoring of cornea elastic properties changes during UV-A/riboflavin-induced corneal collagen cross-linking using supersonic shear wave imaging: a pilot study. Invest. Ophthalmol. Vis. Sci. 53 (9), 5948–5954. Nguyen, T.D., Boyce, B.L., 2011. An inverse finite element method for determining the anisotropic properties of the cornea. Biomech. Model. Mechanobiol. 10 (3), 323–337. Pandolfi, A., Manganiello, F., 2006. A model for the human cornea: constitutive formulation and numerical analysis. Biomech. Model. Mechanobiol. 5 (4), 237–246. Pedersen, I.B., Bak-Nielsen, S., Vestergaard, A.H., Ivarsen, A., Hjortdal, J., 2014. Corneal biomechanical properties after LASIK, ReLEx flex, and ReLEx smile by Scheimpflug-based dynamic tonometry. Graefes Arch. Clin. Exp. Ophthalmol. 252 (8), 1329–1335.

181

Petsche, S.J., 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), 1101–1113. Pin˜ero, D.P., Alco´n, N., 2014. In vivo characterization of corneal biomechanics. J. Cataract Refract. Surg. 40 (6), 870–887. Pinsky, P.M., van der Heide, D., Chernyak, D., 2005. Computational modeling of mechanical anisotropy in the cornea and sclera. J. Cataract Refract. Surg. 31 (1), 136–145. Roberts, C., 2002. Biomechanics of the cornea and wavefrontguided laser refractive surgery. J. Refract. Surg. 18 (S5), S589–S592. Ruberti, J.W., Sinha Roy, A., Roberts, C.J., 2011. Corneal biomechanics and biomaterials. Annu. Rev. Biomed. Eng. 13, 269–295. Scarcelli, G., Pineda, R., Yun, S.H., 2012a. Brillouin optical microscopy for corneal biomechanics. Invest. Ophthalmol. Vis. Sci. 53 (1), 185–190. Scarcelli, G., Yun, S.H., 2012b. In vivo Brillouin optical microscopy of the human eye.Opt. Express. 20 (8), 9197-9202. Scarcelli, G., Kling, S., Quijano, E., Pineda, R., Marcos, S., Yun, S.H., 2013. Brillouin microscopy of collagen crosslinking: noncontact depth-dependent analysis of corneal elastic modulus. Invest. Ophthalmol. Vis. Sci. 54 (2), 1418–1425. Sedaghat, M., Naderi, M., Zarei-Ghanavati, M., 2010. Biomechanical parameters of the cornea after collagen crosslinking measured by waveform analysis. J. Cataract Refract. Surg. 36 (10), 1728–1731. Shen, Y., Zhao, J., Yao, P., Miao, H., Niu, L., Wang, X., Zhou, X., 2014. Changes in corneal deformation parameters after lenticule creation and extraction during small incision lenticule extraction (SMILE) procedure. PLoS One 9 (8), e103893. Sinha Roy, A., Dupps Jr., W.J., 2011. Patient-specific modeling of corneal refractive surgery outcomes and inverse estimation of elastic property changes. J. Biomech. Eng. 133 (1), 011002. Sinha Roy, A., Fant, B., Rocha, K., Dupps, W.J., Jr., 2012. Estimation of modulus change after corneal crosslinking using multiple post-cxl topographies and inverse finite element. ARVO Meeting Abstracts, vol. 53, pp. 6896. Sinha Roy, A., Rocha, K.M., Randleman, J.B., Stulting, R.D., Dupps Jr., W.J., 2013. Inverse computational analysis of in vivo corneal elastic modulus change after collagen crosslinking for keratoconus. Exp. Eye Res. 113, 92–104. Sinha Roy, A., Dupps Jr, W.J., Roberts, C.J., 2014. Comparison of biomechanical effects of small-incision lenticule extraction and laser in situ keratomileusis: finite-element analysis. J. Cataract Refract. Surg. 40, 971–980. Studer, H., Larrea, X., Riedwyl, H., Bu¨chler, P., 2010. Biomechanical model of human cornea based on stromal microstructure. J. Biomech. 43 (5), 836–842. Tejwani, S., Shetty, R., Kurien, M., Dinakaran, S., Ghosh, A., Sinha Roy, A., 2014. Biomechanics of the cornea evaluated by spectral analysis of waveforms from ocular response analyzer and Corvis-ST. PLoS One 9 (8), e97591. Terai, N., Raiskup, F., Haustein, M., Pillunat, L.E., Spoerl, E., 2012. Identification of biomechanical properties of the cornea: the ocular response analyzer. Curr. Eye Res. 37 (7), 553–562. Torricelli, A.A., Ford, M.R., Singh, V., Santhiago, M.R., Dupps Jr., W.J., Wilson, S.E., 2014. BAC-EDTA transepithelial riboflavin-UVA crosslinking has greater biomechanical stiffening effect than standard epithelium-off in rabbit corneas. Exp. Eye Res. 125, 114–117. Touboul, D., Gennisson, J.L., Nguyen, T.M., Robinet, A., Roberts, C.J., Tanter, M., Grenier, N., 2014. Supersonic shear wave elastography for the in vivo evaluation of transepithelial corneal collagen cross-linking. Invest. Ophthalmol. Vis. Sci. 55 (3), 1976–1984.

182

journal of the mechanical behavior of biomedical materials 48 (2015) 173 –182

Ventura, B.V., Machado, A.P., Ambro´sio Jr., R., Ribeiro, G., Arau´jo, L.N., Luz, A., Lyra, J.M., 2013. Analysis of waveform-derived ORA parameters in early forms of keratoconus and normal corneas. J. Refract. Surg. 29 (9), 637–643. Winkler, M., Shoa, G., Xie, Y., Petsche, S.J., Pinsky, P.M., Juhasz, T., Brown, D.J., Jester, J.V., 2013. Three-dimensional distribution of

transverse collagen fibers in the anterior human corneal stroma. Invest. Ophthalmol. Vis. Sci. 54 (12), 7293–7301. Wollensak, G., Spoerl, E., Seiler, T., 2003. Stress–strain measurements of human and porcine corneas after riboflavin-ultraviolet-A-induced cross-linking. J. Cataract Refract. Surg. 29 (9), 1780–1785.