Patient-Specific Computational Modeling of Keratoconus Progression and Differential Responses to Collagen Cross-linking

Abstract

Purpose.: To model keratoconus (KC) progression and investigate the differential responses of central and eccentric cones to standard and alternative collagen cross-linking (CXL) patterns.

Methods.: Three-dimensional finite element models (FEMs) were generated with clinical tomography and IOP measurements. Graded reductions in regional corneal hyperelastic properties and thickness were imposed separately in the less affected eye of a KC patient. Topographic results, including maximum curvature and first-surface, higher-order aberrations (HOAs), were compared to those of the more affected contralateral eye. In two eyes with central and eccentric cones, a standard broad-beam CXL protocol was simulated with 200- and 300-μm treatment depths and compared to spatially graded broad-beam and cone-centered CXL simulations.

Results.: In a model of KC progression, maximum curvature and HOA increased as regional corneal hyperelastic properties were decreased. A topographic cone could be generated without a reduction in corneal thickness. Simulation of standard 9-mm-diameter CXL produced decreases in corneal curvature comparable to clinical reports and affected cone location. A 100-μm increase in CXL depth enhanced flattening by 24% to 34% and decreased HOA by 22% to 31%. Topographic effects were greatest with cone-centered CXL simulations.

Conclusions.: Progressive hyperelastic weakening of a cornea with subclinical KC produced topographic features of manifest KC. The clinical phenomenon of topographic flattening after CXL was replicated. The magnitude and higher-order optics of this response depended on IOP and the spatial distribution of stiffening relative to the cone location. Smaller diameter simulated treatments centered on the cone provided greater reductions in curvature and HOA than a standard broad-beam CXL pattern.

Keratoconus (KC) is a disease that is characterized clinically by a region of abnormally high curvature, reduced corneal thickness, and progressive corneal topographic irregularity. ^1,2 Tensile testing of corneal strips cut from KC keratoplasty buttons has also demonstrated 50% to 60% reductions in the strain-dependent corneal elastic modulus relative to normal cadaveric tissue. ^3,4 These findings, taken with some of the classic clinical signs of KC, such as corneal stress lines (Vogt striae) and rupture of Descemet's membrane (hydrops), support the concept of biomechanical failure as a central pathologic feature of the disease. Although the interactions of genetic, ^5,6 biochemical, ^7,8 and ultrastructural ^9,10 abnormalities that ultimately contribute to material instability and corneal deformation under the ocular loading forces are complex, we propose that differential biomechanical weakening in the area of the cone is a common pathway through which these biological factors influence corneal shape. By corollary, we propose that topographic progression of KC can be simulated as a progressive focal reduction in material elastic modulus without an absolute requirement for specific knowledge of the genetic, biochemical, and ultrastructural processes that contribute to ongoing material failure. 

Corneal collagen cross-linking (CXL) is a relatively new tissue-sparing technique for the treatment of KC. ^{11 –16} In this technique, a photosensitizing agent containing riboflavin is applied to the corneal surface at regular time intervals for a total time of 30 minutes followed by continuous application of a broad ultraviolet (UV-A) beam for another 30 minutes. ^11,13,15 Before riboflavin is applied, the epithelium is typically removed from the central 5- to 9-mm-diameter zone of the cornea to facilitate penetration of riboflavin into the stroma. ^11,13,15 Wollensak et al. ¹⁶ reported a nearly 3.5 fold increase in Young's modulus of human corneas after CXL ex vivo. Histology and optical coherence tomography (OCT) suggest that the morphologic and optical changes associated with CXL are typically restricted to a depth of 200 to 300 μm from the anterior surface with an irradiance of 3 mW/cm². ^14,16,17 In addition to the intended curvature-stabilizing effects of CXL, unexpected reductions in steep corneal curvature and improvements in uncorrected and corrected distance visual acuity have been reported. ^11,13 A recent clinical study has also suggested that optical outcomes in the central cornea can differ significantly after CXL, depending on the location of the cone (defined as the locus of maximum tangential curvature) relative to the optic axis. ¹⁵  

The finite element model (FEM) can be a useful tool for providing insight into the effects of such material property alterations on the geometry and optical performance of the cornea. A small number of studies have explored the effect of corneal thinning and elastic modulus reduction in KC. ^{18 –20} However, these studies have not accounted for major clinical features of KC such as astigmatism and higher order aberrations, nor have they incorporated patient-specific clinical tomography. No previous study has attempted to model the effects of CXL on KC. We recently analyzed the effect of variations in overall corneal elastic modulus in a two-dimensional whole-eye model of myopic laser in situ keratomileusis (LASIK) and demonstrated a greater tendency toward biomechanical steepening of the central cornea (and residual myopia) in a model with weaker hyperelastic corneal properties. ²¹ Conversely, an otherwise identical model with a stiffer cornea demonstrated a greater tendency toward central corneal flattening due to a redistribution of strain from the central to the peripheral cornea. Although the model assumed a uniform stiffening of the entire cornea, the results suggested a potential biomechanical mechanism for the reduction of corneal steepening observed in many patients treated with CXL. ²¹  

The purposes of this study were (1) to develop a method for three-dimensional modeling of KC progression and for modeling the effects of corneal collagen stiffening, (2) to evaluate changes in curvature and optical aberrations of the anterior surface of the KC cornea as a function of magnitude and depth of corneal stiffening, and (3) to compare model results to clinical topographic observations of curvature reduction. Finally, the study also evaluated novel CXL approaches employing reduced treatment diameters centered on the steepest topographic point. During the course of these studies, we hypothesized that a smaller, focal treatment centered on the topographic cone might provide enhanced flattening effects and more favorable changes in higher order aberrations. A smaller treatment zone may expose the patient to less treatment risk by virtue of a smaller epithelial defect with less chance of persistence, associated haze, ^22,23 and UV exposure of the corneal endothelium, limbal stem cells, and intraocular tissues. We also modeled the effect of introducing a spatially graded treatment (which is feasible clinically with a modified ultraviolet intensity profile) incorporating a less abrupt radial transition of elastic properties than might be expected from a more uniform ultraviolet exposure. As a substrate for the model, this study uses data from four clinical corneal tomography scans: (1) two eyes from an asymmetric KC patient whose less-affected eye is progressed to a state resembling the more severely affected eye via simulated changes in elastic properties and corneal thickness and (2) two KC eyes with different cone loci to assess differences in the response to CXL as a function of cone location and depth of CXL. 

Methods

Model Geometry and Applied Pressure

