**Purpose**:
To assess the predictive accuracy of simulation-based LASIK outcomes.

**Methods**:
Preoperative and 3-month post-LASIK tomographic data from 20 eyes of 12 patients who underwent wavefront-optimized LASIK for myopia were obtained retrospectively. Patient-specific finite element models were created and case-specific treatment settings were simulated. Simulated keratometry (SimK) values and the mean tangential curvature of the central 3 mm (*K*_{mean}) were obtained from the anterior surfaces of the clinical tomographies, and computational models were compared. Correlations between *K*_{mean} prediction error and patient age, preoperative corneal hysteresis (CH), and corneal resistance factor (CRF) were assessed.

**Results**:
The mean difference for *K*_{mean} between simulated and actual post-LASIK cases was not statistically significant (−0.13 ± 0.36 diopters [D], *P* = 0.1). The mean difference between the surgically induced clinical change in *K*_{mean} and the model-predicted change was −0.11 ± 0.34 D (*P* = 0.2). *K*_{mean} prediction error was correlated to CH, CRF, and patient age (*r* = 0.63, 0.53, and 0.5, respectively, *P* < 0.02), and incorporation of CH values into predictions as a linear offset increased their accuracy. Simulated changes in *K*_{mean} accounted for 97% of the variance in actual spherical equivalent refractive change.

**Conclusions**:
Clinically feasible computational simulations predicted corneal curvature and manifest refraction outcomes with a level of accuracy in myopic LASIK cases that approached the limits of measurement error. Readily available preoperative biomechanical measures enhanced simulation accuracy. Patient-specific simulation may be a useful tool for clinical guidance in de novo LASIK cases.

^{1}However, LASIK and other forms of corneal refractive surgery that rely on tissue incision or removal induce corneal biomechanical changes.

^{2–4}These changes to the corneal structure and individual differences in the mechanical response to them are important potential drivers of residual postoperative refractive error, postoperative refractive regression, and postoperative corneal ectasia.

^{1}A variety of risk factors for retreatment after primary refractive surgery have been identified, including the type and magnitude of preoperative refractive error, depth of ablation, patient age, and higher degrees of preoperative astigmatism.

^{5–7}In addition to the cost and intrinsic risks of additional surgery, retreatments have been associated with an increased risk of corneal ectasia.

^{8,9}Although historic nomograms have been useful for improving refractive outcomes through regression-based guidance of treatment offsets,

^{10}they are generally agnostic to individual corneal shape variations and other important drivers of the structural response to surgery. A surgical planning tool that can increase the accuracy of surgical algorithms for routine and atypical cases has the potential to reduce cost while increasing patient safety and satisfaction.

^{11}Structural analysis using the FEA method has the potential to account for the biomechanical impact of a complex three-dimensional (3D) surgical intervention while accounting for the patient's particular corneal geometric configuration and biomechanical properties. The sensitive link between the cornea's geometry, biomechanical behavior, and the optical function of the eye presents significant advantages for leveraging simulation-based engineering in keratorefractive applications. Furthermore, the cornea's accessibility with commercially available imaging systems makes creation of patient-specific models more feasible than in other tissues requiring invasive or ionizing imaging techniques.

^{12,13}to more microstructurally representative models

^{14–21}that account for preferred collagen orientations from x-ray diffraction data,

^{22}axially inclined collagen fibers in the anterior stroma characterized by two-photon microscopy

^{23}and corneal hydration behavior.

^{18}Only a limited number of studies have modeled the optical effects of both flap creation and ablation in the context of LASIK,

^{24,25}and no prior publication has performed a direct clinical comparison of computationally predicted refractive outcomes. The goal of the present study was to assess the predictive accuracy of case-specific finite element simulations of LASIK procedures through comparison to clinical topographic and refractive outcomes.

^{14,26,27}To simulate this spectrum of behavior, a microstructural fiber-reinforced material model representing collagen fibers was combined with an isotropic Neo-Hookean solid extrafibrillar matrix to produce a composite corneal stroma.

^{28–30}Collagen crimping behavior was accounted for by representing each collagen fibril as a 3D spring proposed by Freed et al.

^{28,29}where

*H*is the axial distance between two turns of the helix,

*R*is the radius of the helix, and

*r*is the radius of the fiber thickness. The microstructural helical spring model requires three parameters:

*Ho/Ro*(the ratio between initial axial distance between two turns within the helix and the initial radius of the helix),

*Ro/ro*(the initial radius of the helix and the radius of the fiber thickness), and

