Considerations on the Castrop formula for calculation of intraocular lens power

Background To explain the concept of the Castrop lens power calculation formula and show the application and results from a large dataset compared to classical formulae. Methods The Castrop vergence formula is based on a pseudophakic model eye with 4 refractive surfaces. This was compared against the SRKT, Hoffer-Q, Holladay1, simplified Haigis with 1 optimized constant and Haigis formula with 3 optimized constants. A large dataset of preoperative biometric values, lens power data and postoperative refraction data was split into training and test sets. The training data were used for formula constant optimization, and the test data for cross-validation. Constant optimization was performed for all formulae using nonlinear optimization, minimising root mean squared prediction error. Results The constants for all formulae were derived with the Levenberg-Marquardt algorithm. Applying these constants to the test data, the Castrop formula showed a slightly better performance compared to the classical formulae in terms of prediction error and absolute prediction error. Using the Castrop formula, the standard deviation of the prediction error was lowest at 0.45 dpt, and 95% of all eyes in the test data were within the limit of 0.9 dpt of prediction error. Conclusion The calculation concept of the Castrop formula and one potential option for optimization of the 3 Castrop formula constants (C, H, and R) are presented. In a large dataset of 1452 data points the performance of the Castrop formula was slightly superior to the respective results of the classical formulae such as SRKT, Hoffer-Q, Holladay1 or Haigis.


Introduction
Since the introduction of biometry instruments for measuring axial eye length and publication of the first formulae by Fyodorov in 1967 [1] and by Gernet in 1970 [2], many calculation schemes have been proposed for predicting the power of an artificial lens in cataract surgery. The simplest are the SRK and SRK2 formulae [3], which use a linear model to derive the emmetropic lens power from the axial length measurement, the 'keratometry reading' and the A constant used as an intercept. In the SRK2, the A constant is modified by an offset depending on the axial length. The classical theoretical-optical formulae in use today are the SRKT [4,5], Hoffer-Q [6][7][8], Holladay1 [9], and the Haigis formula [3]. The latter is sometimes used in a simplified form with standard values for constants a1 = 0.4 and a2 = 0.1 and optimization of a0 only, or in a version with triple constant optimization a0/a1/a2 which appears to yield superior results in long and short eyes.
In addition to these formulae, published in the 1980s and 1990s, many other calculation concepts have been proposed. Most of these are again based on calculations with paraxial simplifications, and most of them have never been published such as the Holladay2 or the Barrett Universal II formula. In addition, several formulae have been developed to consider the situation after keratorefractive surgery (e.g. the Barrett TrueK, Haigis-L, Masket, Shammas formula or others), where most of the classical formulae fail. Some other calculation strategies are based on full aperture raytracing without simplifications to linear Gaussian optics (e.g. Okulix or PhakoOptics) or on artificial intelligence (e.g. Hill RBF calculator or PEARL formula), where a large dataset replaces an anatomical or optical model [10]. All of these concepts claim to yield superior results to other formulae.
In general, it is difficult to check the performance of different formulae, as it is mostly dependent on the composition of the dataset, and there are no common rules on how to evaluate the performance. Therefore, some formulae are superior in long or short eyes [11,12], or in eyes with a steeper or flatter cornea, or selected combinations of axial length and corneal curvature. Additionally, the performance of a formula mostly depends on the dataset used for constant optimization as well as the technique applied for optimizing the constants, and it can also depend on the biometer settings, quality of the clinical measurement data (e.g. refraction and lane distance), and the labeling tolerance of the lens [13][14][15]. As target criteria for evaluating the performance of a specific formula we use the prediction error (PE) in terms of achieved refraction after cataract surgery (SE) minus the formula prediction (predSE), the absolute value of the prediction error (absPE), or the portion of cases within specific refraction limits (e.g. �±0.25, �±0.5, or �±1.0 dpt) [16]. With the PE we mostly obtain the systematic deviation of the achieved refraction from the formula prediction, e.g. if a study population systematically results in hyperopia or myopia. In contrast, the absPE mostly refers to the variation or scatter of the results, but a proper interpretation of absPE as variation of the results presupposes that the mean or median PE is close to zero [14]. The PE is typically a two-sided distribution, and if close to a normal distribution would be fully characterized by the mean and standard deviation. In contrast, the absPE is a strictly one-sided distribution, where the mean and standard deviation do not have a direct meaning, and in this case the median or confidence levels should be used for characterization [16].
Measurements in previous generations of biometers were restricted to axial length, corneal front surface curvature, anterior chamber depth, and optionally the horizontal corneal diameter and the formulae of that era focused on the parameters available at that time [3]. Most modern biometers provide additional information about corneal thickness, corneal back surface curvature, aqueous depth, central thickness of the crystalline lens, and vitreous depth [17]. Furthermore, some biometers analyze corneal topography or even tomography data at many points, providing information on the asphericity in addition to the curvature of the front or front and back surfaces. In the future, with consistent tomographic data available from most biometers, there might be a quantum leap towards full aperture raytracing in routine intraocular lens power calculations.
The purpose of the present study is: • to introduce and explain the principle and calculation scheme of the Castrop formula; a paraxial vergence formula for intraocular lens power calculation, • to assess the performance of the Castrop formula with an explorative data analysis applied to a large dataset of a cataract population with preoperative biometric data, the power of the implanted lens, and postoperative refraction data using cross validation techniques, and • to compare the results of the Castrop formula with the respective results of classical formulae such as SRKT, Hoffer-Q, Holladay1, and the Haigis formula.