The anterior and posterior corneal surface topography maps for the asymmetric KC patient and the two eyes with central and eccentric cone were obtained from a commercial anterior segment imaging system (Pentacam, ver. 1.61; Oculus Optikgeräte GmbH, Wetzlar, Germany). The x, y, and z coordinates from the elevation maps of anterior and posterior surface were interpolated using linearized orthogonal Zernike polynomials up to 10th order having 66 terms and with a 5.75 mm normalization radius. The mean error, defined as square root of the sum of difference between in vivo and Zernike fit elevation co-ordinate (z) divided by the number of data points, for all eyes and surfaces was 0.29 + 1.29 × 10⁻⁵ μm. The approximated elevation data (z) from Zernike fit were used to obtain a three-dimensional surface with a commercial CAD package (ProEngineer Wildfire; PTC, Needham, MA). Each point along a radius was connected by cubic spline to form a three-dimensional curvilinear edge. Multiple edges that form the shape of a corneal surface were then blended to obtain a three-dimensional surface for the anterior and posterior cornea. These surfaces (anterior and posterior) were then joined at the limbus to form a three-dimensional solid, representing the in vivo cornea. To simulate an unconstrained limbus, the sclera (Figs. 1A, 1B) was modeled as a layer of tissue with an external diameter of 23 mm, an axial eye length of 22 mm and a uniform thickness of 850 μm connected to the corneal three-dimensional structure at the limbus. ²¹  

Figure 1.

 

(A, B) Whole-eye FEM model showing the cornea and sclera. (C) A schematic representation of the model showing the weaker KC zone and the CXL zone in a standard 9-mm-diameter simulated treatment.

The internal surface of the whole-eye model was loaded with a uniform pressure representing the intraocular pressure (IOP). Ex vivo measurements have shown an apparent increase in IOP after CXL related to an increase in corneal elastic modulus. ^{24 –26} Computational simulation also suggests that IOP measurement by applanation is sensitive to variations in corneal elastic modulus. ²⁷ Clinical studies have also suggested a small increase in IOP. ²⁸ Because the slight changes in IOP measured clinically may be an artifact attributable to the stiffening effects of CXL and would likely have little impact on the simulations, even if they represent true changes in IOP, it was kept unchanged at a nominal value of 15 mm Hg from preoperative to postoperative states for most CXL simulations. To address the possibility of a nominal increase in true loading pressure after CXL and any impact on the comparability of these simulations to clinical reports, we investigated the pressure sensitivity of the post-CXL model, as described later in this section. Meshing of the corneoscleral whole-eye model was performed with a mesh generator (TrueGrid; XYZ Scientific Applications, Inc., Livermore, CA), using 8-noded linear, hexahedral elements for an incompressible material (Fig. 1A). The number of elements in the FEM mesh was varied to evaluate the dependency of computed stress and strain on it. Mesh independence was confirmed when a less than 0.5% change in strain and displacement values was obtained. In the clinical setting, corneal topography is measured at IOP (> 0) and is distinct from the unloaded shape that would be obtained at IOP of 0 mm Hg. To solve for the undeformed (0 load) state, a custom inverse model was developed using an FEM analysis package (ABAQUS; Simulia, Inc., Dassault Systèmes, Providence, RI) and Python scripting language. ²⁹ The full Newton's method was used to solve the FEM equations. 

Material Properties

The normal cornea is an anisotropic tissue that, as IOP increases, shows stress stiffening due to a complex collagen fibrillar distribution. ^9,30 Normal corneas have been shown to exhibit a diamond-shaped arrangement of collagen fibers, with orthogonal alignment in the central cornea transitioning to circumferential alignment near the limbus. ³⁰ The fiber packing density follows a definitive pattern in accordance with the fiber orientation. However, in KC corneas, the preferred collagen fibril orientation seen in normal corneas is disrupted, especially in the region of the cone, which could favor a reduced local elastic modulus. ^9,10,30 Since the fiber packing density in KC corneas does not follow the stereotypical spatial distribution of normal corneas and is more random, we adopted a material model for the KC cornea that is fiber independent and thus spatially isotropic. Stress versus strain data derived from ex vivo testing of normal corneas ¹⁶ was fitted to a reduced polynomial material model, W = C ₁₀(I ₁ − 3) = C ₂₀(I ₁ − 3)², where W is the strain energy potential, C ₁₀ and C ₂₀ are constants to be specified later, and I ₁ is the strain invariant. From W, stress tensor (σ) is given by δW/δε, where ε is strain. Elastic modulus can be calculated as δσ/δε. The sclera elastic property was set at four times the elastic modulus of a normal cornea. ³¹  

To simulate the pathomechanics of KC, we implemented a circular zone of reduced elastic modulus in the model (as illustrated in the example in Fig. 1C). ¹⁸ Since the cornea is weaker in this circular zone, the constants C ₁₀ and C ₂₀ for normal cornea were reduced, as a function of radius and independent of corneal thickness, by α equal to  

 

where m ₁ and m ₂ are regressed constants, r is the local radius, and 0 < α ≤ 1. Thus in pre-CXL models, the decrease in elastic modulus in the weak zone was modeled as an exponential decay from the edge to the center of the circular weak zone, such that in the region outside the weak zone, α = 1. The determination of m ₁ and m ₂ is described below for each simulation. 

Simulation of KC Progression

Corneal tomography and tangential curvature of the right and left eyes of an asymmetric KC patient were used to generate patient eye-specific models as described above (Figs. 2A, 2B, respectively). The leftmost columns of Table 1 show the simulated keratometry (SimK), maximum axial curvature (K _max), central corneal thickness (CCT), and coordinates (x, y) of the location of minimum corneal thickness of the two eyes. The left eye (K _max = 47.10 D) showed an early stage of KC with a distinct inferior cone. In the right eye, the tangential curvature map showed irregular astigmatism without a discrete locus of steepening (Fig. 2A), whereas pachymetry demonstrated an inferiorly decentered thickness minimum. Since a distinct zone of steepening was not present in the right eye, we centered the weak zone at the locus of minimum corneal thickness in the right eye. To account for the difference in thickness between the two eyes and to evaluate the effect of a generalized reduction in thickness on KC progression, we shifted the posterior surface of the right eye axially, to reduce the CCT from 563 μm to the 536 μm measured in the left eye. A circular weak zone with center at (x, y) of the right eye was implemented. At r = r ₀, α = 0.5 and dr = 0.075 in equation 1. The diameter of the weak zone was estimated to be 7 mm by extrapolating from ex vivo evidence of similar diameters of collagen fibril disorganization (i.e., for r ≥ 3.5 mm, α = 1). ¹⁰ With these conditions, three pairs of values for m ₁ and m ₂ were calculated, such that the volume-averaged reduction in elastic modulus (1 − α̃) in the weak zone was 0.1, 0.3, and 0.45. Then, an IOP of 15 mm Hg was applied as a uniform pressure to the 0-load model of the right eye with the altered corneal elastic modulus  

 

to investigate the impact of these material property changes on possible progression of KC in the right eye. The value of C ₁₀ and C ₂₀ was 140,938 and 80,134 Pa, respectively. The effect of focal corneal weakening was also investigated without any change in corneal thickness to assess the role of elastic modulus abnormalities independently. 

Figure 2.

 