*Ef*(the elastic modulus of the fiber when it is uncrimped and behaves in its linear region).

^{15,28}The resulting fiber behavior was calculated using the algorithm given elsewhere by Freed et al.

^{29}

^{22,31–33}The material formulation proposed by Freed et al.

^{29}was limited to reflect this behavior in the fiber splay. Therefore, the algorithm was combined with the equation given by Pinsky et al.

^{15}to more realistically simulate the corneal stroma (Equations 1–4). At each integration point within the model, splay of fibers was modeled and angularly integrated to represent a gradual change in fiber orientation from the center to the periphery of cornea. The equations are partially reproduced here to demonstrate how fiber and matrix components are accounted for in the total strain energy density function (Equation 1) and how the fiber component is defined using the probability density function (Φ) to describe the presence of fiber population as a function of rotational angle (

*θ*) and radial distance from the corneal center (

*R*) (Equations 2–4). The separate treatment of fiber orientation for central (Equation 3) and peripheral (Equation 4) cornea allows representation of the transition from predominantly orthogonal to more circumferential/tangential orientation near the limbus. A linear combination of Equations 3 and 4 was used in the transition zone between central and peripheral regions. The details of the components of this equation including

*n*,

*c*

_{1}, and

*c*

_{2}were explained by Pinsky et al.

^{15}

*C*

_{10}, in the material model down to 50% from anterior to posterior regions of the stroma to approximate experimental data in human corneal explants.

^{34}The material model was implemented in an Abaqus v611 UMAT subroutine (Dassault Systemes Simulia Corp., Providence, RI, USA).

^{30}The pretreatment corneal material constants were obtained by performing an inverse FEA

^{30}on previously published experimental inflation and tensile test data from a 50-year-old excised human cornea.

^{35}This analysis produced material constants that were then used in the models to represent the normative preoperative stroma. Preoperative material constants were defined as follows: the shear modulus of isotropic ground matrix

*C*= 0.04 MPa, the compressibility constant

_{10}*D*= 0.0001 MPa

^{−1}, the ratio of the axial length of the collagen fibers spring to the diameter of the collagen fiber

*H*= 30.5, the unitless ratio of radial diameter of the collagen fiber to the diameter of the collagen fiber

_{0}r_{0}*R*= 1.51, and the elastic modulus of the collagen fibers

_{0}r_{0}*Ef*= 32 MPa. Epithelium was modeled as a uniformly thick 50-μm isotropic hyperelastic neo-Hookean material. The shear modulus of the epithelium was set as 1% of the shear modulus of the stroma, and the fiber component was omitted due to absence of continuous collagen elements and the epithelium's negligible contribution to the corneal structural response.

^{20}The preoperative geometry was divided into four regions in depth (Fig. 1B): epithelium (50 μm), stromal component of the flap (50 μm), flap interface wound (10 μm), and posterior stroma (defined as the remaining patient-specific thickness).

^{36}A generic sclera, with an anterior thickness matched to the patient-specific peripheral corneal thickness, was generated to form the posterior portion of the eye in order to capture the effects of limbal displacement within the simulation.

^{20,24}The posterior surface of the corneoscleral globe was then subjected to a homogenous IOP of 15 mm Hg.

^{37}The patients received wavefront optimized treatments. The ablation profile calculation and equations were given by Mrochen et al.

^{37}Finite element analysis were run using a commercial FEA solver Abaqus v6.11 (Dassault Systemes Simulia Corp.). The amount of ablation was calculated for each node in the 3D finite element model and the nodes were repositioned according to their ablation magnitudes in the geometry with a custom Python script. Because the corneal Scheimpflug scans were acquired under physiological IOP, the stress-free geometries were calculated by an iterative technique given elsewhere by Pandolfi et al.

^{14,38}and modified by our group.

^{25,39}The corneal surface coordinates of the unoperated mesh (

*X*

_{0}) depend on the patient specific tomographies acquired under pressure. The nodal coordinates of the mesh within the treatment zone were modified from

*X*

_{pre}to

*X*

_{post}according to the ablation algorithm. The coordinates of the postoperative mesh (including the treatment zone and the peripheral cornea, limbus, and sclera) were taken as a reference (

*X*

_{ref}). A negative pressure with amplitude 90% lower than the physiologic IOP (15 mm Hg) was applied to the posterior surface of the postoperative corneoscleral model (

*X*

_{ref}). The displacements that were calculated from the negative pressure (

*u*

_{deflation}) were added to

