Key points The corneal back surface is a vital parameter for optical calculations, especially in the context of lens power calculation for toric implants or following refractive corneal surgery. The equivalent power and astigmatic vector components of the corneal back surface can be predicted using a multilinear regression model based on front surface measurements and measurements from optical biometers. Deep learning algorithms outperform multilinear regression models in predicting corneal back surface power, even in the case of crossed cylinders between both surfaces. BACKGROUND

Over the last two decades, research has investigated the effect of corneal back surface power on the imaging properties of the eye.1, 2 Especially in biometry and intraocular lens power calculation before cataract surgery, the corneal back surface is widely discussed as a relevant parameter, as well as after keratorefractive surgery.3-5 Modern anterior segment tomographers provide reliable data on the corneal front and back surface curvature1, 6-8 and central corneal thickness (CCT) as well as the dimensions of the anterior chamber and crystalline lens which is used in several newer generation lens power calculation strategies. However, most biometers such as the IOLMaster 500 (Carl-Zeiss-Meditec, zeiss.com), the LenStar 900 (Haag-Streit, haag-streit.com) or the OA-2000 (Tomey, tomey.com) do not provide measurements of the corneal back surface.

Especially in the context of toric lens calculation, it is widely discussed that the corneal back surface curvature does not follow a fixed ratio of anterior to posterior curvature, as assumed by keratometers which measure corneal front surface curvature and convert this to corneal power using an assumed refractive index.3, 4, 9-12 The corneal back surface astigmatism shows some against the rule astigmatism which has to be considered during toric lens power calculations in order to yield the correct refractive results after surgery.13

With the classical Javal's rule14, 15 we already know that there exists some discrepancy between keratometric astigmatism and the total ocular astigmatism, which could be due to lenticular or corneal back surface astigmatism. For clinical purposes, Javal's rule is generalised with Grosvenor's simplification.16 Whereas lenticular astigmatism cannot be currently measured in vivo with sufficient accuracy, high-resolution anterior segment tomographers, mostly based on swept source optical coherence tomography (SS-OCT) do quantify corneal back surface curvature.10, 12, 17 In the last decade, several correction strategies have been proposed to correct keratometric corneal power measurement. Most are based on a comparison of corneal astigmatism and the refractive cylinder after implantation of rotationally symmetric lenses (e.g., the Abulafia-Koch correction9) or the correction strategies proposed by Savini12 or Tutchenko,18 while others use tomographic data from Scheimpflug systems to derive statistical corrections in situations where the corneal back surface cannot be derived.19 All decompose the astigmatism into vector components using trigonometric functions (e.g., Humphrey notation20, 21) and overlay fixed or regression-based offsets to the components to reduce or eliminate the overall trend error.

With LaHood et al.’s17total keratometry (TK) module for the IOLMaster 700 (Carl-Zeiss-Meditec, zeiss.com), which was launched to the European market in 2018, corneal back surface data is captured with a scanning SS-OCT. These data are binned to the corneal front surface data which are measured in classical fashion with telecentric keratometry to maintain distance independence of the measurement.

The purpose of the present study was to use a large dataset of biometric measurements with the IOLMaster 700 in a cataractous population plus the back surface measurements (TK module) to develop and train a deep learning algorithm for prediction of corneal back surface curvature from corneal front surface curvature and other biometric measures such as patient age, axial length, corneal thickness, anterior chamber depth, lens thickness and the horizontal corneal diameter. This prediction model, which could be used for standard biometers where no corneal back surface measurement is available, will be validated and compared to a multivariate linear regression model based on corneal front surface curvature data.

METHODS Dataset for the deep learning algorithm

In total, a dataset with 48,455 measurements from the IOLMaster 700 from two clinical centres (Augenklinik Castrop, Germany and Department of Ophthalmology and Optometry, Kepler-University, Linz, Austria) was considered for this retrospective study. All biometry measurements were performed prior to cataract surgery. This study was registered at the local Ethics Committee (Ethikkommission der Ärztekammer des Saarlandes with the registration number 157/21). The data were anonymised by the source and transferred to a.csv data table using the software module for batch data export. The tables from both clinical centres were merged. Data tables were reduced to the relevant parameters required for our data analysis and finally contained measurements of corneal front surface curvature (flat meridian R1 with axis A1; steep meridian R2 with axis A2), back surface curvature (flat meridian RP1 with axis AP1; steep meridian RP2 with axis AP2), axial length (AL), CCT, phakic anterior chamber depth (ACD) measured from the corneal front apex to the lens front apex, phakic lens thickness (LT), horizontal corneal diameter (W2W) and patient age.