Clinical tangential curvature maps of the anterior corneal surface of (A) the less affected right eye and (B) the more affected left eye of a patient with asymmetric KC before simulation of progression in the right eye. Tangential curvature maps of the anterior corneal surface of (C) a central cone (region of maximum curvature) in a second patient and (D) a more eccentric cone in a third patient before simulation of the CXL effects.

Table 1.

 

Clinical Tomographic and Topographic Features of the Modeled Eyes

	Model Eyes for KC Progression		Central Cone	Eccentric Cone
	Right Eye	Left Eye
SimK, D	42.79/43.38@156°	42.70/43.48@101°	50.99/56.47@117°	44.94/49.13@30°
K max, D	43.76	47.10	66.24	58.13
CCT, μm	563	536	420	500
Location of thinnest point (x, y), mm	−0.41, −0.41	0.54, −1.02	−0.50, −0.57	1.08, −1.08
Corneal volume, mm3	65.3	63.0	57.8	56.3

Simulation of CXL in KC

The location of the cone in KC may influence whether predominantly central or peripheral flattening of the cornea occurs after CXL. ¹⁵ Therefore, two morphologically different KC eyes—one with a relatively central cone (Fig. 2C) and another with a more eccentric cone (Fig. 2D)—were used to construct two whole-eye, three-dimensional models to simulate the effects of CXL. The two rightmost columns in Table 1 show the SimK, K _max, CCT, and minimum thickness coordinates of the two corneas. The diameter of the weak zone was estimated to be 7 mm, according to the rationale presented earlier, and with α = 0.5 and dr = 0.075 at r = r ₀ (equation 1). With these conditions, a pair of values for m ₁ and m ₂ were calculated such that the volume-averaged reduction in elastic modulus (1 − α) in the weak zone was 0.5, in keeping with ex vivo evidence of a 50% to 60% decrease in strain-dependent elastic modulus in KC corneas. ^3,4 The weak zone was centered at the location of peak tangential curvature, the steepest point on the anterior surface. Next, the 0 load configuration of the model for each eye was calculated at an IOP of 15 mm Hg with the corneal elastic property, as given by equation 2. In the region outside the weak KC zone, α was 1 (indicating no decrease in elastic modulus). 

Clinical UV sources for CXL exhibit an intensity decrease of 10% along the radius in the standard 9-mm treatment diameter.³² Further, attenuation of UV intensity with increasing depth within the cornea has been experimentally demonstrated to follow the Beer-Lambert law.^33,34 Thus, incident UV intensity at any zone within the cornea is a function of radius and depth, assuming that the beam intensity is radially symmetric (independent of meridian) and the beam is aimed at the geometric center of the cornea. To model CXL, we assumed that the magnitude of increase in elastic modulus is linearly proportional to the beam intensity and that the peak magnitude of elastic modulus increase is at the geometric center of the cornea. Therefore in the weak zone, the factor α was modified to α′ such that  

 

where f (f > 1) is the magnitude of increase in elastic modulus due to CXL at the anterior-most cross-linked portions of the cornea, and (1 − g) is the decrease in the magnitude of elastic modulus increase f due to beam attenuation in deeper layers of the cornea. In equation 3, g < 1. In the region outside the weak zone and within the treatment zone  

 

The value of f was varied to account for the decline in beam intensity away from the beam center up to a treatment diameter of 9 mm. The value of g at any depth within the cornea was estimated by using an exponential decay function similar to the Beer-Lambert law. Two maximum depths of CXL-induced stiffening were simulated in each eye model—one at 200 and one at 300 μm—to calculate g at a given depth. Henceforth, this simulated CXL procedure will be referred to as the standard broad-beam CXL protocol: a 9-mm diameter treatment centered on the topographic map and with a spatial stiffening profile that accounts for loss of UV intensity as a function of radial distance and corneal depth. 

As an example, assume that the increase in elastic modulus, f, is 2 (a 200% increase) at the center of the topographic map and that the stiffening effect is limited to the anterior 200 μm of the cornea. In calculating g, it is also assumed that the concentration of the CXL agent is uniform up to the maximum depth within the cornea. Thus, g was modeled as an exponential decay function caused by beam attenuation only. In the anterior-most cross-linked portions of the cornea, α′ + (2 × 1 + 1) × α = 3 × α, where g is 1, which implies that in the anteriormost cross-linked portions of the cornea, the post-CXL elastic modulus is three times the pre-CXL elastic modulus at the center of the cornea. At 200 μm depth, g = 0 which implies that α′ = α, which in turn implies that at corneal depths greater than or equal to 200 μm, the post-CXL elastic modulus is unchanged from the pre-CXL state. In summary, the post-CXL elastic modulus was  

 

The zero-load model of each eye was then stressed with an IOP of 15 mm Hg with the post-CXL material property values specified. For each eye model (central and eccentric cone), two values of f at the center, 2 and 3, representing a maximum of 200% and 300% increase in elastic modulus were simulated and compared. Since a stiffer cornea would be more resistant to IOP-induced curvature changes, the IOP sensitivity of corneal curvature after CXL was analyzed by also simulating an IOP of 18 mm Hg with a 300% increase in elastic modulus in both central and eccentric cone CXL models. 

Alternate CXL Methods: Variable-Intensity Beam and Topography-Guided Focal CXL

Although UV sources are designed to maximize the uniformity of the beam, intensity drops approximately 10% toward the edges of the beam. ³² Even with some decrease in intensity, this pattern of irradiation may induce a steplike spatial transition in elastic modulus at the transition zone. For example, if at a 9-mm diameter the intensity of a UV source suddenly drops from 100% to 0%, a similar step change in elastic modulus at the 9-mm diameter could occur. To explore an alternative approach to irradiation and potential differences in the response to CXL, we simulate a UV beam with a graded change in intensity such that the slope of intensity change is continuous at the edge of the UV beam. We measured the intensity of the UV beam of a prototype UV-A source (PriaVision Inc., Menlo Park, CA) using a digital SLR camera and analyzed the blue channel scalar values (MatLab, ver. 2010; Mathworks, Inc., Natick, MA). All images were taken normal to the plane of the beam. The measured data are shown in Figure 3 as normalized intensity versus unit diameter. The beam diameter was 9 mm. The spatial–average intensity of the beam was 3 mW/cm². A significant variation in the intensity as a function of beam radius was observed. We chose to model these measurements as a representative form of a variable-intensity beam, with a smoother intensity transition near the edge of the beam compared with the standard broad-beam protocol. 

Figure 3.

 

Experimentally measured UV intensity normalized with its peak value at the center versus the unit diameter measured from a prototype UV-A source. This profile was used in simulations of variable-intensity broad-beam and cone-localized, variable-intensity treatment patterns.