The Castrop IOL calculation concept
The Castrop formula is a paraxial vergence formula based on a pseudophakic model eye with 4 refractive surfaces and 3 formula constants (C, H, and R): a refractive correction at the spectacle plane, a thick cornea with front and back surfaces, and a thin intraocular lens (IOL).
Castrop-Rauxel is a city in the middle west of Germany in the federal state of North Rhine Westphalia. The respective refractive powers of these elements are: predSE for the predicted spherical equivalent refraction, P CA = (n C -1)/R CA and P CP = (n-n C )/R CP for the corneal front and back surface power, and P IOL for the lens implant. R CA and R CP refer to the corneal front and back surface radii of curvature, the vertex distance d VD refers to the interspace between the spectacle plane and the corneal front vertex plane, d C to the central corneal thickness, d AQD to the aqueous depth, and d V to the vitreous, respectively. For the refractive index we used n C = 1.376 (cornea), and n = 1.336 for aqueous humour and vitreous humour derived from a common schematic model eye, respectively [17].
To account for the effect of refraction measurement lane distance or systematic error in refraction measurement, an offset of R was superimposed on the refraction term ta achieve predSE. The IOL is assumed to be located within the envelope of the crystalline lens, similar to the recommendation of the lens haptic plane concept of Norrby [18] or the Olsen formula [19,20]: The effective lens position ELP, defined as the distance of the corneal front vertex plane to the plane of the IOL, is derived from the preoperative phakic anterior chamber depth (ACD), the central thickness of the crystalline lens (LT), and an axial offset (H) as: where C refers to the (anterior) portion of the thickness of the crystalline lens LT to the lens equator. C is typically in a range of 0.36 to 0.44.
With AL as a measure for the axial length of the eye, the vergence formula for the IOL power reads: In total the Castrop formula is characterized with 3 formula constants, C, H, and R. Solving the equation for the predicted refraction predSE yields: In cases where measurements of R CP or d C are not available or not reliable, the respective data can be derived from a schematic model eye [17] maintaining a fixed ratio of corneal front to back surface curvature:

Formula constant optimization
There are different potential strategies for optimizing C, H and R. The constants could be optimized in turn (e.g. first the C, then the H, and then R), or en bloc using an enhanced optimization strategy as commonly used in engineering. As the effects of the 3 constants are not independent (similar to the Haigis 3 constant optimization) we decided in this study to use a nonlinear iterative optimization strategy. As target measures for optimization we considered the mean absolute prediction error (mean absPE) and the root mean squared PE (rmsPE). The Levenberg-Marquardt algorithm [21,22] was initialized with the formula constant triplet C/ H/R = 0.4/0/0 with boundary conditions 0.35 to 0.45/-0.8 to 0.8, and -0.8 to 0.8. We limited iterations to a maximum of 100, with a damping of 1e-2, a termination step size of 1e-14 (optimization ends when last step is smaller than 1e-14), and a function tolerance of 1e-16 (optimization ends if rmsPE improvement is smaller than 1e-16).

Measurement data
In this retrospective study we analysed a dataset with 1601 clinical data points from a cataract population from Augencentrum Rosenheim which was transferred to us (802 left eyes and 799 right eyes; 961 female and 640 male). Mean age was 74.2±8.1 years, median: 75 years, range: 47 to 97 years). According to the statement from the 'Bayerische Landesärztekammer', an ethics approval or patient informed consent was not required for this study. This study is a retrospective evaluation of data which were collected during routine examinations. No extra examinations or measurements were performed.
The data were transferred to us in an anonymized fashion, which precludes back-tracing of the patient. The anonymized data contained preoperative biometric data derived with the IOL-Master 700 (Carl-Zeiss-Meditec, Jena, Germany) including axial length (AL), central corneal thickness d C , external phakic anterior chamber depth ACD measured from the corneal front apex to the anterior apex of the crystalline lens, lens thickness LT, corneal front surface radius measured in the flat (R1) and in the steep meridian (R2). In all cases a Sensar 1 piece intraocular lens (Johnson & Johnson, Brunswick, USA) was inserted. In addition to the refractive power of the inserted lens (P IOL ), the postoperative refraction (sphere and cylinder) 6 to 8 weeks after cataract surgery was recorded. We restricted the dataset to postoperative data with a visual acuity of 0.6 or higher to ensure that the postoperative refraction was reliable. It is shown in the literature that with a low visual acuity e.g. due to retinal pathologies refractometry might be unreliable [16]. From the original dataset, a subset of N = 1452 complete data points were used for formula constant optimization. The spherical equivalent of postoperative refraction (SE) was calculated as sphere + ½�cylinder, and the mean corneal front surface radius was calculated as R CA = ½ (R1+R2). The corneal back surface curvature was calculated from the corneal front surface curvature as shown above using a fixed ratio of front to back surface radius. Axial length was pre-processed using the Cooke correction [23,24]. The descriptive data on pre-cataract biometry, P IOL and postoperative refraction are summarized in Table 1.

