Customized eye modeling for optical quality assessment in myopic femto-LASIK surgery

Refractive surgery is recognized as an effective method for myopia treatment, but it can induce night vision disturbances such as glare. We present an eye modeling method for the optical quality assessment in response to the structural changes in the eyes by femto-LASIK surgery. Customized eye models were built from the measurements of 134 right eyes pre- and post-operatively. Optical performance was evaluated using spot diagrams, point spread functions (PSFs), modulation transfer functions (MTFs), and chromatic aberrations at various fields (0°–30°), different pupil diameters (2–6 mm), and initial myopias (− 1.25 to − 10.5 D). Pupil size and initial myopia are the two major factors that affect visual performance of post-operative eyes. The results of spot diagrams, PSFs, and MTFs indicated that post-operative visual performance deteriorated as the visual field and pupil size increased, and it was significantly influenced by initial myopia. Post-operative chromatic aberrations were also affected by initial myopia. As pupil size increased, the post-operative longitudinal chromatic aberrations tended to decrease slightly, while the transverse chromatic aberrations remained similar. The use of eye modeling for refractive surgery assessment could possibly provide a more personalized surgical approach, could improve the prediction accuracy of refractive surgery outcomes, and promote the invention and development of better surgical methods.

cause of night vision disturbances is mainly considered to be that the surgery can induce high-order aberrations when the magnitudes of the low-order sphero-cylindrical aberrations are reduced. When the pupil is enlarged to adapt to a low light illumination condition, the marginal rays of vision at the transition zone of surgery tend to be diverged, so visual aberrations are produced and vision quality is deteriorated 14 . Therefore, understanding the changes in visual performance caused by refractive surgery is important to reduce the night vision disturbances and to improve patient satisfaction 15,16 .
Zernike polynomials have been used as standards to represent the low-and high-order monochromatic wavefront aberrations of the eyes 17 . It has been commonly assumed that the vision quality decreases as the root mean square (RMS) values of the high-order ocular wavefront aberrations increase, and wavefront-guided refractive surgery solutions have been developed to reduce the residual RMS values of high-order aberrations in myopic eyes [18][19][20][21] . However, the decrease or increase in the RMS values of the Zernike wavefront aberrations cannot be simply linked to the decrease or increase in the visual performance 22,23 . It has been demonstrated that one term of Zernike aberration can compensate another and improve the visual quality (e.g., the modulation transfer function, MTF) and image quality 24 . In addition, specific combinations of Zernike modes can interact to improve acuity despite an increase of the total RMS wavefront errors 22,23 . The RMS values of ocular aberrations are too rough and simple to judge optical performance, whereas the analysis of all the Zernike terms is too complicated to make the assessment. In addition, with the Zernike analysis method it is hard to describe any field-dependent aberrations or chromatic aberrations.
The eye modeling method has been used as a computational tool in visual science since the mid-nineteenth century to establish a conceptual framework for explaining optical phenomena in vision, to predict how the changes in ocular biometry affect refraction and aberrations, to explore the limitations imposed on vision by the optical system of the eye, or to provide a physical standard used to design or test the ocular imaging instruments 25 . To meet each of the aforementioned purposes, a variety of optical models of the human eye have been developed with different levels of complexity 25 -ranging from single-surface models (reduced eyes) 26 to sophisticated finite-element models that describe layered features of each component 27,28 , from models with all spherical surfaces to aspherical 29 or even high-order Zernike surfaces 30,31 , from a symmetrical structure 27 to asymmetrical structures with tilted and decentered ocular components 32 , from models with constant refractive index elements to models with gradient refractive index 33 or continuously accommodable lenses 31,34 . The use of a suitable eye model is mainly based on the principle whereby this eye model can provide valid prediction results for a particular purpose with enough simplicity, but does not necessarily rely on whether the model is anatomically or mechanistically right 25,29 . The goal of the present study was to develop an analytical eye model for accurately evaluating the optical quality change in response to the structural changes in patients' eyes by femto-LASIK surgery. To balance the functional accuracy and anatomical simplicity, we prioritized the modeling accuracy of the anterior segment topography (especially cornea) using sixth-order Zernike polynomials, while we simplified other parts of the eye (i.e., the lens, vitreous body, and retina) as symmetrical components with aspherical surfaces and constant refractive index values. Ocular geometries were clinically measured and calculated before and after the surgery from 134 right eyes, and the customized eye models were built in Zemax (Zemax, LLC). We performed the optical quality analysis for the whole eye, including the spot diagram, point spread function (PSF), MTF, and chromatic aberrations. We evaluated the imaging quality at various pupil diameters (2-6 mm), different fields of view (0°-30°), and different initial myopic diopters (− 1.25 to − 10.5 D). The use of customized eye models in the prediction of refractive surgery can relate the effect of anterior segment change (especially corneal shape change) on the retinal imaging quality performance more directly using the standard optical metrics (e.g., PSF and MTF). We aimed to provide a reliable evaluation tool with the potential to aid in the development of better surgery methods and strategies to reduce surgery-induced eye symptoms and increase post-operative visual performance.
Levofloxacin eye drops were used three days before the surgery. In femto-LASIK, an 8.5-mm corneal flap was cut at a 100-μm depth by an IntraLase FS 60 Femtosecond laser (IntraLase, Irvine, California), and the stroma was removed using an EX500 (Alcon) excimer laser that was guided by Q-algorithm 35 . The optical region was 6.5 mm, the residual corneal stroma thickness was no less than 300 µm, and the postoperative corneal thickness was no less than 400 µm. The flap was carefully repositioned after laser treatment. Patients were required to wear bandage contact lenses (Bausch) to avoid corneal flap displacement and to promote epithelial repair.  29 . The anterior segment was built from the measurement of Pentacam using a rotating Scheimpflug camera with 138,000 points in 2 seconds 36 . The measurement parameters included the radii of curvatures and asphericity values of anterior and posterior corneas, central corneal thickness, and the anterior chamber depth. Particularly, the heights of the anterior and posterior corneas were fitted to a serious of Zernike polynomials in a 6-mm diameter according to the OSA standard 17 : where z is the sag of cornea surface; R is the average corneal radius; r is a radial ray coordinate; and Q is a conic constant (asphericity). The conic constant Q is less than − 1 for hyperbolas, − 1 for parabolas, between − 1 and 0 for ellipses, 0 for spheres, and greater than 0 for oblate ellipsoids. Zernike polynomials are defined with polar coordinates (ρ, φ), where ρ is a radial coordinate ranging from 0 to 1, and φ is an azimuthal component ranging from 0 to 2π. Z i j represents a Zernike polynomial, and α i j is the coefficient for Z i j , where i represents a radial order (the degree of the polynomial), and j represents a meridional frequency (the number of cycles of sinusoidal variation across 360° of meridian). The value of i (i = 0 to n) is a positive integer or zero, and the value of j (j = 0 to m) can only be − i, − i + 2, − i + 4, … i. The negative value of j represents that the meridional component is in a sine phase, and the positive value of j indicates a cosine phase with respect to the horizontal direction. Each Zernike polynomial Z i j can be represented as 17 where δ j0 is the Kronecker delta function, δ j0 = 1 for j = 0, and δ j0 = 0 for j ≠ 0. Currently, clinical measurements of lens biometric parameters-such as radii of curvature, asphericities, and refractive index distribution-are challenging to obtain 37 . Previous studies have shown very weak relationships between the lens biometric parameters and refraction, such as P-value = 0.867 between the lens thickness www.nature.com/scientificreports/ and refraction, and P-value = 0.848 between lens power and refraction 25 . In this study, we assumed that the lens biometric parameters were unaffected by refraction, and we adapted previous estimation methods 25,[38][39][40][41][42] to build the parameters of the crystalline lens in each eye model. We set the tilt and decentration of the lens as zero to simplify the estimation of the lens parameters. The equivalent lens power P L can be calculated based on the measured parameters using Bennett's method 39 : where S is the spherical refraction of the eye (unit: D); P C is the power of the cornea (unit: D); T AC , T L , and T V are the thickness/depth of the anterior chamber, lens, and vitreous body, respectively (unit: mm); n h = 1.337 is the refractive index of the aqueous humor; c 1 T L is the distance between the anterior lens surface and the first lenticular principal plane, and c 2 T L is the distance between the posterior lens surface and the second lenticular principal plane, where the constants are c 1 = 0.571 and c 2 = − 0.378 39,40 . We used the T AC values from Pentacam measurements and the T L and T V values from Lenstar measurements. We estimated the radii of curvature of the arterial lens surface (R La ) and the posterior lens surface (R Lp ) using a linear regression fitting result provided by Rozema et al. 41 : The values of asphericities (conics, or Q-values) of the unaccommodating lens vary over a wide range in different eye models and literature due to the difficulties of measurement 25,27,29,34,37,[42][43][44] . Here we estimated the asphericities of the arterial and posterior lens surface using the fixed values of − 3.13 and − 1, respectively, based on the Navarro eye model 34 .
The retinal shapes of emmetropic and myopic eyes were described by Atchison et al. using a non-rotationally symmetrical ellipsoid fitting method 45 . They found the shapes of retinas were oblate in most of emmetropic eyes. As myopic diopters increased, the ellipsoid dimensions increased, with the axial dimension increasing more than the vertical dimension, which in turn increased more than the horizontal dimension 45 . The mean values of the vertex radii of curvature (R R ) and asphericities (Q R ) of retina are described as 45 In all cases, an object was placed at infinity so that the entered lights had plane wavefronts without any aberrations, and the crystalline lenses were unaccommodated. The object was ideally focused on the retinas at the visual field of 0° for all eyes. For the pre-operative eye models, ideal spectacles were used to achieve the optimized visual acuity since most myopic subjects wear spectacle lenses in their daily lives to compensate for the refraction errors. For the post-operative eye models, the retina positions were shifted to focus the light beams. Zemax simulations were performed on the pre-and post-surgical eyes to compare their optical performances, and these included the metrics of the geometry-optics-based spot radius diagram and chromatic aberration, and the fast Fourier transform (FFT)-based polychromatic PSF and polychromatic MTF.
Statistics. The pre-and post-operative ocular diopters and corneal geometry parameters were compared.
Linear regressions were performed to illustrate the surgically induced changes in these parameters in response to the initial myopia. The Pearson correlation coefficient (r) and the level of significance (P value) were used to investigate the strengths of the correlations. The value of r ranges from − 1 to 1, where r = 1 indicates a perfect positive correlation, r = − 1 indicates a perfect negative correlation, and r = 0 indicates no linear correlation. The statistics were considered significant when P < 0.05, and extremely significant when P < 0.01. Optical performance was evaluated on eye models based on the initial myopic diopters, which were mild myopia group (0 to − 3 D), medium myopia group (− 3 to − 6 D), and high myopia group (over − 6 D). The optical quality metrics of PSF, MTF, and chromatic aberrations were analyzed at different pupil diameters (2-6 mm) and at different fields of view (0°-30°), and they were compared based on the mean and standard deviation (SD) values for each initial myopia group.

