An innovative approach for determining the customized refractive index of ectatic corneas in cataractous patients

The aim of this study is to determine the customized refractive index of ectatic corneas and also propose a method for determining the corneal and IOL power in these eyes. Seven eyes with moderate and severe corneal ectatic disorders, which had been under cataract surgery, were included. At least three months after cataract surgery, axial length, cornea, IOL thickness and the distance between IOL from cornea, and aberrometry were measured. All the measured points of the posterior and anterior parts of the cornea converted to points cloud and surface by using the MATLAB and Solidworks software. The implanted IOLs were designed by Zemax software. The ray tracing analysis was performed on the customized eye models, and the corneal refractive index was determined by minimizing the difference between the measured aberrations from the device and resulted aberrations from the simulation. Then, by the use of preoperative corneal images, corneal power was calculated by considering the anterior and posterior parts of the cornea and refractive index of 1.376 and the customized corneal refractive index in different regions and finally it was entered into the IOL power calculation formulas. The corneal power in the 4 mm region and the Barrett formula resulted the prediction error of six eyes within ± 1 diopter. It seems that using the total corneal power along with the Barrett formula can prevent postoperative hyperopic shift, especially in eyes with advanced ectatic disorders.


Results
In this study, 7 eyes of 5 patients with PMD or KCN were evaluated. Six eyes with KCN were categorized into moderate and severe groups based on Amsler-Krumeich classification (Based on this classification, the average keratometric readings of ≤ 48 D is considered as mild, > 48 D and ≤ 55 diopters is considered as moderate, and > 55 is considered as severe 13 ). These eyes had been under phacoemulsification surgery and intraocular lens implantation. Demographic information, ectasia type, keratometry (measured by Pentacam), and biometry values before each patient's surgery are shown in Table 1. In addition, Table 2 shows the type and power of the implanted IOL and the subjective refraction of each patient measured at least 3 months after surgery.
Determining the eyes customized refractive index. After designing the intraocular lens and creating a customized model for each patient's eye, the ray tracing method was implemented by using optical analysis software. The aberrations of each eye were calculated using this method. Then by applying the optimization process, we tried to minimize the difference between the aberrations of simulation and those measured from the  Table 2. Type and power of the implanted IOLs and refraction of the patients after surgery. www.nature.com/scientificreports/ device. Table 3 demonstrates the absolute value of the difference between simulated aberrations and measured ones in different eyes for different Zernike coefficients. By minimizing the difference between aberrations, the unknown parameter (i.e., refractive index of the cornea), was calculated by the optimization method. These two parameters were specifically obtained for the studied eyes, which are shown in Table 4. As you can see, the refractive indexes of the ectatic corneas deviated from the standard value (1.376), so that the refractive index increased in four of seven eyes. According to the small number of samples, it cannot be stated with certainty that the refractive index in these corneas is higher than the standard corneal refractive index. But it seems that increased stress and strain caused by ectasia can lead to more tightly packed collagen lamellae and changes in refractive index 14 .
Determine the corneal power and ioL power of the eyes. The corneal power of each eye was calculated by the use of keratometry values measured by IOLMaster and keratometry values obtained from the customized keratometric index. Moreover, corneal power was calculated by the use of our proposed method in different areas of 0.5, 1, 2, 3 and 4 mm for normal and customized refractive index of each eye. The corneal power values from different methods were entered into the computational formulas for determining the power of the IOL. For the power of implanted IOL for patients, the predicted refraction value of each formula was determined. Then, the PE was calculated based on the refraction after surgery and the predicted refraction of each formula. Table 5 provides PE values for different values of corneal power and different formulas. As you can see, the calculated keratometry values with the customized index have fewer errors than the index of 1.3375. However, the minimum error is related to the calculated corneal power based on our proposed method in the 4 mm region by applying the Barrett formula. The calculated power values based on the proposed methodology, led to have six eyes out of seven with the refractive PEs in the range of ± 0.5 D.