Calculation strategy
The anonymized Excel data (.xlsx-format) were imported into MATLAB (Matlab version 2019b, The Math Works, Natick, USA) for further processing. For the Castrop formula we implemented the calculation strategy described above. To compare the formula performance to classical formulae we also implemented the SRKT formula [4,5], the Hoffer-Q formula [6][7][8], the Holladay1 formula [9], the simplified Haigis with one formula constant a0 and standard values for a1 = 0.4 and a2 = 0.1, and the Haigis formula [3] with optimized constant triplet a0/ a1/a2 [13]. For all calculations, the deviation of the achieved postoperative refraction SE from the formula predicted refraction predSE was quoted as the prediction error PE. The absPE refers to the absolute value of PE [16].
For cross validation, the N = 1452 measurements were split randomly into a training set (70%, 1017 cases) and a test set (30%, 435 cases) [25]. For the training set, the constants for the Castrop formula (C, H, R), the SRKT (A), the Hoffer-q (pACD), the Holladay1 (SF), the simplified Haigis (a0) and Haigis formula with triplet optimization (a0/a1/a2) were optimized in terms of minimizing the rmsPE using the nonlinear Levenberg-Marquardt algorithm [21,22]. In the next step we calculated the mean and median absolute prediction error, and the root mean squared prediction error for all lens calculation formulae based on the optimized formula constants using the test dataset. Then for all formulae we derived the fraction of cases with absPE � 0.25 dpt, � 0.5 dpt, � 1.0 dpt, and � 2.0 dpt.
As an overall metric for formula performance we calculated the formula performance index FPI as described by Hoffer & Savini [16]. For estimation of the slope of the correlation of prediction error PE vs. axial length AL as one component of this index, we used a robust estimator with Studentized residuals.
The Wilcoxon signed rank test was used to compare the absPE of the Castrop formula to the absPE of other formulae. A p-value of <0.05 was considered statistically significant.  Table 2 summarizes the descriptive statistics for the prediction error PE for all formulae. The prediction error is derived using the constants optimized on the training data and crossvalidated on the N = 435 test data. PE is listed in terms of mean, standard deviation, median, range, and 5% / 95% quantiles. Fig 1 displays the scatter of the formula prediction error PE for the test data for all 6 formulae under test, together with the approximated kernel distribution, the median, the quartiles and the 90% confidence interval. Table 3 summarizes the descriptive statistics for the absolute value of the prediction error absPE for all formulae. The absolute value of the prediction error is calculated using the constants optimized on the training data and cross-validated on the N = 435 test data. absPE is listed in terms of mean, standard deviation, median, range, and 5% / 95% quantiles. Fig 2 displays the scatter of the absolute value of the formula prediction error absPE for the test data for all 6 formulae under test, together with the approximated kernel distribution, the   . The upper part of the plot shows the approximated kernel density distribution for the 6 formulae under test, and the lower part the scatter of the PE together with the median (circle), the quartiles (black line), and the 90% confidence interval (5% and 95% quantile, blue line). We see that after optimization of the formula constants on the training data none of the formulae shows a systematic offset in the prediction error. The blue lines in the lower plot representing the 90% confidence interval show that the variation on the test data is largest with the Hoffer-Q formula and lowest with the Castrop formula.

Results
https://doi.org/10.1371/journal.pone.0252102.g001 The formula constants were optimized on the training dataset and cross validated on the test set.

PLOS ONE
https://doi.org/10.1371/journal.pone.0252102.t003 In Fig 4 the ratio of cases with an absolute prediction error within limits of �0.25 dpt, �0.5 dpt, and �1.0 dpt is shown for the SRKT, the Hoffer-Q, the Holladay1, the simplified Haigis, the Haigis, and the Castrop formula. The bars correspond to the results of the test data, and the constants were optimized for the root mean squared prediction error rmsPE on the training data.