Results
Ocular parameters. Each eye model was composed of two Zernike surfaces representing the anterior and posterior corneas, a stop at the pupil plane to represent the iris, two aspherical surfaces representing the crystalline lens, and an aspherical retina, as shown in Fig. 1. The ocular parameters (mean ± SD) from 134 right eyes are shown in Table 1. The change in corneal geometry is the major outcome of the femto-LASIK surgery and has been considered as the major source of the production of high-order aberrations and night vision disturbances [14][15][16] ; while the parameters of the lens, vitreous body, and retina for each patient were assumed to be unchanged after surgery. Corneal parameters before and after the femto-LASIK surgery were compared in response to the initial myopia (Fig. 2). The pre-surgical corneal parameters were not correlated with the initial myopia. The reductions of corneal central thickness (0.54 ± 0.03 mm to 0.45 ± 0.04 mm) and corneal diopter (42.44 ± 1.38 D to 37.84 ± 2.08 D) were directly caused by the surgery [P-values < 0.01, Figs 2a,b]. The anterior corneal geometries were greatly altered by the surgery, and were correlated with the intimal myopic diopters with P-values smaller than 0.01 [ Fig. 2c-f]. Due to the femto-LASIK surgery, the anterior corneas were flatter, while the cornel radii increased from 7.77 ± 0.25 mm to 8.56 ± 0.43. The asphericity (Q) values of anterior cor-  Optical performance analysis. Optical performance was evaluated and compared between the preoperative eyes with spectacles and the postoperative eyes without spectacles at various fields (0°-30°), pupil diameters (2-6 mm), and initial myopias (− 1.25 to − 10.5 D). Ocular chromatic aberration is an optical distortion in which the focusing positions, in either the longitudinal (axial) or the transverse (lateral) direction, are wavelengthdependent [46][47][48] . Figure 3 shows . Before the femto-LASIK surgery, neither the longitudinal nor the transverse chromatic aberrations were dependent on pupil size or initial myopia. Both aberrations increased slightly after the surgery, and increased as the initial myopic diopter increased. As the pupil size increased, the longitudinal chromatic aberrations after refractive surgery tended to decrease slightly, while the transverse chromatic aberrations remained similar. The degree of spreading (blurring) of a point source at the image plane can be described by spot diagrams using the ray-tracing method and by polychromatic PSF using the FFT method. The Strehl ratio is defined as the ratio of the light at the peak of the diffraction pattern of an aberrated image to that at the peak of an ideal image. Figure 4a-c show the spot size evaluation results using the RMS radius values, and Fig. 4d-f evaluate the Strehl ratio for the polychromatic PSF. In general, the spot sizes were enlarged and the Strehl ratios were reduced as the field of view increased from 0° to 30°, indicating that the optical performance decreased from the central retina to the peripheral retina. When the pupil size was small (i.e., 2 mm), the optical performance was similar for all eyes before and after surgery. As the pupils enlarged from 4 to 6 mm, the pre-operative eyes had similar optical performance for all the myopic groups, but the post-operative eyes were significantly influenced by the initial myopia. After surgery, the optical performance of eyes that were initially more myopic tended to be much worse at larger pupil sizes, as demonstrated in Fig. 4c,d. Figure 5 presents the MTFs (mean values between the tangential and sagittal directions) for the pre-and postsurgical eyes at a maximum spatial frequency of 60 cycles per millimeter. In general, MTF declined as spatial frequency increased, and from central retina to peripheral retina. The post-surgical MTF was significantly affected Table 1. General geometry of eye models for 134 right eyes before and 6 months after femto-LASIK. *The magnitudes of anterior and posterior corneal surfaces were fitted as sixth-order Zernike polynomials. [ www.nature.com/scientificreports/ by the initial myopia, and declined as the pupil size increased. MTFs of more myopic eyes tended to worsen at larger pupil sizes, as demonstrated in Fig. 5i. This result agreed with the spot size and PSF analysis (Fig. 4).

Discussion
We have built customized eye models for 134 eyes before and after the femto-LASIK surgery, and we have provided quantitative and objective evaluations of changes in retinal imaging quality associated with the structural changes in the human eye induced by the surgery. The results of RMS spot size, Strehl ratio, and MTF analysis have shown that the optical performance is reduced from the central retina to the peripheral retina (e.g., from 0° to 30°). When the pupil size is small (e.g., 2 mm in diameter), the optical performance is similar for all eyes before and after surgery. As the pupil enlarges (e.g., up to 6 mm in diameter), the pre-operative eyes have similar optical performance for all the myopic groups, but the post-operative eyes are significantly influenced by the initial myopia. After femto-LASIK surgery, the optical performance with higher initial myopia tends to deteriorate at a larger pupil size. These results are consistent with previous studies on the influence of pupil size 14,49 and initial myopia 50,51 on the post-surgery optical performance. The post-operative chromatic aberrations also increase as the initial myopia increases after the surgery, but have different trends with the pupil size compared with other aberrations, such as monochromatic coma and spherical aberrations. The longitudinal chromatic aberrations after refractive surgery tend to decrease slightly as the pupil size increases, while the transverse chromatic aberrations remain similar for different pupil sizes. To the best of our knowledge, the finding of the post-surgical longitudinal chromatic aberrations decreasing with pupil size has not been presented previously.
In general, when we consider the contributions to overall aberrations, pupil size and initial myopia are the two major factors that induce higher-order aberrations and decrease visual performance. The accuracy of the eye model is dependent upon the precision of the measurements, especially the measurement of the crystalline lens. The human crystalline lens is a complex and inhomogeneous optical component with gradient refractive index distribution and aspheric surfaces 33 . We know that the lens shape and gradient index distribution may help to reduce the spherical aberration and coma that originate in the cornea 32 , but how the lens curvature and index gradient contribute to lens aberrations is still poorly understood due to a lack of reliable measurements on the lens shape outside of the central zone 37 . To date, the in vivo measurement of the posterior lens remains a longstanding challenge in the clinic. Optical methods, such as Scheimpflug photography 42 and anterior-segment optical coherence tomography (OCT) 52,53 , are incapable of penetrating the iris tissue and imaging the peripheral region of the lens beneath the iris. Ultrasound 54 and magnetic resonance imaging (MRI) 44 can   44 . The estimation differences of Q-values usually arose from the limited methods of measurement and the fitted diameter zone 37 . Simulation of the unaccommodated lens is already a great challenge, and enabling the eye model to account for the accommodation of the crystalline lens is far more difficult 34,44,56 . The measurement of retina shape is also limited by the optical distortions of the optical elements before the retina, particularly the gradient index distribution of the crystalline lens 57 . Therefore, better ocular measurement methods are required to improve the modeling accuracy of the human eye and to improve the simulation precision for the optical ray tracing method. It should be noted that there is no perfect optical model of the eye that is best for every purpose 29 . Each of the previous eye models could be either a simple single-surface model to represent the on-axis chromatic and spherical aberration 26 or a sophisticated model for anatomical accuracy 27 . An appropriate eye model is the one that gives valid results for a particular purpose, so a more complicated model is not necessarily better 29 . The main purpose of this study was to describe the effect of the post-surgical corneal structure change on visual performance. To meet this purpose, we prioritized the modeling of the anterior segment (especially cornea) while simplifying other parts of the eye, such as the lens, vitreous body, and retina. As is well known, the human eye is a decentered optical system with a non-rotationally symmetric structure 17 , where each element (e.g., cornea, pupil, and lens) can be decentered and tilted, and the photoreceptors (i.e., cones and rods) are not uniformly distributed at the retina 25,32,58 . We started the simulation of the human eye from a rotationally symmetric schematic eye model in an unaccommodated condition 29 . The asymmetricity of the anterior and posterior corneal surfaces can still be sufficiently represented by the surface heights in Zernike terms (e.g., tilt: Z 1 ±1 , astigmatism: Z 2 ±2 , and coma: Z 3 ±1 ), but neither the tilt nor the decenter values of other ocular elements were presented in this model. We considered this as a valid simplification process. If we had used the empirically estimated (rather than www.nature.com/scientificreports/ trustworthily measured) decenter or tilt values for all the ocular elements of the model, it would have led to more system complexity and slower computing time, but it would not have improved the optical analysis accuracy. For the same reason, we assumed the gradient-index crystalline lens 25,27 as homogeneous, and we assumed that all of the elements had homogeneous refractive indexes and Abbe-numbers 29 . Another issue is that we only evaluated the optical performance in each eye model, and ignored the neural factors that can also affect the visual performance. The combination of an individual's optical and neural transfer functions could likely be better to predict the actual visual performance 23 . A neural transfer function expresses the loss of modulation as spatial frequency increases in the process of converting the optical signals to neural impulses and transferring the information through to the neural system and out as a percept 23 . Visual perceptual sensitivity is field-dependent due to the uneven distribution of photoreceptors and optic nerve cells from the fovea to the peripheral retina. This factor can be modeled optically using an apodization filter at the pupil plane 59 . Post-receptoral neural processing of the unevenly sampled retinal image affects the processing of blurred retinal images in a manner that can also be constructed as a form of apodization, or as a mathematical convolution of the optical PSF with a neural PSF 60 . The main focus of this study was the optical performance, but we will consider the effect of the interactive optical-neural performance on visual perception in our following studies. Although designed for Femto-LASIK, this eye modeling method could be utilized in other types of myopic surgeries 4 and corneal surgeries, such as treatments of hyperopia 61 and keratoconus 62 .
Corneal biomechanics, such as elasticity and viscosity, is another factor that may also contribute to the induction of high-order aberrations and the deterioration of optical performance. Integrating biomechanical factors into a more sophisticated eye model could potentially improve the accuracy of the refractive surgery prediction 63 . Corneal biomechanics maintain the ocular rigidity and are inherently tied to corneal health and visual performance 64 . Any morphological changes of the cornea are usually accompanied by biomechanical changes, and vice versa. Myopic refractive surgery may weaken the corneal strength and cause long-term instability 8,65 and may rarely lead to several complications, such as corneal ectasia 66,67 ; corneal biomechanics can, in turn, influence the predictability, stability, and safety of refractive surgery 68,69 . The pre-operative evaluation of corneal biomechanics is key to preventing ectasia as well as predicting and evaluating treatment outcomes 4 , but the clinical methods to precisely access corneal biomechanics are still limited. Developing more advanced imaging and quantitative methods, such as high-sensitivity optical coherence tomography (OCT)-based elastography methods, to more reliably access corneal biomechanics in clinic is an active research area in vision science and is of great interest to our research group [70][71][72] . Once the new corneal elastography tool and the measurable data of corneal biomechanics (such as natural frequency 73,74 and Young's modulus [75][76][77] ) are provided, we could possibly provide an updated eye model that involves comprehensive information of ocular structure, biomechanics, optical performance, and neural functions. The application of such an eye model for refractive surgery could provide better pre-operative prediction of surgery outcomes, promote better surgical methods with improved visual performance, prevent post-surgical complications (e.g., corneal ectasia) more successfully, and reduce post-surgical patient complaints more effectively 78 .