Duplicate measurements of eyes, missing data or data with a ‘Failed’ or ‘Warning’ in the internal quality check of the IOLMaster 700 for keratometry, AL, CCT, ACD, LT or W2W were excluded. After checking for ‘Successful’ measurement for corneal front and back surface curvature data a dataset of N = 19,553 measurements was used for training, validation and test of our prediction algorithm. The data were transferred to MATLAB v.2019b (MathWorks, mathworks.com) for further processing.

Pre-processing of the data Custom software was written in MATLAB 2019b. From the corneal front surface curvature data (R1/A1 and R2/A2) measured with the IOLMaster 700 the three vector components CFSPMEQ, CFSPM0° and CFSPM45° were derived using vector decomposition according to Alpins20 For the measured corneal back surface curvature the respective 3 vector components CBSPMEQ, CBSPM0° and CBSPM45° were derived using with the refractive index nC = 1.376/nA = 1.336 for the cornea and aqueous humour, respectively.22 Setup of the prediction algorithm

To identify the relevant input parameters, a stepwise fit strategy was applied. With this procedure, a stepwise linear regression model is fitted to the input data (corneal front surface power [CFSPMEQ, CFSPM0° and CFSPM45°]), AL, CCT, ACD, LT, W2W and patient age to read out the impact of the potential effect sizes. This stepwise fit begins with an initial constant model and takes forward or backward steps to add or remove variables, until a stopping criterion (maximum number of iterations: 12, significance level for a variable to be entered: 0.05) is satisfied. ACD and lens thickness were not considered as significant effect sizes in this stepwise fit model, therefore both variables were discarded.

A feedforward shallow multi-layer multi-output neural network23, 24 was set up for predicting the corneal back surface power (output: CBSPPEQ, CBSPP0° and CBSPP45°) from the input parameters corneal front surface power (CFSPMEQ, CFSPM0° and CFSPM45°), AL, CCT, W2W and patient age. For the training function we used the Levenberg-Marquardt algorithm,25, 26 as this algorithm is known to provide good performance in terms of convergence and stability. Based on the squared prediction error derived from the three vector components of the measured and the predicted back surface power values optimisation was performed in terms of minimizing the mean squared prediction error, which refers to a metric for performance P of the prediction:

In a preconditioning step we tested several options using the entire dataset with several choices for the number of hidden layers (1/2/3/4 hidden layers) and 1 to 25 neurons per layer with a maximum number of 50 iterations (epochs) for optimisation in order to form an impression of the required number of layers and the number of neurons in our perceptron. Figure 1 shows the performance of the feedforward network for some selected combinations of hidden layers 1 to 4 and number of neurons per layer of 1, 2, 4, 8, and 16 using the entire dataset for training with the cumulative histograms of the squared prediction errors. From the graphs we see that the squared prediction error sPE becomes systematically smaller with an increasing number of neurons and increasing the number of hidden layers. However, in contrast, the complexity of the network increases rapidly with the number of hidden layers. As we noticed a large improvement only from 1 to 2 hidden layers and a very small improvement from 2 to 3 or 3 to 4 hidden layers, we decided to use 2 hidden layers in our final model. For the number of neurons per layer the behaviour is similar; there is a large improvement from 1 to 2, from 2 to 4 and from 4 to 8 neurons, but with increasing numbers of neurons the performance only improves marginally. However, an increasing number of neurons per layer makes the perceptron structure much more complex. Therefore, we decided to use 12 neurons per hidden layer in our final model as a compromise between performance and complexity.

The result of the preconditioning to determine the proper number of hidden layers and number of neurons per layer. The figure displays a series of cumulative histograms of the squared prediction error for variations of number of neurons per layer (to the right, 1st column: 1 neuron per layer, 2nd column: 2 neurons per layer, 3rd column: 4 neurons per layer, 4th column: 8 neurons per layer, and 5th column: 16 neurons per layer) and number of hidden layers (1st row: 1 hidden layer to 4th row: 4 hidden layers). All histograms consider a range for the squared prediction error between 0 and 0.05 dpt2. For comparability, all calculations in this preconditioning step are performed with 50 epochs