To implement this form of intensity variation, f was also modeled as a function of radius similar to the mathematical function, α. The function returned the value of f such that an exponential decay away from the center of the beam was modeled with a maximum value of f at the center of the incident beam and f = 0 at the edge of the beam. The exponential decay function for f was obtained by curve-fitting the normalized intensity data from the UV-A source. Then, the normalized intensity values were linearly scaled to the desired magnitude of increase in elastic modulus f. In other words, f was similar in shape to the UV intensity profile shown in Figure 3. Thus, similar to the standard broad-beam protocol, the post-CXL elastic property of the cornea within the treatment zone was given by equation 5, where α′ was described by equations 3 and 4, within and outside the weak zone, respectively. For each eye model (central and eccentric cones), two values of f at the center, 2 and 3, representing a maximum of 200% and 300% increase in elastic modulus at two maximum CXL depths of 200 and 300 μm, were simulated. The simulated treatment diameter was 9 mm. The remaining modeling features were the same as the standard broad-beam protocol. Henceforth, these simulated CXL procedures are referred to as the variable-intensity broad-beam protocol. 

To simulate a variable-intensity and cone location–specific treatment, we simulated a smaller diameter UV beam equal to the diameter of the weak zone and with the beam centered on the weak zone. Similar to the variable-intensity broad-beam protocol, f was represented by the mathematical function α, such that the maximum increase in stiffness was simulated at the center of the weak zone (the locus of maximum tangential curvature), and no increase in stiffness was simulated beyond the weak zone. Thus, within the weak zone, the post-CXL elastic property was modeled by equation 5. Outside the weak zone, the post-CXL elastic property was modeled by equation 2. For each eye model (central and eccentric cone), two values of f at the center, 2 and 3, representing a maximum of 200% and 300% increase in elastic modulus at two limiting CXL depths of 200 and 300 μm were simulated. Henceforth, these simulated CXL procedures are referred to as the cone-localized, variable-intensity protocol. 

Analysis of Tangential Curvature and Wavefront Aberrations of Anterior Corneal Surface

The x, y, and z coordinates of the anterior surface of the cornea were fitted to 10th-order Zernike polynomials with a normalization radius of 5.75 mm. The local curvatures of the anterior surface cornea along meridians from 0° to 355° at 5° intervals and along the radii at 0.05-mm intervals. Tangential curvature was calculated as  

 

where R _t is the tangential radius of curvature in meters, and n = 1.3375 is the keratometric refractive index of the cornea. R _t was calculated using the first and second derivatives of Zernike polynomials fitted to the nodal co-ordinates of anterior corneal FEM surface. Maximum axial curvature is referred to as K _max. SimK, calculated at the central 3 mm of diameter, is given by flattest/steepest meridian@axis. The wavefront aberrations of the anterior surface of the cornea with a pupil diameter of 5 mm were calculated by using Zernike polynomials. The lower order root mean square (LO RMS) was calculated by using the second-order Zernike coefficients. The higher order root mean square (HO RMS) was calculated by using the 3rd- to 10th-order Zernike coefficients. We also analyzed the change induced by different simulated CXL procedures on independent Zernike terms, namely, defocus (C _2,0), coma (C _3,−1,C _3,1), and spherical aberration (C _4,0). The effects of KC progression and of each CXL protocol on the surface area of the anterior cornea from the FEM results were also analyzed across a central 5-mm-diameter zone. 

Results

Simulated Progression of KC

Figures 4A–C show the tangential curvature maps of the anterior corneal surface (Fig. 2A: right eye) with focal 10%, 30%, and 45% reductions in the elastic modulus centered on the steepest anterior corneal point. Figures 4D–F show the corresponding difference in tangential curvature between the weakened cornea and the measured in vivo state (Fig. 2A). With further weakening of the cornea, the tangential curvature increased, producing a steep cone surrounded by an annular zone of progressive flattening (decreased tangential curvature). Two distinct regions were observed in all the difference maps: a region of steepening in the weak zone and a region of flattening surrounding the weak zone. A comparison of SimK values calculated at a diameter of 3 mm before and after a 45% average reduction in elastic modulus (Table 2, case 1) demonstrated increases of 7.30 and 5.59 D along the steepest and flattest meridian, respectively. With the same 45% reduction in elastic modulus, K _max increased by 12.23 D in the right eye. K _max in the more clinically affected left eye (Fig. 2B) was 47.1 D. If similar predisease geometries are assumed, a 32% difference (Fig. 5) in the initial elastic moduli of the two eyes can be estimated from the decrease in elastic modulus required to produce a similar K _max. In Figure 5, K _max is seen to increase exponentially as a function of decreasing elastic modulus. 

Figure 4.

 

FEM-generated tangential curvature maps of the anterior corneal surface of the less affected eye of an asymmetric KC patient after elastic modulus reductions of (A) 10%, (B) 30%, and (C) 45%. (D–F) Associated tangential curvature difference maps for each elastic modulus decrement.

Table 2.

 

Simulation of KC Progression

Reduction in Elastic Modulus (%)	SimK (D)	K max (D)	LO RMS (μm)	HO RMS (μm)	DefocusComaSpherical				Surface Area (mm2)
					C2,0	C3,−1	C3,1	C4,0
Case 1: Reduced Thickness by 27 μm and Reduced Elastic Modulus
Right Eye	42.79/43.38@156°	43.76	0.69	0.29	0.457	−0.111	−0.129	0.128	89.186
10	43.08/43.71@156°	44.22	0.63	0.32	0.337	−0.165	−0.161	0.100	89.228
30	44.36/45.24@159°	46.50	0.75	0.64	−0.305	−0.500	−0.338	−0.042	89.369
45	48.35/50.68@159°	55.99	3.51	2.57	−3.023	−2.132	−1.129	−0.647	89.705
Case 2: Reduced Elastic Modulus Only
10	43.10/43.61@162°	44.01	0.69	0.34	0.826	−0.122	−0.157	0.177	89.227
30	44.27/45.05@162°	46.02	0.73	0.58	0.273	−0.430	−0.313	0.052	89.365
45	47.91/50.00@162°	54.38	2.64	2.23	−2.038	−1.869	−0.998	−0.468	89.684

 Changes in topographic and optical characteristics such as SimK, K _max, lower and higher order aberrations, and surface area of the anterior corneal surface with progressive hyperelastic weakening of the cornea.

Figure 5.

 

K _max as a function of the decrease in hyperelastic modulus in a simulation of KC progression in an eye with subclinical KC.

The LO RMS and HO RMS of the right eye at various stages of disease progression are presented in Table 2. Both LO and HO RMS increased with progressive simulated weakening of the cornea. With a decrease in average elastic modulus by 45% in the weak zone, both LO RMS and HO RMS increased significantly to 3.51 and 2.57 μm, respectively. The initial decrease in LO RMS in the right eye (Table 2) with a 10% reduction in corneal elastic modulus is most likely attributable to redistribution of the more diffusely irregular original curvature distribution (Fig. 2A) to a more paraxial steep pattern (Figs. 3B, 3C). Further weakening of the weak zone resulted in significant increases in defocus and astigmatism. Table 2 also lists the changes in the individual Zernike terms, with all terms becoming progressively more negative as a function of corneal weakening. The area of the anterior surface increased with progression of the disease, although the magnitude of increase was very small. 