*X*

_{ref}to obtain the initial set of coordinates (

*X*

_{0}) to start the iterative no-load calculation (Equation 5). Following the deflation, the physiologic IOP (15 mm Hg) was iteratively applied and error (

*e*) was calculated as the root mean square (RMS) of the distances between the nodes in the whole corneoscleral structure from the reference and the calculated geometries (Euclidian norm). The iterations were stopped once RMS error was less than 1 × 10

^{−5}, which typically required six or seven iterations.

^{40}The side cut was simulated with a 50% weakening in the isotropic component and omitting the fiber component of the strain energy equation (Equation 6). The flap was simulated by omitting the fiber component of the strain energy equation (Equation 7).

^{36}All patients were operated on using a circular flap with 9-mm diameter and 90° side cut angle. Hinge width was reported as 55° for all patients at 12 o'clock position. All simulations conformed to these reported surgical parameters for each patient.

^{41}The material properties of the posterior cornea remained unchanged for the LASIK simulation.

*K*values including mean tangential curvature from the central 3-mm region (

*K*

_{mean}), the average curvature of the steepest meridian (

*K*

_{1}), and the average curvature of the flattest meridian (

*K*

_{2}) of the anterior corneal surface were calculated. Tangential curvature at each point (

*K*

_{tan}) on the cornea was calculated using the first and second derivative of an 8

^{th}order Zernike fit (Z) of the anterior corneal surface with respect to distance from the center of the cornea (

*r*

_{cornea}) (Equation 9). The constants of Zernike polynomials (Z) were calculated with a linear regression analysis using a custom Python script. The calculated tangential curvature values were converted to diopters (D) using Snell‘s law with a corneal refractive index of 1.3375 (the keratometric index, which estimates the total corneal power assessed at the anterior corneal surface). Pretreatment clinical K

_{1}and K

_{2}values were directly taken from Pentacam tomographies for comparison. Pretreatment model, posttreatment clinical, and simulated K

_{1}and K

_{2}values were calculated using SpecifEye. The pairwise differences of these variables for clinical and simulated geometries were calculated. Statistical significance was determined by paired Student's

*t*-tests for nonindependent samples with

*P*less than 0.05 indicating significance (MINITAB, Minitab, Inc., State College, PA, USA).

*K*

_{mean}from postoperative

*K*

_{mean}for actual and simulated cases, respectively. Treatment and prediction errors were calculated as differences between actual – programmed changes and predicted – actual changes, respectively.

*K*

_{1},

*K*

_{2,}corneal astigmatism (

*K*

_{1}–

*K*

_{2}), steep axis, and

*K*

_{mean}values and comparisons of actual and simulation values for the preoperative and post-LASIK states, respectively. Preoperative model fidelity was verified by comparing each value from clinical tomography to that of the associated finite element model, and none of the differences were statistically significant. The mean difference in

*K*

_{mean}(Δ

*K*

_{mean}) between actual outcomes of the surgeries and the predictions produced by their corresponding LASIK simulations was not statistically significant (−0.13 ± 0.36 D,

*P*= 0.1 by paired Student's

*t*-test). The only statistically significant difference between clinical and simulated values across the preoperative and postoperative comparisons was for

*K*

_{2}(+0.13 ± 0.28 D,

*P*= 0.03). Errors in the agreement of the axis of steep astigmatism were also insignificant and were higher in individual cases with very low magnitudes of astigmatism where axis identification (both by refraction and by keratometric analysis) is anticipated to be error prone but of minimal clinical consequence.

*K*

_{mean}) between right and left eyes of the eight patients who had both eyes enrolled. The analysis suggested no significant intraclass correlation (

*r*= 0.11,

*P*= 0.95). The minimum detectable difference for

*K*

_{mean}was 0.23 D for a paired

*t*-test using values

*n*= 20,

*α*= 0.05, and

*β*= 0.2 (80% statistical power) by post-hoc calculation. All predicted

*K*

_{mean}values were within 1.00 D of the clinical postoperative value, with 55%, 85%, and 95% of the cases demonstrating prediction errors within 0.25, 0.50, and 0.75 D, respectively (Fig. 2). Figure 2 also presents a prediction error frequency based on an adjustment for preoperative biomechanical measurement and is further described below.

*K*

_{mean}expressed as a function of CRF, CH, and age. Forty percent of the variance in

*K*