In the next step, the dataset with N = 19,553 measurements was divided using a random selection into a training set (70%, N = 13,687), a validation set (N = 2933) and a test set (N = 2933). 23 The neural network was trained with the training dataset and later validated with the validation dataset. A final check was performed using the N = 2933 test dataset.

For reference we defined the corneal back surface curvature from a fixed ratio of front to back surface radius based on the Gullstrand schematic model eye (7.7 mm/6.8 mm) with the three vector components CBSPGEQ = 6.8/7.7·CFSPMEQ, CBSPG0° = 6.8/7.7·CFSPM0° and CBSPG45° = 6.8/7.7·CFSPM45°, respectively. In addition, we defined a multivariate linear regression model (with a 3×3 transformation matrix and a 3×1 intercept vector) which uses the three vector components CFSPMEQ, CFSPM0° and CFSPM45° as input values to predict the three vector components CBSPmvEQ, CBSPmv0° and CBSPmv45° as output parameters. For a direct comparison of the prediction performance of the neural network based setting and the multivariate linear regression model we also implemented a smart version of a neural network, using the same input parameters as the multivariate linear model (CFSPMEQ, CFSPM0° and CFSPM45°) to predict the three vector components CBSPsmartEQ, CBSPsmart0° and CBSPsmart45° as output parameters. For this purpose, we did not perform a preconditioning, but instead used a neural network structure with two hidden layers, 12 neurons per layer and 16 epochs as done for the neural network based on the full set of input parameters. All three reference models were also validated against the corneal back surface measurements CBSPMEQ, CBSPM0° and CBSPM45° in terms of performance with the same metric (squared prediction error) using the test data.

Statistics

The input data, output data and prediction errors are shown with mean, standard deviation, median, minimum – maximum and 90% confidence intervals. The squared prediction error which shows in general a one-sided distribution is described with median and upper 90% and upper 99% confidence intervals. For the multivariate linear prediction model we used maximum likelihood estimation with iterative ECM algorithm,27, 28 and the respective results are described with the regression coefficient matrix and LogL as the value of the log likelihood objective function after the final iteration.28

RESULTS

Table 1 shows the explorative statistics of the seven input parameters used to train the deep learning algorithm (AL, CCT, W2W, patient age, CFSPMEQ, CFSPM0°, CFSPM45°) in terms of mean, standard deviation, median, minimum and maximum and the lower and upper limits of the 90% confidence intervals for the 19,553 measurements.

TABLE 1. Input data for our neural network for prediction of corneal back surface power N = 19,553 Age (years) AL (mm) CCT (mm) W2W (mm) CFSPMEQ (D) CFSPM0° (D) CFSPM45° (D) Mean 71.11 23.68 0.55 11.98 48.82 0.19 −0.02 Standard deviation 10.44 1.34 0.04 0.40 1.70 1.13 0.68 Median 73 23.48 0.55 11.98 48.79 0.14 −0.01 Minimum 25 18.71 0.42 9.71 40.35 −9.36 −9.29 Maximum 99 31.92 0.76 13.92 55.92 8.94 5.70 5% quantile 53 21.90 0.49 11.33 46.07 −1.43 −0.99 95% quantile 85 26.12 0.61 12.64 51.63 2.02 0.94 Note Axial length (AL), central corneal thickness (CCT), and horizontal corneal diameter W2W are directly used from the IOLMaster 700 biometer. The corneal radii and axes of the front surface curvature measured with the IOLMaster 700 are converted to vector components CFSPMEQ, CFSPM0°, and CFSPM45°. Abbreviation: D, dioptres.