Table 2 also presents the predicted optical characteristics of KC progression in the right eye with the same reductions in elastic modulus as case 1, but with no thickness change from the in vivo state (case 2). Between cases 1 and 2, the difference in CCT was 27 μm. Incremental reductions in the elastic modulus without a reduction in thickness resulted in progressive steepening and increases in aberrations and surface area qualitatively similar to those observed in case 1. However, compared to case 1, the magnitudes of observed increases were lower (e.g., LO RMS: 2.64 vs. 3.51 at 45% reduction), which suggests some additive effect of reduced thickness and reduced elastic modulus on steepening and aberrations. 

Modeling of CXL of a Central Cone

Figures 6A–C are tangential curvature maps of a cornea undergoing simulated CXL with the standard broad-beam, variable-intensity broad-beam, and cone-localized, variable-intensity protocols, respectively, with a 300% maximum increase in elastic modulus limited to the anterior 200 μm. Figures 6D–F show the corresponding changes in tangential curvature between the CXL state and the pre-CXL state. (The latter is presented in Fig. 2C). In all the tangential curvature difference maps (Figs. 4D–F), the area within and around the cone flattened, and an annular region around the cone steepened. With a change in the treatment protocol from standard broad-beam to cone-localized, variable-intensity, the flattening effect became more focal and was accompanied by a discrete ring of steepening around the cone (Fig. 4F). With the standard treatment, this ring of steepening was less continuous and superior to the cone (Fig. 4D). K _max decreased by 2.61, 4.36, and 5.62 D with the standard broad-beam, variable-intensity broad-beam, and cone-localized variable-intensity protocols, respectively (Table 3). Focal treatment centered on the cone with spatial variation of UV-A intensity therefore yielded the greatest flattening effect. Whereas LO RMS was relatively insensitive to the magnitude of the elastic modulus increase within the range explored here, HO RMS, defocus, coma, and spherical aberration all decreased significantly with the incremental increase in stiffness (Table 3). These stiffness-dependent decreases in aberrations, such as K _max, were at their maximum in the cone-localized, variable-intensity protocol. The surface area of the anterior cornea also decreased with increased stiffening of the cornea, although as before, the magnitude of change was very small. No differences between the post-CXL cone locations were found when the results of the three simulated CXL protocols were compared to each other. 

Figure 6.

 

Tangential curvature maps of a cornea with a central cone after simulated CXL to a maximum depth of 200 μm in (A) a standard broad-beam stiffening protocol, (B) a variable-intensity broad-beam protocol, and (C) a cone-localized, variable-intensity protocol. (D–F) The corresponding tangential curvature difference maps.

Table 3.

 

CXL of a Cornea with a Central Cone

Protocol	% Increase in Elastic Modulus	SimK (D)	K max (D)	LORMS (μm)	HORMS (μm)	DefocusComaSpherical				Surface Area (mm2)
						C2,0	C3,−1	C3,1	C4,0
200 μm
Pre-CXL	—	50.99/56.47@117°	66.24	4.98	3.67	−3.116	−3.215	−1.232	−0.722	90.627
Standard broad-beam	200	49.06/55.33@114°	64.30	4.89	3.40	−1.936	−2.904	−1.192	−0.493	90.418
	300	48.42/54.98@114°	63.63	4.93	3.30	−1.555	−2.790	−1.181	−0.421	90.344
Variable-intensity broad-beam	200	49.39/54.39@114°	62.96	4.61	3.31	−1.651	−2.850	−1.164	−0.414	90.453
	300	46.50/53.73@114°	61.88	4.61	3.19	−1.222	−2.708	−1.138	−0.325	90.389
Cone-localized, variable-intensity	200	48.85/54.61@114°	62.11	4.26	3.04	−1.055	−2.557	−1.151	−0.312	90.566
	300	48.06/53.94@114°	60.62	4.24	2.83	−0.309	−2.321	−1.122	−0.164	90.543
Standard broad-beam	300% with IOP = 18 mm Hg	49.88/55.76@117°	65.19	4.90	3.53	−2.516	−3.059	−1.201	−0.616	90.347
300 μm
Standard broad-beam	200%	48.56/55.03@114°	63.74	4.92	3.32	−1.648	−2.811	−1.183	−0.439	90.353
	300%	48.11/54.63@114°	62.99	4.94	3.22	−1.183	−2.684	−1.172	−0.352	90.264
Variable-intensity broad-beam	200%	47.68/53.87@114°	62.10	4.61	3.21	−1.293	−2.741	−1.145	−0.341	90.394
	300%	46.64/53.11@114°	60.87	4.70	3.08	−0.766	−2.585	−1.121	−0.233	90.313
Cone-localized, variable-intensity	200%	48.29/54.14@114°	61.13	4.25	2.90	−0.578	−2.403	−1.128	−0.216	90.546
	300%	47.38/53.38@114°	59.46	4.33	2.68	0.255	−2.141	−1.095	−0.051	90.518

 Changes in topographic and optical characteristics as a function of simulated treatment depth, spatial stiffening profile, and treatment localization.

Figures 7A–C are tangential curvature maps of a central cone undergoing simulated CXL with the same three protocols and a 300% maximum increase in elastic modulus, but localized to the anterior 300 μm. Figures 7D–F show the corresponding changes in tangential curvature between the CXL state and the pre-CXL state. Although the trends in SimK, K _max, aberration metrics, and surface area were similar to those observed with the 200-μm stiffening depth, greater magnitudes of all changes were observed with a deeper simulated stiffening effect. K _max decreased by 3.25, 5.37, and 6.78 D with the standard broad-beam, variable-intensity broad-beam, and cone-localized, variable-intensity protocols, respectively (Table 3). Even with the increased depth of stiffening, the cone-localized, variable-intensity protocol yielded the greatest flattening effect and aberration decreases of the three protocols. The effect of a modest IOP increase on model geometry was also analyzed in the post-CXL model using the standard broad-beam protocol with a treatment depth of 200 μm and a 300% increase in hyperelastic modulus (Table 3). An IOP increase from 15 to 18 mm Hg reduced the flattening effect from 2.61 to 1.05 D as measured by the change in K _max. Aberrations (LO RMS, HO RMS, and individual terms) were lower after CXL at both pressures compared to the pre-CXL state, but were higher at an IOP of 18 mm Hg than at 15 mm Hg. No change in the pretreatment cone location was observed with any of the simulated CXL protocols. 

Figure 7.

 

Tangential curvature maps of a cornea with a central cone after simulated CXL to a maximum depth of 300 μm in (A) a standard broad-beam stiffening protocol, (B) a variable-intensity broad-beam protocol, and (C) a cone-localized, variable-intensity protocol. (D–F) The corresponding tangential curvature difference maps.