_{mean}prediction error was accounted for by interindividual differences in preoperative CH. Twenty-nine percent was accounted for by CRF, and 25% by age (

*P*= 0.016, 0.004, and 0.02, respectively). Δ

*K*

_{mean}values were confirmed to be normally distributed, and a one-sided outlier analysis using Grubbs' test for the null hypothesis that all data points come from the same normally distributed population did not suggest that the largest prediction error, −1.02 D, was a statistical outlier (

*P*= 0.08). When the linear regression equation for the strongest preoperative predictor of prediction error (CH) was used as a linear adjustment to the model-predicted curvature, the magnitude and frequency of observed error improved (Fig. 2).

*K*

_{mean}) and assesses treatment error and prediction errors. The treatment error, expressed as the difference between programmed refractive change calculated at the anterior corneal plane and actual change in

*K*

_{mean}, was statistically significant (+0.40 ± 0.30 D,

*P*< 0.05). However, the model prediction error in

*K*

_{mean}was not statistically significant (−0.11 ± 0.34 D,

*P*= 0.2).

*P*= 0.99) and the scatterplot graph (Fig. 6) also demonstrates a strong correlation between predicted changes in corneal curvature and corneal vertex–corrected changes in actual SE manifest refraction.

^{1}Surgeons frequently use personal nomograms or commercial software-based equivalents as a guide for modifying the laser input to reduce the risk of outcome variability. Although planning systems vary in complexity, the process involves entry of a small number of patient- and surgery-specific variables, interrogation of a regression formula developed from historic outcome data as a function of some subset of those inputs, and generation of a predicted refractive outcome and/or suggested treatment entry for generating the desired correction. These tools have been an important part of outcomes analysis and optimization since the inception of corneal refractive surgery, and analogous methods are also widely used in incisional keratotomy, intraocular lens calculation, intracorneal ring implantation, and other forms of ocular surgery.

^{42}The variables represented in nomograms are sparse representations of the patient and the procedure and do not leverage potentially predictive 3D corneal tomography datasets. The statistical associations produced by regression models do not explicitly reflect the biophysical mechanisms of variability and are therefore of limited utility for de novo inferences in complex or atypical cases. Also, aside from surrogate associations related to patient age, nomograms do not explicitly account for the biomechanical impact of the intervention on corneal geometry and optical outcome.

^{2,4,26}

^{24,36,43}the current study is the first to assess the predictive performance of a model in a series of keratorefractive surgery patients. This is a critical step in assessing computational modeling as a potential clinical decision support tool in the surgical planning process.

*K*

_{mean}) derived from each computational model from the same metric measured in the actual post-LASIK tomographies—an absolute comparison of predictive accuracy—was insignificant (−0.13 ± 0.36 D,

*P*= 0.1). In relative terms, the mean difference between the surgically induced clinical change in

*K*

_{mean}and the model-predicted change was also insignificant (−0.11 ± 0.34 D,

*P*= 0.2). The keratometric variable

*K*

_{mean}was selected as the primary outcome measure because it was thought to represent the most optically important changes within a typical photopic entrance pupil and because the data directly represent the input and output of the modeling process. However, to relate the predictive performance of the model more directly to clinical refractive outcomes, a direct comparison of predicted changes in model

*K*

_{mean}to the vertex-adjusted SE manifest refractive change was performed and produced close agreement across the range of corrections (Table 4; Fig. 6) with a mean error of +0.23 ± 26 D. While

*K*

_{mean}considers only corneal first-surface changes in curvature, simulation results compared favorably with subjective manifest refraction results that depend on whole-eye optics. Ray tracing analysis is being incorporated into the modeling software to account for additional sources of refractive change such as axial length and posterior corneal surface changes.

*K*

_{mean}change (0.40 ± 0.30 D) was compared with the prediction error in the same metric (−0.11 ± 0.34 D), the ranges of clinical and prediction error were comparable. However, the model prediction was closer in magnitude to the actual

*K*

_{mean}change than the programmed correction by a factor of 4. Individual comparisons of these differences in Table 3 suggest that preoperative adjustments to the surgical plan based on case-specific model predictions would potentially have improved treatment accuracy in 16 of 20 cases. The actual programmed correction for each case was used to drive each associated simulation and included any spherical and/or cylindrical treatment offsets incorporated by the surgeon after reviewing commercial nomogram software output.

*K*

_{mean}(an apparent 1.18-D myopic undercorrection) and also the largest prediction error in

*K*