After our check for proper conditioning of the neural network with variations of number of hidden layers and number of neurons per layer, we decided to use two hidden layers with 12 neurons per layer for our final version of the prediction algorithm. With this setup (combination of number of hidden layers and number of neurons per layer) the median squared prediction error could be kept below 0.0064 D² (prediction error of 0.08 D, corresponding to the typical repeatability of keratometers or topographers in measuring corneal power). After splitting our dataset into training, validation and test data we checked the performance of the feedforward network. Figure 2 (left side) shows the performance plot of our final version of the feedforward neural network with two hidden layers and 12 neurons per layer, which shows the mean squared prediction error (msPE) of the network applied to the training data, validation data and test data for variations of iterations (epochs). The vertical cyan line indicates proper results with 16 epochs, and the respective msPE for the training data, validation data and test data with 16 epochs are shown in the legend. Figure 2 (right side) displays the performance of the network with the cumulative histogram of the squared prediction error (sPE) for the N = 2933 test data using the feedforward neural network with two hidden layers and 12 neurons per layer optimised with 16 iterations (epochs). The median sPE is provided in the figure title.

Left side: Performance plot of our final version of the feedforward neural network with 2 hidden layers and 12 neurons per layer which shows the mean squared prediction error msPE of the network applied to the training data, validation data and test data for variations of iterations (epochs). The vertical cyan line indicates proper results with 16 epochs, for more iteration cycles no further improvement could be reached. The respective msPE values for 16 epochs of optimisation are shown in the legend for the training data, validation data, and test data. Right side: Cumulative histogram of the squared prediction error sPE for the N = 2933 test data using the feedforward neural network with 2 hidden layers and 12 neurons per layer optimised with 16 iterations (epochs)

Using the multivariate linear model for prediction of corneal back surface power vector components (CBSPmvEQ, CBSPmv0°, CBSPmv45°) from corneal front surface power components measured with the IOLMaster 700 (CFSPMEQ, CFSPM0°, CFSPM45°) we get a transformation matrix with a likelihood objective function value after the last iteration LogL = 5.28e4. From the coefficient matrix structure we see that the dominant elements are located in the diagonal which means that the three vector components of corneal back surface power are mostly determined by the respective component of the corneal front surface vector component, with a scaling factor of around −0.115 (for the vector component in 0°/90°) to −0.137 (for the vector component of equivalent power). In addition we can see that especially for the vector component 0°/90° there is an offset of nearly −0.2 D, whereas for the vector component in the oblique axis (45°/135°) there is only a slight offset of around 0.04 D. Table 2 lists the values of the vector components of the corneal back surface measured with the IOLMaster 700, together with the prediction of the feedforward neural network, the data derived from a fixed ratio of front to back surface radii for both cardinal meridians extracted from the Gullstrand schematic model eye and the predicted vector components for the corneal back surface power from the multivariate linear model as described above (the respective data for the smart version neural network are not shown here). We can see from the table that for the vector component at 0° (with the rule) and 90° (against the rule), the Gullstrand model eye cannot properly represent the corneal back surface power vector components, which leads to the discrepancy between the total corneal astigmatism and keratometric corneal astigmatism based on a corneal front surface curvature measurement only. In contrast, both other prediction strategies (the neural network approach or the multivariate linear regression model) properly resample the mean value and median of all three corneal back surface vector components. TABLE 2. Descriptive statistics of the corneal back surface power in vector components in terms of mean, standard deviation, median, minimum and maximum, 5% and 95% quantile N = 19,553 Measured corneal back surface power (D) Predicted corneal back surface power with neural network (D) Corneal back surface power from Gullstrand model (D) Corneal back surface power from multivariate regression model (D) CBSPMEQ CBSPM0° CBSPM45° CBSPPEQ CBSPP0° CBSPP45° CBSPGEQ CBSPG0° CBSPG45° CBSPmvEQ CBSPmv0° CBSPmv45° Mean −5.83 −0.23 −0.01 −5.83 −0.23 −0.01 −5.88 −0.02 0.00 −5.83 −0.23 −0.01 Standard deviation 0.24 0.16 0.11 0.21 0.12 0.08 0.21 0.14 0.08 0.21 0.12 0.07 Median −5.83 −0.22 −0.01 −5.83 −0.23 −0.01 −5.88 −0.02 0.00 −5.83 −0.23 −0.02 Minimum −7.12 −1.15 0.91 −6.97 −1.16 −0.81 −6.74 −1.08 −0.69 −6.68 −1.13 −0.62 Maximum −5.00 1.28 1.18 −5.18 1.26 0.91 −4.87 1.25 1–12 −4.80 0.84 0.96 5% quantile −6.23 −0.48 −0.18 −6.18 −0.43 −0.12 −6.22 −0.24 −0.11 −6.17 −0.42 −0.12 95% quantile −5.45 0.00 0.15 −5.49 −0.06 0.09 −5.55 0.17 0.12 −5.50 0.06 0.09 Note The first data block (columns 2–4) refers to the corneal back surface power measured with the IOLMaster 700, the next block (columns 5–7) to the corneal back surface power predicted with the feedforward neural network, the next block (columns 8–10) to the back surface power derived from a fixed curvature ratio (front to back surface in both cardinal meridians) extracted from the Gullstrand model eye, and last block (columns 11–13) to the prediction of corneal back surface power using a multivariate linear model with corneal front surface power vector components as predictors (the data for the smart version neural network are not shown in this table). Abbreviation: D, dioptres.