Modeling of CXL of an Eccentric Cone

Figures 8A–C are tangential curvature maps of an eccentric cone undergoing simulated CXL with the standard broad-beam, variable-intensity broad-beam, and cone-localized, variable-intensity protocols, respectively, with a 300% maximum increase in elastic modulus limited to the anterior 200 μm. Figures 8D–F show the corresponding changes in tangential curvature from the pre- and post-CXL states (The former is shown in Fig. 2D). As in the case of a central cone, the eccentric cone flattened after CXL (Figs. 8D–F). However, with the standard broad-beam protocol, the difference map is marked by a notable flattening slightly peripheral to the location of the inferior cone and compensatory steepening along the opposite hemimeridian superior to the cone (Fig. 8D). This differs from the annular ring pattern seen in the central cone with the standard broad-beam CXL protocol (Figs. 6D, 7D). With the variable-intensity broad-beam protocol, the flattening occurred at the location of the cone and was surrounded by steepening in nearly all other areas of the anterior corneal surface. This topography could have contributed to the increase in spherical aberration observed with this CXL protocol (Table 4). The pre-CXL K _max (58.13 D) was reduced to a lower magnitude with the variable-intensity broad-beam protocol (55.86 D, a 2.27-D decrease) than with the standard broad-beam protocol (57.16 D, a 0.97-D decrease). With the cone-localized, variable-intensity protocol, the magnitude of flattening was greatest (55.66 D, a 2.47-D decrease with f = 3; Table 4), and the zone of maximum flattening was localized to the cone and surrounded by a distinct annular ring of steepening. Although the CCT differed by 80 μm and the location of minimum thickness differed in the central KC and eccentric KC cases, the two model eyes had similar corneal volumes of 57.8 and 56.3 mm³ (Table 1). 

Figure 8.

 

Tangential curvature maps of a cornea with a more eccentric cone after simulated CXL to a maximum depth of 200 μm in (A) a standard broad-beam stiffening protocol, (B) a variable-intensity broad-beam protocol, and (C) a cone-localized, variable-intensity protocol. (D–F) The corresponding tangential curvature difference maps.

Table 4.

 

CXL of a Cornea with an Eccentric Cone

Protocol	% Increase in Elastic Modulus	SimK (D)	K max (D)	LORMS (μm)	HORMS (μm)	DefocusComaSpherical				Surface Area (mm2)
						C2,0	C3,−1	C3,1	C4,0
200 μm
Pre-CXL		44.94/49.13@30°	58.13	3.25	3.54	1.642	−2.822	1.989	0.349	91.104
Standard broad-beam	200	45.35/48.58@27°	57.39	2.74	3.34	1.671	−2.679	1.881	0.349	90.928
	300	45.50/48.41@27°	57.16	2.55	3.28	1.664	−2.632	1.845	0.345	90.861
Variable-intensity broad-beam	200	44.36/47.99@30°	56.47	2.90	3.17	1.667	−2.541	1.773	0.374	90.976
	300	44.26/47.60@27°	55.86	2.74	3.04	1.634	−2.439	1.696	0.374	90.925
Cone-localized, variable-intensity	200	45.75/48.74@30°	56.33	2.06	3.02	0.871	−2.449	1.708	0.219	91.062
	300	46.09/48.59@27°	55.66	1.63	2.84	0.601	−2.310	1.604	0.172	91.046
Standard broad-beam	300% with IOP = 18mmHg	44.99/48.72@30°	57.84	3.15	3.51	1.905	−2.798	1.975	0.379	90.839
300 μm
Standard broad-beam	200	45.45/48.44@27°	57.14	2.58	3.28	1.631	−2.634	1.847	0.340	90.868
	300	45.64/48.26@27°	56.86	2.35	3.20	1.597	−2.576	1.803	0.330	90.790
Variable-intensity broad-beam	200	44.29/47.68@27°	55.97	2.77	3.07	1.633	−2.460	1.712	0.372	90.929
	300	44.14/47.26@27°	55.25	2.56	2.91	1.564	−2.340	1.622	0.365	90.865
Cone-localized, variable-intensity	200	45.98/48.60@27°	55.85	1.76	2.89	0.706	−2.347	1.633	0.191	91.047
	300	46.30/48.43@27°	55.07	1.29	2.69	0.408	−2.186	1.512	0.138	91.026

 Changes in topographic and optical characteristics as a function of simulated treatment depth, spatial stiffening profile, and treatment localization.

With all three CXL protocols, significant decreases in aggregate aberration RMS values as well as individual aberrations occurred (Table 4). As in the case of the central cone, the cone-localized, variable-intensity protocol produced the greatest decreases in aberrations in the eccentric cone. The decrease in surface area (Table 4) was maximum with the cone-localized stiffening protocol. Similar to CXL of the central cone with standard broad-beam protocol, an IOP increase partially counteracted the flattening effect of CXL (Table 4). However, defocus and spherical aberration increased with IOP increase (Table 4), possibly due to the differential location of the cones in the two eye models. Another effect of CXL on the eccentric cone was that the cone location changed with stiffening. With the standard broad-beam protocol, the cone was located 1.9 mm from the central cornea along the 306° radius in the pre-CXL state, 1.825 mm @ 306° in the post-CXL state with a 200% increase in modulus, and 1.775 mm @ 306° in the post-CXL state with a 300% increase in modulus. Thus, the cone shifted toward the corneal center along the 306° radius as a function of increasing CXL effect from a standard broad-beam treatment simulation. By contrast, with the variable-intensity broad-beam protocol and the same pre-CXL cone location, the post-CXL cones were located at 1.95 mm @ 306°, with a 200% increase in modulus, and 1.975 mm @ 306°, with a 300% increase, indicating an outward radial shift of the cone. With the cone-localized, variable-intensity protocol, the cone location did not change. 

Figures 9A–C are tangential curvature maps of the CXL cornea with the standard broad-beam, variable-intensity broad-beam, and cone-localized, variable-intensity protocols respectively, with the increase in elastic modulus of 300% localized to the anterior 300 μm. Figures 9D–F show the corresponding differences in tangential curvature between the pre- and post-CXL states. With the increased depth of simulated treatment, the trends (SimK, K _max, aberrations, and surface area) were similar to the 200-μm depth results, except that the magnitudes of the changes were higher. K _max decreased by 1.27, 2.88, and 3.06 D with the standard broad-beam, variable-intensity broad-beam, and cone-localized, variable-intensity protocols, respectively (Table 4). The cone-localized, variable-intensity protocol yielded greater reductions in K _max (to 55.07 D with the 200-μm depth and 55.66 D with the 300-μm treatment) and aberrations than the other two protocols. The CXL-mediated decreases in aberrations and K _max were lower when IOP was increased to 18 mm Hg (Table 4). Similar trends in the movement of the cone location were also found with increased depth of stiffening. With the standard broad-beam protocol, the cone location was 1.9 mm @ 306°, 1.8 mm @ 306°, and 1.725 mm @ 306° in the pre-CXL state, the post-CXL state with 200% increase in modulus, and the post-CXL state with 300% increase in modulus, respectively. With the variable-intensity broad-beam protocol, the post-CXL cone locations were at 2.0 mm @ 306° and 2.05 mm @ 306°, respectively. As before, with the cone-localized, variable-intensity protocol, the cone location did not change. 