_{mean}(−1.08 D) also had the lowest preoperative CH value. While CH is not a measure of the corneal elastic modulus, the trend toward undercorrection of myopia in eyes with lower CH and overcorrection in those with higher CH reflects earlier FEA-based sensitivity analyses suggesting greater anterior displacement of the residual stromal bed and thus myopic undercorrection in eyes with lower corneal stiffness properties.

^{44}For simulation purposes, every patient was assumed to have identical, normative initial material properties at the anterior stroma. The isotropic component of these properties were then reduced as a function of depth as described in the methods, and thus produced different depth-dependent property profiles due to differing patient-specific corneal thicknesses. However, the absence of knowledge regarding actual patient-specific material properties is a potential source of prediction error in these simulations, and the finding that a linear adjustment for CH could improve prediction fidelity is important evidence that incorporation of individual measures of corneal properties in the future could enhance model performance. On average, the models predicted slightly flatter corneas than were observed in actual cases, and this could be due to a slightly stiff corneal material assumption.

^{45,46}coupled with the ability of emerging OCT and high-frequency ultrasound devices to generate patient-specific epithelial thickness maps

^{47}will provide even greater opportunities for improved model accuracy.

*. 2009; 116: 691–701.*

*Ophthalmology**. 2005; 31: 175–184.*

*J Cataract Refract Surg**. 2007; 23: 800–807.*

*J Refract Surg**. 2001; 17: 658–669.*

*J Refract Surg**. 2003; 110: 748–754.*

*Ophthalmology**. 2009; 25: 273–276.*

*J Refract Surg**. 2016; 35: 607–612.*

*Cornea**. 2002; 109: 1991–1995.*

*Ophthalmology**. 2003; 110: 267–275.*

*Ophthalmology**. 2010; 117: 1236–1244. e1.*

*Ophthalmology**Simulation-Based Engineering Science: Revolutionizing Engineering Science through Simulation*. Arlington, VA, May 2006. Available at: https://www.nsf.gov/pubs/reports/sbes_final_report.pdf.

*. 1996; 118: 473–481.*

*J Biomech Eng**. 2006; 34: 1628–1640.*

*Ann Biomed Eng**. 2008; 130: 061006.*

*J Biomech Eng**. 2005; 31: 136–145.*

*J Cataract Refract Surg**. 2013; 12: 1101–1113.*

*Biomech Model Mechanobiol**. 2012; 53: 873–880.*

*Invest Ophthalmol Vis Sci**. 2015; 12: 20150241.*

*J R Soc Interface**. 2013: 7: 0409131–0409132.*

*J Med Dev**. 2014; 40: 943–953.*

*J Cataract Refract Surg**. 2015; 42: 76–87.*

*J Mech Behav Biomed Mater**. 2001; 20: 95–137.*

*Prog Retin Eye Res**. 2013; 54: 7293–7301.*

*Invest Ophthalmol Vis Sci**. 2009; 25: 875–887.*

*J Refract Surg**. 2011; 52: 9174–9187.*

*Invest Ophthalmol Vis Sci**. 2006; 83: 709–720.*

*Exp Eye Res**. 2005; 31: 136–145.*

*J Cataract Refract Surg**. 2005; 127: 587–593.*

*J Biomech Eng**. 2005; 4: 100–117.*

*Biomech Model Mechanobiol**. In press.*

*J Refract Surg**. 2009; 28: 369–392.*

*Prog Retin Eye Res**. 2011; 52: 1243–1251.*

*Invest Ophthalmol Vis Sci**. 2006; 47: 901–908.*

*Invest Ophthalmol Vis Sci**. 2008; 24: S85–S89.*

*J Refract Surg**. 2007; 32: 11–19.*

*Curr Eye Res**. 2014; 40: 971–980.*

*J Cataract Refract Surg**. 2004; 30: 775–785.*

*J Cataract Refract Surg**. 2013; 35: 211–216.*

*Med Eng Phys**. 2013; 7: 409131–409132.*

*J Med Device**. 2005; 123: 741–756.*

*Arch Ophthalmol**. 2005; 21: 433–445.*

*J Refract Surg**. New York, NY: Chapman & Hall; 1993.*

*Predictive Inference: An Introduction**. 2014; 40: 905–917.*

*J Cataract Refract Surg**. 2009; 25: 875–887.*

*J Refract Surg**. 2015; 122: 677–686.*

*Ophthalmology**. 2016; 123: e5–e6.*

*Ophthalmology**. 2012; 28: 195–201.*

*J Refract Surg*