The distributions of the three vector components of the prediction error are displayed in Figure 3 with respect to the vector components derived from the measured corneal back surface curvature CBSPMEQ, CBSPM0° and CBSPM45°. The three vector components of the neural network prediction error (left graph) show a very good performance with a narrow distribution and no systematic offset from the respective components of the measured corneal back surface power. The smart version neural network with the same input parameters as used in the multivariate linear model (second graph from left) also shows a good performance, but the prediction is slightly worse compared to the neural network as shown on the left graph. With the multivariate prediction model (third graph from left) the distributions of the three vector components of the prediction error are much broader, but with no systematic offset. In contrast, with the fixed anterior to posterior curvature ratio extracted from the Gullstrand model eye (graph on the right), the distributions of the prediction error are much broader and show a systematic offset from zero for the vector components equivalent power and astigmatism in 0°/45°.

Normalised histogram together with the kernel probability density distribution plot of the prediction error with respect to the vector components of the measured corneal back surface power equivalent power: CBSPMEQ, astigmatism projected to the 0°/90° meridian: CBSPM0° and astigmatism projected to the 45°/135° meridian: CBSPM45°. The 3 vector components of the neural network prediction error ([left graph] show a very good performance with a narrow distribution and no systematic offset from the respective components of the measured corneal back surface power. The smart version neural network based on the input parameters used in the multivariate linear prediction model [2nd graph from left] shows a slightly lower performance compared to the network based on all input parameters, but performs better compared to the multivariate prediction model. With the multivariate prediction model [3rd graph from left]) the distributions of the 3 vector components of the prediction error are much broader, but without a systematic offset. In contrast, with the fixed anterior to posterior curvature ratio extracted from the Gullstrand model eye (graph on the right) the distributions of the prediction error are much broader and show a systematic offset from zero mainly for the vector components equivalent power and astigmatism in 0°/45°

To check for any trend, Figure 4 displays scatterplots with the prediction error using the neural network approach (in red), the smart version of the neural network (in cyan), the multivariate linear prediction model (in blue) and the prediction using a fixed anterior to posterior corneal curvature ratio extracted from the Gullstrand model eye (in green). In the upper/middle/lower graph the situation for the vector component equivalent power/astigmatism in 0°/90°/astigmatism in 45°/135° are shown, respectively. The prediction error resulting from the neural network shows no systematic offset or effect on the respective component of the measurement. The smart version of the neural network also shows no offset or visible effect on the respective component of the measurement, but a slightly broader scatter compared to the neural network considering all seven input parameters. The prediction of the multivariate linear regression also shows no systematic offset, but an inverse effect on the respective component of the measurement, and the prediction based on the Gullstrand model eye shows both a systematic offset for equivalent power and astigmatism component in 0°/90° as well as an inverse effect in the respective component of the corneal back surface measurement.