Discussion
Although the number of eyes entered into this study is not high, based on to our best knowledge, this is the first study that focuses on determining the refractive index of the cornea in patients with keratoconus. Moreover, this is the only study that tries to evaluate the ray tracing technique and finding a method for determining the corneal  Table 4. The results of optimization customized refractive index of corneas. www.nature.com/scientificreports/ power of these patients in order to obtain a suitable IOL power. In previous studies, it has been stated that choosing the IOL power for keratoconus patients is a challenge, especially in patients with corneal keratometry of more than 48 diopters 11,13 . Therefore, only the moderate and severe keratoconus corneas were studied in this study. Watson et al. reviewed the capabilities for predicting IOL power in patients with KCN. They used actual keratometry for moderate eyes and also some severe eyes. 16 eyes out of 40 moderate eyes were placed within ± 1D from the target refraction. The number of severe eyes in this range has not been reported 6 . Kamiya et al. estimated corneal refractive power in the region of 2.4 mm based on the standard keratometry index of 1.3375. When this power was entered into different computational formulas, 8 eyes out of 25 moderate eyes and 0 eyes out of 7 severe eyes were within ± 1D 11 . Savini et al. also used the resulted keratometry from different optical biometry devices measurements and stated that, for different IOL power calculation formulas, 1 to 6 eyes out of 13 eyes of stage II and 0 eyes out of 7 eyes of stage III were in the range of ± 1D. It should be mentioned that stage III in this article has been described eyes with an average keratometry of over 53 diopters and those with keratometry values over 55 diopters have not been included in the study 13 .  www.nature.com/scientificreports/ The achievements from previous studies 6,11,13,15 show that the actual keratometry values usually result in overestimation of corneal power, which leads to underestimation in IOL power and finally terminate with high hyperopic shift especially in eyes with severe and posterior KCN. It can be said that one of the most important causes of incorrect IOL calculation in these eyes is that this process is usually done by measuring the anterior part of the cornea and estimation of total corneal refractive power is based on the ratio between the anterior and posterior parts of the cornea. It is obvious that such a ratio can be changed due to corneal ectasia process 9,13,16 . Rojas et al. 9 calculated the average keratometry index of 30 eyes with KCN, stated that this index was significantly higher in subjects with KCN than normal subjects. It should be noted that the keratometry index of normal people in this study was also 1.3282, which is different from the standard number of 1.3375. The customized KI calculated in our study was in the range of 1.3210 to 1.3307. The calculation of corneal power by the use of customized KI showed that in comparison with using standard KI, this method reduced the value of PE but did not lead to an acceptable level of error reduction in all eyes.

