**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.

^{ 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.

^{ 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.

^{ 16 }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

^{2}.

^{ 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.

^{ 15 }

^{ 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.

^{ 21 }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.

^{ 21 }

^{ 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.

*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

^{−5}μ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.

^{ 21 }

^{ 24 –26 }Computational simulation also suggests that IOP measurement by applanation is sensitive to variations in corneal elastic modulus.

^{ 27 }Clinical studies have also suggested a small increase in IOP.

^{ 28 }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.

^{ 29 }The full Newton's method was used to solve the FEM equations.

^{ 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.

^{ 30 }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

^{ 16 }was fitted to a reduced polynomial material model,

*W*=

*C*

_{10}(

*I*

_{1}− 3) =

*C*

_{20}(

*I*

_{1}− 3)

^{2}, where

*W*is the strain energy potential,

*C*

_{10}and

*C*

_{20}are constants to be specified later, and

*I*

_{1}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.

^{ 31 }

^{ 18 }Since the cornea is weaker in this circular zone, the constants

*C*

_{10}and

*C*

_{20}for normal cornea were reduced, as a function of radius and independent of corneal thickness, by α equal to where

*m*

_{1}and

*m*

_{2}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*

_{1}and

*m*

_{2}is described below for each simulation.

*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}, α = 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).

^{ 10 }With these conditions, three pairs of values for

*m*

_{1}and

*m*

_{2}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*

_{10}and

*C*

_{20}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.

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, mm^{3} | 65.3 | 63.0 | 57.8 | 56.3 |

^{ 15 }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*

_{0}(equation 1). With these conditions, a pair of values for

*m*

_{1}and

*m*

_{2}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).

^{32}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.

*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.

^{ 32 }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

^{2}. 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.

*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.

*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.

*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.

*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.

Reduction in Elastic Modulus (%) | SimK (D) | K _{max} (D) | LO RMS (μm) | HO RMS (μm) | DefocusComaSpherical | Surface Area (mm^{2}) | |||
---|---|---|---|---|---|---|---|---|---|

C_{2,0} | C_{3,−1} | C_{3,1} | C_{4,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 |

*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.

Protocol | % Increase in Elastic Modulus | SimK (D) | K _{max} (D) | LORMS (μm) | HORMS (μm) | DefocusComaSpherical | Surface Area (mm^{2}) | |||
---|---|---|---|---|---|---|---|---|---|---|

C_{2,0} | C_{3,−1} | C_{3,1} | C_{4,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 |

*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.

*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

^{3}(Table 1).

Protocol | % Increase in Elastic Modulus | SimK (D) | K _{max} (D) | LORMS (μm) | HORMS (μm) | DefocusComaSpherical | Surface Area (mm^{2}) | |||
---|---|---|---|---|---|---|---|---|---|---|

C_{2,0} | C_{3,−1} | C_{3,1} | C_{4,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 |

*outward*radial shift of the cone. With the cone-localized, variable-intensity protocol, the cone location did not change.

*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.

^{ 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.

^{ 19 }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.

^{ 38 }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.

^{ 39 }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

^{ 40 }and no tomographic corneal volume differences in any other than the most severe group of KC patients.

^{ 41 }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.

^{ 39 }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.

^{ 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,

^{ 42 }less lamellar interweaving,

^{ 9 }and fewer lamellar insertions into Bowman's layer.

^{ 43 }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.

^{ 43 }the presence of the corneal epithelium,

^{ 14,44 }and the integrity of the epithelial tight junctions.

^{ 45 }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.

*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,

^{ 15 }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,

^{ 46 }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)

^{ 47 }may provide a more complete characterization of the topographic effect of treatment.

^{ 48 }Similar reduction in curvature was observed after IOP reduction in post-LASIK ectasia eyes, which are also biomechanically weaker than normal.

^{ 49 }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.

^{ 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.

^{ 16 }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,

^{ 29 }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.

^{ 54 }Characteristic nonuniformities in epithelial thickness such as those demonstrated in normal eyes

^{ 55 }and KC eyes

^{ 56 }can mask underlying stromal curvature and introduce a source of modeling error. Because the keratoconic epithelium is typically thinnest over the cone,

^{ 56 }epithelial debridement during CXL may unroof a steeper stroma, which could explain the appearance of a paradoxical steepening response in the first months

^{ 53 }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,

^{ 57 }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).

^{ 58 }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

^{ 59 }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.