Figure 9.

 

Tangential curvature maps of a cornea with an eccentric cone after simulated CXL to a maximum depth of 300 μm in (A) a standard broad-beam stiffening protocol, (B) a variable-intensity broad-beam protocol, and (C) a cone-localized, variable-intensity protocol. (D–F) The corresponding tangential curvature difference maps.

Discussion

Through these patient-specific, three-dimensional, whole-eye, FEM experiments, we showed that focal mechanical weakening and thinning of the cornea produce topographic characteristics consistent with those in clinical KC. The simulations predicted changes in curvature with progressive elastic weakening of the mildly affected eye that closely mimicked the more advanced disease of the contralateral eye. Lower and higher order aberrations increased with progression of disease in a manner consistent with previous clinical studies. ^{35 –37} We also found that K _max increased sharply as the reduction in elastic modulus progressed from 30% to 45% (Fig. 5), whereas a generalized thickness reduction produced no significant additional increase in K _max (Table 2). The nonlinear nature of the relationship shown in Figure 5 allows for a disease progression model where microstructural events driving a linear decrease in hyperelastic properties could manifest as an accelerating topographic decompensation in later stages of the disease. 

In a previous FEM study, ¹⁹ the investigators concluded that thinning, one of the clinical hallmarks of KC, is the most influential factor in the progression of KC. However, these conclusions were based on extreme reductions in thickness (minimum thicknesses of 200 μm) much more representative of end-stage disease than early KC. The authors also concluded that thinning alone was insufficient to cause marked degradation in every modeled case. The relative sensitivity of different models to geometric and mechanical variables is likely to depend on differences in model assumptions, differences in the range of geometric and mechanical abnormalities modeled, and the specific case geometries modeled. Although geometry and constitutive mechanical properties can be modeled as independent features of a structure and can thus vary independently, in the clinical setting, thinning and mechanical property changes are apt to be closely linked. However, the current simulations, as well as clinical studies demonstrating the absence of net increases in corneal thickness after CXL, stand as evidence that corneal elastic properties can change, independent of corneal thickness in the clinical domain. ³⁸ We believe that more extensive sensitivity analyses should be performed to explore the relative contributions of thickness reduction and reduced elastic modulus for different patient-specific initial geometries. Such simulations could also explore the differential effect of a uniform reduction in thickness versus more focal thickness reductions. The observation that the topographic and optical features of KC progression were generated without any requirement for corneal thinning suggests that a progressive reduction in elastic modulus is a critical driving force in the early development of the topographic phenotype of KC. 

Very small increases (<0.6%) in anterior corneal surface area were measured during simulated progression of KC. This relative conservation of surface area is consistent with the clinical analysis of Smolek and Klyce, ³⁹ demonstrating no statistically significant placido-derived surface area differences between keratoconic, control, and post–refractive surgery patients and no measurable change in surface area in one patient who exhibited clear topographic progression over 6 years of follow-up. Other clinical studies have shown no differences in corneal cross-sectional mass between KC and control patients ⁴⁰ and no tomographic corneal volume differences in any other than the most severe group of KC patients. ⁴¹ Our finding that the topographic phenotype of mild to moderate KC can be replicated in simulations without thinning or significant increases in surface area is consistent with this body of evidence and supports the hypothesis that abnormal bending or warpage rather than true ectasia (distention or elongation) of tissue predominates the topographic phenotype. ³⁹ If corneal thinning is not a prerequisite for KC progression and is a later or less consistent feature of KC than material property abnormalities, then methods for directly measuring biomechanical abnormalities may offer important sensitivity and specific advantages for screening and early intervention. 

The precise nature of the material property deficits in KC are not known and remain an area of great interest as well as a potential source of error in modeling experiments. In modeling the spatial mechanopathology of KC and its progression, we started with evidence from formal extensometry experiments that the KC cornea has a lower elastic modulus than do normal corneas. ^3,4 The same experiments describe a percentage reduction in whole-specimen corneal elastic modulus that was applied to the model for simulating KC progression. It should be noted that the cited experiments were not capable of differentiating the spatial distribution of the elastic properties; therefore, a 55% reduction in elastic modulus across an entire penetrating keratoplasty specimen may represent an aggregate measure of higher and lower local moduli across the specimen. In the absence of any available peer-reviewed literature on the actual spatial profile of corneal elastic property reductions in KC, an assumption was made that the biomechanical abnormalities are most profound in the area of the topographic cone, an approximation that is consistent with observations of more concentrated ultrastructural and histologic abnormalities in the region of the cone, including fragmentation of Bowman's layer, ⁴² less lamellar interweaving, ⁹ and fewer lamellar insertions into Bowman's layer. ⁴³ The assumption that property variations occur in nature as a continuous gradation of properties rather than a discrete step-change is consistent with biological systems in general and is also consistent with the histopathological features of KC. 

We simulated CXL procedures that incorporated the physical effects of radial variation and axial attenuation of the intensity of typical clinical UV-A sources. We also explored the performance of alternate CXL protocols with combinations of a smoother peripheral spatial intensity transition, a deeper CXL effect, and a cone-centered treatment approach. The optical response to CXL as measured by curvature and first-surface aberrations was sensitive to each of these variables. Simulation of a deeper CXL effect resulted in enhanced flattening of the cone in all cases modeled in this study. The depth of effect may be influenced by the presence and composition of the precorneal riboflavin film, ⁴³ the presence of the corneal epithelium, ^14,44 and the integrity of the epithelial tight junctions. ⁴⁵ Attempts at deeper CXL with the current UV approach could be limited by anterior absorption of UV energy by riboflavin-saturated stroma and potential risks to the keratocytes, limbal stem cells, endothelium, and deeper intraocular structures that are attributable to intensified UV exposures. Our simulations also suggest that greater flattening of the cone can be achieved with a reduction in the effective treatment diameter, more gradual intensity transition zones, and centration of treatment on the cone. An additional theoretical benefit of these approaches is a reduction in the total energy delivery requirement. These simulations suggest that novel stiffening treatments specifically targeting areas of lowest material strength may produce more optimal shape changes than a standard broad-zone treatment approach. It is possible that in certain forms of transepithelial cross-linking, an epithelium that is consistently thinner over the cone could favor an intrinsic localization of the CXL treatment by maximizing penetration and exposure in the most severely affected regions of the stroma. 