Scatterplots comparing the prediction error of four different models with the corresponding parameter as measured using the IOLMaster. In order, from top to bottom, the three plots compare the prediction error against the measured value for the three vector components: equivalent power, astigmatism in the 0°/90° orientation, and astigmatism in the 45°/135° orientation respectively. Each plot is colour coded with the colour of the plotted data points indicating the model as follows: the neural network approach with 7 effect sizes (red), the smart version neural network approach with 3 effect sizes (cyan), the multivariate linear prediction model with 3 effect sizes (blue) and the prediction using a fixed anterior to posterior corneal curvature ratio extracted from the Gullstrand model eye (green). The prediction errors resulting from the two neural network approaches show no (visible) systematic offset or trend compared to the respective component of the measurement. The prediction of the multivariate linear regression shows no vertical offset but does in each case show an inverse trend (negative gradient) against the respective measured vector component, and the prediction based on the Gullstrand model eye shows both a systematic vertical offset mainly for equivalent power and astigmatism component in 0°/90°, in addition to an inverse trend against the respective measured vector component. The scatter of the smart version neural network based model is slightly larger compared to that of the neural network based on all input parameters

DISCUSSION

It is well known that the keratometric measures, especially keratometric astigmatism, do not fully reflect corneal power.1, 10, 15 Keratometers (or ophthalmometers) measure the anterior corneal curvature by projecting keratometer marks onto the cornea reflected off the tear film that are analysed by the examiner. In automatic keratometers or topographers, multiple marks or a Placido pattern with concentric rings are projected onto the cornea, and the reflected mires are evaluated. Based on a schematic model eye – in traditional keratometers, mostly the Gullstrand model eye – the power of the cornea as a meniscus lens is estimated using a keratometer refractive index, which refers to an assumed value characterised by the schematic model eye.2 In ophthalmology, two different keratometer indices are commonly used: the Zeiss index with nK = 1.332, and the Javal index with nK = 1.3375. Based on the Gullstrand model eye, the Zeiss index describes the front vertex power of the model, while the Javal index represents the back vertex power. However in general, such a conversion from corneal front surface curvature to corneal power is not valid, as the corneal back surface is not in a fixed proportion to the front surface, and corneal thickness may vary.2

Javal's rule describes the discrepancy between keratometric astigmatism and the refractive cylinder, and this discrepancy might be due to the effect of corneal back surface curvature or lenticular astigmatism.14, 16 In the last two decades, many attempts have been made to measure the corneal back surface curvature using tomography rather than estimating the effect. Several slit projection systems fulfilling the Scheimpflug measurement condition or otherwise were launched to market that allowed the examiner to look behind the corneal front surface. In the last decade, the development of swept source OCT has made high resolution three-dimensional measurement of the anterior eye segment possible and the corneal back surface can be measured directly.1, 7, 8 Especially in the case of biometry prior to cataract surgery with implantation of toric lenses, which has gained popularity in the last two decades, the effect of corneal back surface astigmatism is highly relevant to achieve optimal refractive results after surgery. Therefore, several prediction models have been proposed,3, 4, 9, 10, 12 either adjusting the keratometric measurement with a fixed offset or based on regression analysis using the vector components of the keratometric astigmatism and in some cases considering the orientation of corneal astigmatism before surgery (i.e., with the rule, oblique or against the rule).

Since 2018, a biometer has become commercially available which allows measurement of corneal thickness and corneal back surface curvature in addition to classical keratometry.17 As the corneal back surface curvature can be measured directly with these biometers, such prediction models for the contribution of corneal back surface astigmatism to the total astigmatism should be based on the measurement instead of using indirect prediction strategies, e.g., using keratometry and refraction data after cataract surgery with implantation of a rotationally symmetric or toric intraocular lens.

Based on a large dataset of biometric data in a population prior to cataract surgery, corneal back surface measurements were used to set up a prediction model which estimates corneal back surface power (in vector components) from the corneal front surface power (in vector components) and several distances within the eye, plus the patient age. A deep learning algorithm with a feedforward neural network23, 24 with three output parameters (predictions of vector components of the corneal back surface power) from the input parameters patient age, AL, CCT, W2W and the three vector components of the corneal front surface power was set up and trained. The potential input parameters ACD and LT were discarded as a stepwise linear regression fit prior to setting up the model did not show a significant improvement with these variables included. A smart version of the neural network was set up using the same input parameters as the multivariate linear prediction model in order to make a direct comparison to the multivariate model. The entire dataset was split into training, validation and test data to make cross-validation and to avoid overfitting of the algorithm. Design parameters such as number of hidden layers or number of neurons per layer were evaluated before the final design of the algorithm was defined. The mean squared prediction error as the Euclidian norm of the prediction errors of the three vector c