S1-OD S2-OD S2-OS S3-OS S4-OS S5-OD S5-OS
In addition to assessing the customized keratometry index, we have focused on another method to determine the corneal refractive power. In this method, corneal power was calculated by considering both the anterior and posterior parts of the cornea and by implementing the ray tracing technique. We determined corneal power in two different conditions, including customized refractive index and specific corneal refractive index (1.376), in different areas of the cornea including 0.5, 1, 2, 3 and 4 mm distances from the corneas' apex. The customized refractive index in each of these areas did not reduce the error to an acceptable level despite expectations. But when the calculated power of the ray tracing method with the refractive index of 1.376 in the 4 mm region was entered into the Barrett formula, PE of all six eyes with KCN dropped dramatically in the range of ± 1D. Consequently, the PE of our two severe eyes in our study, meaning the left eye of our second and third subjects, was calculated as 0.81 and − 0.09, respectively. However, in previous studies, entering the resulted keratometry values from the anterior corneal segment measurement in any IOL computational formulas failed to detect the refractive PE of severe eyes within ± 1D 11,13 . According to the results of our study, the error of our first subject's eye with PMD was the only one that was not similar to other eyes and it did not place within the range of one diopter. Since the toric lens had been implanted for our first subject, we also used the Barrett toric formula in order to assess the prediction of this computational formula in condition of toric lens implantation. Entering the corneal power of the 4 mm region in this formula, led PE to be zero. A study on patients with PMD and cataract eyes, believes that the refractive outcomes of these cases are highly accurate after cataract surgery 17 and this is because of the paracentral ectasia region in these corneas, which has less effect on the keratometry values of visual axis zone and thus the IOL is calculated with higher accuracy 18 . However, the PMD eye of our study still has a post-operative hyperopic shift, which requires further studies in this area.
In previous studies, suggestions have been made to determine the best formula for calculating the IOL power in patients with keratoconus, but there is no similar methodology in them. Therefore, it is not possible to compare their results and reach a specific rule for power determination. Leccisotti et al. 19 have measured AL in their study by ultrasound, but did not specify whether they used the contact or immersion biometry method. Also, the method used to select the constant values of formulas has not been expressed. Watson et al. 6 have used the optical biometry for the most eyes, the contact ultrasound biometry was used for the rest of the eyes and for severe keratoconus eyes, the standard value of 43.25 D was used for corneal power. Moreover, PE is considered as the difference between the target refraction and the measured postoperative refraction, which in such studies target refraction does not have any meaning. Hashemi et al. 20 have used different methods for determining corneal power in their study, but he did not state whether he used optimum constants or not. Alio et al. 18 have calculated the IOL power of the short-length eye with HofferQ formula and long-length eyes with the SRK/T formula but he did not provide the results of both formulas in all eyes. Kamiya et al. 11 have evaluated various formulas for determining the IOL power. By examining the anterior corneal power, they have reported a high hyperopic shift, especially in advanced KCN eyes and eventually suggested the SRK/T formula as the most correct one. They did not examine the Barrett Universal II formula, which is considered to be the most correct formula in normal eyes [21][22][23] . Also, in half of the eyes, the effect of the total corneal power in the 3-mm region was evaluated in the formulas, but they did not state the result by dividing the keratoconus staging in order to determine if this condition contributes to predictive error reduction of the advanced eyes or not. By entering the keratometry of anterior corneal segment and total corneal power into the formulas, the PE of 63% and 55% of eyes, dropped into ± 1 D from target refraction, respectively. They have stated that corneal power is causing a negligible but meaningful myopic shift. Savini et al. 13 has evaluated different formulas, even the Barrett formula, but did not assess the effect of the total corneal power by considering the curvature of both anterior and posterior segments of the cornea. In this study, different values of corneal power by applying standard and customized keratometry index and the posterior corneal segment with standard and customized refractive index in five different formulas were evaluated. Moreover, the basis of calculating the PE was determined as the difference between the postoperative sphere and the predicted refraction of each formula. When we set the postoperative SE as the basis, there was no eye for any of the corneal power in any formula within the prediction error range of ± 1 diopter. But by considering the postoperative sphere, the corneal power of 4 mm region and the Barrett formula the best result was achieved. It seems that if the predicted refraction of each formula is considered to be the estimation of sphere value after surgery, by using our proposed methodology it is possible to select an appropriate IOL power and prevent the occurrence of hyperopic shift.
Finally, despite the few numbers of eyes in this study, it can be said that determining corneal power plays an important role in determination of the correct IOL power in keratoconus patients. Ectasia causes significant changes in corneal geometry and these changes are ignored in standard keratometry and paraxial optics. But in ray tracing, both anterior and posterior surfaces, thickness and asphericity play the role and corneal power can be obtained in different areas. It seems that considering the keratometry index of these patients, results in an overestimation of corneal power, but the resulted power from measuring both posterior and anterior corneal www.nature.com/scientificreports/ surfaces has better outcomes. Another important factor in determining the IOL power is the ELP estimation that each formula somehow tries to predict it. Based on our results from different formulas, it may be possible to say that the Barrett formula uses a more accurate method for estimating ELP. In fact, regarding appropriate corneal power, ELP estimation has been calculated more precise and a more appropriate IOL power has been obtained. More researches are required to assess the impact of corneal power and prove this study's solution for reaching correct IOL power.