Discussion
Numerous formulae for calculation of intraocular lens power have been proposed in the last 20 years. In contrast to the basic formulae of Fyodorov [1] or Gernet [2] or the classical  [3][4][5][6][7][8][9], most of the formula authors nowadays do not disclose or publish the calculation strategy. At best they offer web based applications or software solutions for calculating the lenses. Such software tools do not allow batch calculations on a large set of patient data. Today, in many countries of the world classical formulae are increasingly being replaced by 'modern' calculation strategies such as the Barrett Universal II, Kane, Pearl, EVO, VRF or T2 formula [10,11,26]. To compare the prediction performance with other formulae it is necessary to enter the data manually, introducing a large risk of transcription errors. Additionally, a systematic optimization of constants is not possible for undisclosed formulae. Reading out the appropriate formula constant is mostly trial and error e.g. to eliminate the mean or median prediction error by varying the constants in the calculation scheme.
Withholding disclosure of formulae gives authors the option to adjust internal factors and offsets over the years in order to enhance the formula performance with more and more clinical data. This however bears the risk of different versions of the software yielding different outputs being available.
In the scientific world it is most appropriate to publish a new formula and to outline all options for adjusting the calculation. Such a publication includes the calculation strategy for the intraocular lens power, all the input parameters required for the calculations, and the adjustment parameters such as formula constants. Last but not least the formula author should advise how to optimize the formula constants based on a clinical dataset. But even in the classical formulae-with the exception of the SRKT formula-there is no guidance on how to optimize the formula constants. And in the original publication of the SRKT formula, the instruction is mostly based on a simplified back-calculation of the A constant similar to the empirical SRK2 formula.
In the present paper we describe the strategy of intraocular lens power calculation using the Castrop formula as described in Materials and Methods. The large dataset of 1452 eyes was split into training and test subsets [25]. The training data were used for formula constant optimization, and the test data for cross validation. In contrast, using the same dataset for constant optimization and validation risks 'overfitting' which may lead to an overestimation of an algorithm's performance. Our results as presented in an explorative data analysis show that the Castrop formula yielded a slightly better performance compared to the classical SRKT, Hoffer-Q, Holladay1, or Haigis formula in terms of prediction error and absolute prediction error. In Fig 1 we see that the kernel distribution approximated to the prediction error shows the highest peak and a better concentration of the data with the Castrop formula. In contrast the Hoffer-Q formula seems to have a poorer performance. However, the distribution of all formulae is more or less symmetric about a prediction error of 0. In Table 2 the respective numbers are listed; in this table the mean describes any systematic offset in the formula predicted refraction and the standard deviation gives us a measure for the scatter of the individual values for the prediction error. With the absolute prediction error shown in Fig 2 the outcome is very similar: with the Castrop formula the entire distribution is closer to 0 which indicates a good performance compared to other formulae. The peak of the kernel distribution approximated to the absolute prediction error is larger and the distribution appears slimmer compared to the other formulae. In Table 3 the respective descriptive values are listed. As the distribution of the absolute prediction error values is strictly one-sided the mean and standard deviation are of minor interest. In contrast, the median value and the 95% confidence level are more relevant for interpretation of the results. From the data we find that 95% of all cases in the test data are within limits of 0.9 dpt of absolute prediction error. It is unfortunate that the Castrop formula cannot be compared directly with modern formulae as these formulae are not disclosed and therefore even with manual data entry in calculation software available online the standards for optimizing the formula constants are different.

In conclusion
This study describes the calculation principle of the Castrop vergence formula for calculation of intraocular lens power based on biometric data of the eye. The concept is based on a pseudophakic model eye with 4 refractive surfaces: a refractive correction at spectacle plane, a meniscus lens for the cornea described with corneal front and back surface curvature and central thickness, and an intraocular lens implant. If tomographic data for the corneal back surface curvature and central corneal thickness are available and reliable, they can be directly used in the formula. Where they are not available or not reliable, the corneal back surface is derived from a fixed anterior to posterior corneal curvature ratio taken from the Liou-Brennan model eye. The Castrop formula uses 3 constants: the C constant in accordance with the Olsen formula, the Offset value H for a systematic shift of the IOL plane (which is mostly characterized by the optics and haptics design), and R for an offset in predicted refraction.
The Castrop formula together with 5 other classical lens power calculation formulae was applied to a large dataset of 1452 clinical measurements treated with one lens model. The dataset was split into a training set for constant optimization and a test set for cross-validation. Constant optimization for all formulae was performed using nonlinear optimization techniques aiming for a minimum root mean squared prediction error. In our explorative data analysis, the Castrop formula shows slightly superior results compared to the classical formulae in terms of prediction error and absolute prediction error. Further investigation of this calculation concept is needed with different datasets of different lenses to assess the advantages and limitations of this calculation concept.