Another key finding of the study was the sensitivity of the CXL outcome to cone location. The eccentric cone centralized slightly with the standard CXL simulation, decentered slightly with a smaller effective treatment zone, and remained in its original location with the cone-centered, variable-intensity treatment. The dependence of standard broad-zone CXL results on cone location are consistent with the clinical trend, ¹⁵ but the latter two observations are novel and point to the importance of the stiffening distribution relative to cone location. Specifically, a smaller zone of stiffening that is centered on the cornea and does not fully encompass an eccentric weak zone differentially stiffens the central cornea, reduces the arc length of central corneal collagen, and pushes the bending maximum (K _max) away from the central cornea. Conversely, a broader central zone of treatment encompasses more of the paracentral cone and favors a central shift by normalizing the ratio of material properties between the cone and the less affected areas. Finally, a cone-centered treatment is neutral with respect to impact on final cone location because the stiffening effect is centered on and geometrically matched to the zone of weakness. In addition to illustrating the importance of considering cone location in future CXL treatment designs, these results suggest the need for a more nuanced approach to measuring the outcomes of CXL procedures in clinical trials. Specifically, spatially fixed clinical outcome measures such as standard keratometry, simulated K values, or central corneal curvature could underestimate the effect of CXL or even suggest paradoxical steepening in cases of centralization of eccentric cones. Furthermore, axial curvature algorithms are rooted in paraxial optics assumptions and provide less accurate representations of the shape and location of steep topographic features, ⁴⁶ particularly outside of the central cornea. Although measures with a central bias do provide important information about the optical impact of topographic changes across the pupillary zone, metrics of cone severity such as the cone location and magnitude index (CLMI) ⁴⁷ may provide a more complete characterization of the topographic effect of treatment. 

CXL outcomes were sensitive to another intrinsic patient factor, IOP, as summarized in Tables 3 and 4. This sensitivity was found to a lesser degree with the more effective treatment protocols. Even a modest 3-mm Hg increase in IOP reduced the topographic effect of CXL, increased surface area nominally, and mimicked the topographic appearance of KC progression. This model behavior is consistent with clinical observations of topographic steepening with forced IOP elevation in KC and potential reductions in corneal axial curvature with IOP reduction. ⁴⁸ Similar reduction in curvature was observed after IOP reduction in post-LASIK ectasia eyes, which are also biomechanically weaker than normal. ⁴⁹ Medication-induced reduction of IOP has been suggested as a treatment for KC. The present study supports the theoretical utility of IOP reduction for at least a temporary reduction of topographic severity. However, data on the ability of IOP reduction to reduce the long-term risk of progression are not available, and any long-term disease-modifying effect would invoke viscoelastic behavior and long-term modifications of mechanobiology (i.e., interruption of a cycle of progressive hyperelastic weakening and thinning) that are beyond the scope of this study. This study does suggest that short-term reduction of IOP during and for some period after CXL until a stiffening effect is achieved could be of benefit. 

The quantitative results of these computational simulations compare favorably with those of published clinical studies of CXL for KC with long-term follow-up. ^{50 –53} The cited studies did not analyze outcomes as a function of cone location, size of treated zone, or depth of effect, but have reported mean decreases in K _max between 2.5 and 3 D. ^{50 –53} The model results for the case of a central cone showed that corneal stiffness increases in the range of 200% to 300% would be necessary to achieve this degree of flattening. These estimates are comparable to measurements of the stiffening caused by UV-A/riboflavin CXL in ex vivo corneas. ¹⁶ Case-specific estimates of material property changes can be obtained through an inverse three-dimensional FEM process, demonstrated previously by our group in laser in situ keratomileusis (LASIK) patients, ²⁹ and is currently being explored as a method for characterizing the material property changes typically produced by a specific CXL treatment. The FEMs described herein required only a few minutes of computational time to reach convergence, which lends to viability in large-scale clinical modeling. 

Whether applied prospectively as a predictive tool or retrospectively as a tool for estimating property changes after treatment, the modeling process has limitations related to the accuracy of its assumptions. Although the curvature maximum is a reasonable surrogate for the weakest point of the cornea, the exact dimensions and distribution of weakness are not known and have been assumed to be circular and radially symmetric as a first approximation. Our models are derived from clinical topography and tomography maps based on anterior and posterior corneal surface tracings and therefore take into account the aggregate clinical shape of the corneal stroma and epithelium. However, the epithelium was not modeled as a separate entity due to its minimal contribution to corneal elastic behavior relative to the biomechanically predominant stroma. ⁵⁴ Characteristic nonuniformities in epithelial thickness such as those demonstrated in normal eyes ⁵⁵ and KC eyes ⁵⁶ can mask underlying stromal curvature and introduce a source of modeling error. Because the keratoconic epithelium is typically thinnest over the cone, ⁵⁶ epithelial debridement during CXL may unroof a steeper stroma, which could explain the appearance of a paradoxical steepening response in the first months ⁵³ after CXL that is followed by progressive flattening. Our post-CXL model assumes that the epithelium has healed with the same thickness distribution that was present in the pre-CXL state. A recent study found a slight difference in the epithelial thickness profile before and after CXL in post-LASIK ectasia, ⁵⁷ but similar data in KC patients after CXL are not yet available. More accurate modeling of KC progression and CXL outcomes may be possible if the geometric effects of the epithelium are considered in future studies. Future modeling efforts could also incorporate case-specific estimates of depth of treatment effect rather than idealized depths using the optical demarcation line visible in post-CXL corneas with optical coherence tomography (OCT). ⁵⁸ In this study, we modeled KC progression as a progressive reduction in the hyperelastic material properties without appealing to viscoelastic properties. However, viscous dissipation functions could be incorporated ⁵⁹ and investigated as another means of simulating progression or regression. Future studies also should address the impact of a gradient concentration of riboflavin within the cornea and corresponding material property gradients, which could be incorporated through the inclusion of another function in equation 3. 

In summary, this simulation-based study demonstrates a useful patient-specific methodology for evaluating the proposed role of corneal hyperelastic property abnormalities in the progression of KC. We also provide in silica proof of concept that focal treatment with a variable-intensity UV-A source may enhance the CXL effect and optimize corneal first-surface optics while reducing energy exposure relative to standard broad-beam approaches. The simulations also suggest that clinical tools for resolving spatial biomechanical properties are needed not only for better understanding the early pathogenesis of KC but also for better defining the therapeutic target of customized cross-linking treatment algorithms. Modifications to clinical protocols would require validation in clinical studies, and tissue-based experiments are under way in our laboratory to investigate these concepts. 

Footnotes

 Supported in part by National Institutes of Health Grants K12RR023264 and 1KL2RR024990; challenge and unrestricted grants from Research to Prevent Blindness to the Department of Ophthalmology of the Cleveland Clinic Lerner College of Medicine of Case Western Reserve University; and the National Keratoconus Foundation/Discovery Eye Foundation. WJD is a recipient of a Research to Prevent Blindness Career Development Award.

Footnotes

 Disclosure: A.S. Roy, P; W.J. Dupps, Jr, P