Methods
Our Study includes patients which their corneas had PMD or KCN and they had been under cataract surgery and IOL implantation between October 2015 and November 2018. Patients who had surgeries such as CXL, Intracranial ring implantation, corneal transplantation, and refractive surgeries were excluded from the study. Seven eyes from five patients entered into the study. All surgeries were performed using the standard technique of small-incision phacoemulsification without suture. After creating a 5 or 5.25 mm capsulorhexis and phacoemulsification, an IOL is inserted into the capsular bag through the main incision. The inserted IOLs were include: Acrysof IQ SN6AT9 toric (Alcon Laboratories, Inc., Rochester, NY, USA), Envista MX60 (Bausch & Lomb, Inc., Rochester, NY, USA), Ophthalight ML FLEX 2AS (Mehr Davar, Co., Tehran, Iran), Tecnis ZCB00 (Abbott Medical Optics, Santa Ana, CA), PreciSAL 302AC (Millennium Biomedical, Inc., Cal, USA).
At least 3 months after cataract surgery, patients were examined using slit lamp and auto-refractometer (KR.800, Topcon, Tokyo, Japan). Necessary parameters to simulate patients' eyes were measured by using devices such as tomography (Pentacam HR, Optikgerate GmbH, Wetzlar, Germany), aberrometry (iTrace, Tracey Technologies, Houston, TX, USA), optical biometry (OLMaster700, Carl Zeiss Meditec, Jena, German) and AS-OCT (Casia, SS1000, Tomey, Nagoya, Japan). A trained operator with high precision carried out all measurements. This study was approved by the Institutional Review Board at Tehran Science and Research University and conducted in accordance with the tenets of Declaration of Helsinki. Informed consent was obtained from all participants. eye modeling. By using achieved measurements from biometry, tomography, aberrometry and AS-OCT devices and information of implanted IOL, pseudophakic eye model of each patient was created. Optical analysis was performed on created models with the application of ray tracing method and by minimizing the difference between the aberrations resulted from ray tracing analysis and aberrations obtained from the aberrometry device; we tried to calculate the corneal refractive index as the only missing parameter in each patient's eye model. In the following steps, we explain how the model is built and how the simulation process is optimized.
Determine the required distances for model. Axial length (AL) was measured by partial coherence interferometry (PCI) optical biometry 24 . PCI has significantly improved the measurement accuracy, so that it can measure AL compared to ultrasound. The first commercialized device in the market with PCI technique was IOLMaster. For an accurate interpretation of AL, we should mention that the output from the IOLMaster is not the actual length of the eye's optical path. The achieved values from the IOLMaster device against the immersion ultrasound are calibrated by the use of Eq. (1) 25 : Which the Ax Zeiss is the output from the PCI and OPL is the optical path length measured by the PCI. Therefore, the length of the optical path can be calculated by the Eq. (2): The Eq. (3) was used to calculate the eye length from the OPL 25 : To determine the distance between the posterior part of the cornea to the anterior part of the IOL as well as the thickness of the IOL, the AS-OCT device was used. These measured values were used to determine the location of the IOL and its design. customized modeling of patient's cornea. We used the anterior segment imaging system of Pentacam based on Scheimpflug technology. With the help of this technology, the spatial coordinates from multiple points of both anterior and posterior surfaces of the cornea can be obtained. The most reliable data for building the corneal geometrical model is the raw elevation data of the posterior and anterior parts of the cornea that has not been manipulated by any software algorithm in the topography device 26 . Hence, in the first step, the elevation data of the anterior and posterior parts of cornea were extracted from the device in the form of comma-separated values (CSV format). Then, the data was converted to the Cartesian format by using an algorithm, which was written in MATLAB software version 2018 (Mathworks, Natick, USA). The next step is to produce the posterior and anterior surfaces of the cornea from the geometry of points cloud of the previous step, in which these points cloud, demonstrate the Cartesian coordinates of the points scanned by the topographer in a three-dimensional space. The points cloud of both posterior and anterior parts of the cornea were entered into the SolidWorks software Version 2018 (SolidWorks Corp., MA, USA). By the use of this software, a mesh was generated on the points cloud and the best surface that was fitted into the mesh was selected. www.nature.com/scientificreports/ To do this, deviation analysis was used to evaluate the difference between the surface and the mesh and based on the smallest difference, the most accurate and possible surface, was created. This process was performed to create both posterior and anterior surfaces. The created surface was then cut to the size of the cornea and finally the gap between the two surfaces was filled and they were knitted together and became like a solid shape. Therefore, with the help of this software, a solid model was built that represents the actual while customized geometry of each patient's cornea. Figure 1 shows the resulted cornea of the described method. The final model was stored in a format that is compatible with the Zemax software version 2014 (Zemax LLC, WA, USA).
patient's customized ioL modeling. In order to model the implanted IOL for each patient, the specific optical parameters of each IOL provided by the manufacturer companies were used. These parameters included the refractive index and the spherical aberration of the lens and the abbe number (if they exist in the lens specification) 27 . Table 6 shows these specifications for modeled lenses in the study. Since the thickness values of the lenses were not available in the tables provided by manufacturers, the AS-OCT device was used to determine them. For modeling the lenses, the optical properties of each lens including refractive index, spherical aberration, abbe number, thickness, and refractive power of the lens were used and the radius of anterior, posterior and asphericity conic constant were considered unknown. These unknown parameters were determined by considering the conditions of lens production (i.e., considering the refractive index of the aqueous fluid and vitreous around the lens to achieve the required power) and optimization process in Zemax optical design software. To model Acrysof IQ toric lens, due to the toricity of the posterior side of the lens, the unknown parameters of posterior section consist two parameters including curvature radius of flat and steep axes. After determining the missing parameters, the lens was evaluated for the power and aberration (Table 7). exact ray tracing and optimization process. In geometrical optics, it is assumed that the wavelength of the light is small enough and the light emission can be described in terms of the beam. The beam path is determined by reflection and refraction, which is defined by the Snell law. Based on Snell law, tracing a large number of beams with finite size in an optical system with defined pupil size is called exact ray tracing or real ray tracing (RRT), which is a well-known approach to analyze the system's optical characteristics. The use of RRT method leads to a comprehensive description of the optical system, in particular, systems in which there are effects of irregular surfaces. With the help of this approach, any deviations from the ideal optical system can be measured as optical aberrations 27 .
To implement the RRT method, an ophthalmic model must be made firstly. For this purpose, the information obtained from the previous steps, including the cornea (anterior and posterior surfaces), designed IOL and  www.nature.com/scientificreports/ customized interval measurements (including eye length and distance between the cornea and the lens) for each patient were used to build each patient's eye model. In Zemax software, we determined the eyeball using eye length and corneal thickness of the patient. The interior environment of sphere was considered as an environment with refractive index of aqueous humor (1.336 28 ). Then the created file of cornea's surfaces was entered to the software as CAD part, so that all cornea's measured points were fully entered into the eye model. In the next step, using the pupil center coordinates, pupil size and scan size measured by the iTrace device, we considered the precise location and dimensions of input pupil in the model. After defining the pupil, the individual's designed IOL was determined in the software by entering its anterior, posterior and its conic constant curvature radii. We determined the vitreous environment in the ocular model by defining an enclosed environment and tangent to the both posterior surface of the IOL and the surface of retina with vitreous refractive index (1.337 28 ).
After building a pseudophakic ophthalmic model, it is time to define the aberrations measured by the iTrace device. By studying the order of considering the Zernike coefficients in the iTrace device and checking the standard Zernike coefficients defined in the software, each measured coefficient was compared with its counterpart coefficient in the software. The values of 18 Zernike coefficients from the order of 2 to 5 were entered into the optical design software with the definition of merit function. The coefficients of order 0 and 1, meaning tilt and piston, were eliminated because the piston is only a constant offset and shows absolute distance from the object. The tilt also places the image according to the fovea center, which can be compensated by adjusting the eye relative to the fixation angle 27 . Then we implemented the RRT technique to simulate the process of light propagation and its passage through the optical system of the eye and its arrival at the retina. In this simulation, a wavelength of 785 nm was used to achieve eye aberrations. The optimization approach using Damped Least-Squares method 29 was used to adapt the aberrations derived from the RRT method with predefined merit functions. In fact, by minimizing the difference between the values of aberrations derived from RRT and the measured ones from the iTrace device, by applying the optimization approach, we tried to achieve value of the refractive index of the cornea as an unknown parameter in the eye model. corneal power measurement. The eyes that were included in the study had Pentacam tomography measurements before cataract surgery. As the method described in the related section of customized corneal modeling of patients, their corneal surfaces were built based on the images taken before surgery and were entered into the optical analysis software. This time, we considered the customized model of patient's cornea that was followed by aqueous humor and with dimensions larger than the individual's eye length, and the rest of the elements were deleted. Light passes through the cornea and enters to the aqueous humor, after that it goes forward until it focuses in an area and then diverges.
In order to calculate the minimum and maximum corneal power and their corresponding axes, we need to determine the occurrence location of the focusing spot minimum diameter. The minimum and maximum power of the cornea refer to the places that the focus of beams has the maximum and minimum distances from the cornea, respectively. To find these locations, we first determined the location of beam's focal point in the radial form. In the radial form, the location where the diameter of the spot is minimal in both x and y directions is specified as the focal point. Then we moved the image screen once to 5 mm before this location (to find the www.nature.com/scientificreports/ maximum power) and once after this location (to find the minimum power). We defined the merit function to find the spot's smallest diameter and the tilt of cornea around the z-axis considered unknown. By implementing the Ray Tracing method and optimization process, the cornea was rotated around the Z-axis and the defined merit function was optimized. Thus, the angles for which the smallest diameter of focal spot is placed at the minimum and maximum distance from the cornea, were determined. Then, at the specified angles, we calculated the actual distance from the focus of light to the posterior surface of the cornea. Finally, the refractive index of 1.336 divided to the distance between the focal point and the posterior surface of the cornea, determined the amount of corneal power (Power (D) = 1.336 / f (meters)). Inserting the minimum and maximum distance from the corneal surface as the f value in this formula resulted in calculating the maximum and minimum corneal power.
The described process was repeated in areas with diameters of 0.5, 1, 2, 3 and 4 mm from the cornea, in order to calculate the corneal power at different intervals from the cornea's apex. In addition, all of the steps mentioned above were performed in two different modes in order to determine the corneal power; once by taking into account 1.376 as the refractive index of the cornea and another time by considering each individual's eye refractive index obtained from the previous stage.
calculating the ioL power. In the previous section, the total power of the cornea was calculated in different regions. In this section, we used the calculated values of corneal power and other preoperative measurements including eye length, anterior chamber depth (ACD) (from epithelium to lens), crystalline-lens thickness, and WTW in order to determine the IOL power. These values were entered into the IOLMaster 700 device to calculate the IOL power by using the SRK/T, Holladay 2, HofferQ, and Haigis formulas. Moreover, the Asia Pacific Association Cataract and Refractive Surgeons website was used for the Barrett Universal II formula v1.05. It should be mentioned that because of the implanted toric lens in the first patient, the Barrett Toric Calculator formula v2.0 was also evaluated for that patient. For the constants of each IOL, the optimized constants were also used which is available on the User Group for Laser Interference Biometry (ULIB) website. For the Ophthalight lens, there were no optimal constants available on this website, therefore we used the constants provided by the manufacturer.
To evaluate the effect of total power of the cornea on the IOL power calculation, the anterior keratometry values that had been measured by the IOLMaster device before surgery were entered into the formulas. These keratometry values were calculated by the use of keratometric index (KI) of 1.3375. Since it has been stated that the values of KI differ in the patients with keratoconus 9 , we decided to also calculate the KI for each eye. According to Eq. (4) 9 , we measured KI by the help of true net power and anterior radius obtained from the Pentacam results before surgery. Then, using the customized KI of each eye, we converted the measured radius values from IOLMaster into keratometry numbers in terms of diopter (Eq. (5)). www.nature.com/scientificreports/ In formula 5, r is the radius measured in meter and K is the keratometry value in diopter. The keratometry numbers derived from customized KIs are also included in the IOL power calculation formulas.
Patients were examined after at least 3 months from cataract surgery until their refraction riches to a stable condition. Subjective refraction outcome was used to evaluate the optical results in patients. The distance between the patient's eye and Snellen chart was about 20 feet (6 m). To evaluate each formula, the Prediction Error (PE) was calculated. This error was obtained by subtracting the predicted refraction of each formula (based on the power of implanted IOL for patient) from the actual refraction. The negative PE value indicates that the resulting refraction is more myopic than the predicted refraction and the positive PE value shows a more hyperopic result. Figure 2 demonstrates a general overview of implementing our proposed methodology.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.