**Purpose**:
Stochastic eye models are a method to generate random biometry data with the variability found in the general population for use in optical calculations. This work improves the accuracy of a previous model by including the higher-order shape parameters of the cornea.

**Methods**:
The right eye biometry of 312 subjects (40.8 ± 11.0 years of age) were measured with an autorefractometer, a Scheimpflug camera, an optical biometer, and a ray tracing aberrometer. The corneal shape parameters, exported as Zernike coefficients, were converted to eigenvectors for dimensional reduction. The remaining 18 parameters were modeled as a sum of two multivariate Gaussians, from which an unlimited number of synthetic data sets (SyntEyes) were generated. After conversion back to Zernike coefficients, the data were introduced into ray tracing software.

**Results**:
The mean values of nearly all SyntEyes parameters were statistically equal to those of the original data (two one-sided *t*-test, *P* > 0.05/109, Bonferroni correction). The variability of the SyntEyes parameters was similar to the original data for most important shape parameters and intraocular distances (*F*-test, *P* < 0.05/109), but significantly lower for the higher-order shape parameters (*F*-test, *P* > 0.05/109). The same was seen for the correlations between higher-order shape parameters. After applying simulated cataract or refractive surgery to the SyntEyes model, a very close resemblance to previously published clinical outcome data was seen.

**Conclusions**:
The SyntEyes model produces synthetic biometry that closely resembles clinically measured data, including the normal biological variations in the general population.

^{1}Navarro and colleagues,

^{2,3}and Liou and Brennan,

^{4}that have played an essential role in many great scientific developments of the last century. However, despite these contributions, such generic eye models do not take the wide variety in ocular biometry into account that was observed in epidemiologic biometry studies. For example, assuming a tolerance of

*±*0.25 diopters (D) in keratometry and ±0.10 mm in axial length, the Gullstrand eye model would correspond with only 0.51% of the 1178 Western European eyes in our previous work on the epidemiology of ocular biometry.

^{5}Similarly, the more recent Navarro and Liou and Brennan models both correspond with 1.02% of the same cohort. To study eyes that deviate from these average values, these generic models are often customized, for example, by inserting measured biometry values into a chosen eye model while keeping other unmeasured (or unmeasurable) values unaltered.

^{6–8}Alternatively, customization can also be achieved by considering the unmeasured parameters as free variables in an optimization process until, for example, the ocular wavefront of the individual eye model matches the measured values.

^{9,10}Although this approach leads to a data set much closer to the patient's actual biometry, one has to make sure that the inserted biometry parameters are very complete to avoid introducing systematic errors resulting from a lack of correlation between the measured data and the unaltered parameters of the original eye model. Moreover, access to measured data is still required for customized modeling.

^{11}based on ideas proposed earlier by Sorsby et al.,

^{12}Thibos et al.,

^{13}and Zhao.

^{14}However, although this model successfully reproduced the distributions of the ocular biometry in West European subjects, it was restricted to producing only lower orders of wavefront aberrations, thus limiting its applications.

^{15}more accurate methods to estimate the crystalline lens shape

^{16}and power,

^{17}and a linear combination of multivariate Gaussians

^{5}instead of refractive filtering. After verification of the model, examples will be given of possible applications.

*T*and axial length

*L*were measured with a noncontact optical biometer (Haag-Streit Lenstar, Koenitz, Switserland), and the ocular wavefront was determined using a ray tracing aberrometer (Tracey iTrace, Houston, TX, USA). Furthermore, the refraction was determined using an autorefractometer (Nidek ARK 700, Nidek, Gamagori, Japan) and the intraocular pressure with a pneumotonometer (Reichert ORA, Depew, NY, USA). A detailed overview of the parameters used in the modeling is given in Table 1.

**Table 1**

^{2}from which the basic structure and some of the clinically inaccessible parameters are taken. It comprises two Zernike surfaces representing the cornea, two aspherical surfaces representing the crystalline lens, a stop at the pupil plane, and a spherical retina (Table 2). The anterior lens surface also contains an additional Zernike phase correction to account for wavefront aberrations of noncorneal origins, such as the intrinsic aberrations of the crystalline lens and those resulting from the lenticular alignment relative to the cornea.

**Table 2**

^{15}The remaining 0.5% of variability lies mostly in the higher-order aberrations and is of no real consequence to the model (Supplementary Material A). These eigencorneas are used during the stochastic process and converted back to Zernike coefficients and CCT afterward.

*P*

_{L}is calculated from the ocular biometry using the Bennett equation with optimized parameters,

^{17,18}and the lens radii

*r*

_{La}and

*r*

_{lp}are derived from regressions based on lens power

*P*

_{L}and lens thickness

*T*.

^{16}Similarly, a regression was used to estimate missing lens thickness values in eyes where it is missing.

^{16}Once lens power and radii are available, the equivalent refractive index of the lens may be derived using equation 3 in Reference 16. Finally, the asphericity of both lens surfaces is taken from the Navarro model.

^{2}

*Z*

_{n}

^{m}[

*int*]; diameter, 5.50 mm) to account for the aberrations of noncorneal origins, which is equivalent to including a Zernike phase plate at the iris as both the anterior lens apex and the pupil lie in the same plane. These constant coefficients were determined by matching the wavefront calculated using the measured biometry of the original subjects with the measured wavefront.

*LT*

_{x}

*, LT*

_{y}) and shift (

*LS*

_{x}

*, LS*

_{y}) could not be determined in the eye of the original subjects, these parameters are not included in the multivariate Gaussian model. Instead these are generated independently based on mean and standard deviation (SD) values found in the literature.

^{19,20}

^{21}to add a fixed value for retinal thickness

*RT*of 0.200 mm to the measured axial length

*L*prior to any calculations. Vitreous depth

*V*is then calculated as

*V*=

*L*−

*T*−

*ACD*−

*CCT*−

*RT*.

*ACD*, axial length

*L*, lens power

*P*

_{L}, anterior lens radius

*r*

_{La}, posterior lens radius

*r*

_{Lp}, and lens thickness

*T*. These parameters form an 18-dimensional point cloud that may be fitted with a linear combination of multivariate Gaussian functions,

^{5}which forms the core of the stochastic model. The fit was performed in Matlab 2011b (The Mathworks, Natick, MA, USA) using an expectation-maximization (EM) algorithm using the procedures described in References 5 and 11. Early trials indicated that a fit of two multivariate Gaussians gave a good fit accuracy, which could not be improved significantly by increasing the number of Gaussians.

^{22}During the development of both the ray tracing module and the statistical model, the ray tracing results were verified against results obtained with commercial software (Zemax, Zemax LLC, Kirkland, WA, USA) to ensure compatibility.

*t*-test (TOST) is used, which defines certain thresholds of equivalence between which the average of both sets may be considered equal.

^{23,24}Although some of the parameters involved may not be normally distributed, one can still use parametric testing such as

*t*-tests, provided the populations involved are sufficiently large. Lumley et al.

^{25}stated that in most cases this threshold lies below 100 cases for a mild degree on nonnormality and below 500 for extremely high nonnormality. Given that the worst distributions in this work are only moderately nonnormal (as verified with Kolmogorov-Smirnov tests), this approach is warranted by the standards set by Lumley et al. All tests are performed at a confidence level of 0.05, adjusted with a Bonferroni correction in case of multiple simultaneous comparisons (indicated by 0.05/N, with N the number of comparisons performed).

*J*

_{45}(Fig. 1c), the histogram shows a very good match with those of the autorefractometer and wavefront measurements, whereas for

*J*

_{0}(Fig. 1b), the distributions were a bit flatter, with a slightly higher occurrence of with-the-rule astigmatism in the calculated wavefronts compared with the measured wavefronts. The average Zernike coefficients presented in Figure 1d are all significantly equal to each other (TOST,

*P*> 0.05/45, Bonferroni correction for 45 repeated comparisons), but in many of their SDs, significant differences are seen (

*F*-test,

*P*< 0.05/45; Fig. 1e). For the wavefront calculated based on the full corneal elevation, the SD is often significantly larger than that of the measured wavefront, whereas for the eigencornea compressed data, the SD is often significantly smaller (

*F*-test,

*P*< 0.05/45). Only for the second-order aberrations were the SDs similar to those of the measured wavefronts. Although this variability could be improved by including more than 12 eigencorneas, such a measure would also increase the dimensionality of the model, thus increasing the risk of overfitting the original data. We therefore opted to accept this as a limitation of the proposed model. More information on how well different numbers of eigencorneas reproduce the variability of the original wavefront can be found in Supplementary Material A.

**Figure 1**

**Figure 1**

*P*> 0.05/109; Table 3). Using the statistical eye model to generate 1000 SyntEyes and comparing their mean parameter values to those of the original data, both are significantly equal (TOST,

*P*> 0.05/109). For the SD of the 18 parameters used by the eye model, no significant differences are seen between SyntEyes and the original data (

*F*-test

*P*> 0.05/109), but for most corneal Zernike coefficients

*Z*

_{n}

^{m}(

*CA*) and

*Z*

_{n}

^{m}(

*CP*), the SDs are significantly lower than those of the original data (

*P*< 0.05/109). This is most likely the result of the eigencornea compression. A complete version of Table 3 is available in Supplementary Material D.

**Table 3**

**Figure 2**

**Figure 2**

*P*< 0.05/48; Table 4). Often this is due to the presence of several outliers, but for

*SE*,

*Z*

_{0}

^{0}(

*WF*), and

*Z*

_{2}

^{0}(

*WF*), the distributions were obviously leptokurtic and skewed (Fig. 3). Based on the criteria by Lumley et al.,

^{25}our population size is sufficiently large for parametric statistics, however, so the TOST and

*F*-test could be performed for these parameters.

**Table 4**

**Figure 3**

**Figure 3**

*P*> 0.05/48; Table 4), but the SDs of the SyntEyes is in most cases significantly smaller than that of the original data (

*F*-test,

*P*< 0.05/48). This is also illustrated in Figure 3, where the histograms of the calculated Zernike coefficients up to the fourth order are shown for 1000 SyntEyes (black lines) and compared with the measured wavefronts of the original data (gray bars). Overall, the agreement between both is good, except for the trefoil coefficients (

*Z*

_{3}

^{±3}), where the SyntEyes display considerably less variation compared with the original data. To a lesser degree, this is also seen for

*Z*

_{4}

^{−4}and

*Z*

_{4}

^{+2}. A complete version of Table 4, as well as a comparison of the wavefront metrics proposed by Marsack et al.,

^{26}is available in Supplementary Material D.

_{0}

^{0}), tilt (

*Z*

_{1}

^{±1}), defocus (

*Z*

_{2}

^{0}), astigmatism (

*Z*

_{2}

^{±2}), coma (

*Z*

_{3}

^{±1}), and spherical aberration (

*Z*

_{4}

^{0}). For the other third- and fourth-order Zernike coefficients, the correlations between the SyntEyes parameters were higher than between the parameters of the original data, associated with the SyntEyes' lower degree of randomness (Fig. 4).

**Figure 4**

**Figure 4**

^{27}As this IOL has a 0% rate for posterior capsular opacification

^{27}and a very high degree of postoperative centration and rotation stability,

^{28,29}it is an ideal lens to verify the optical performance of the SyntEyes with clinical data. All lens powers were calculated using the Sanders-Retzlaff-Kraff/Theoretical (SRK/T) formula,

^{30}aiming at postoperative emmetropia. Information on the radii of curvature and lens thickness were obtained directly from the manufacturer (Morcher GmbH, Stuttgart, Germany).

**Figure 5**

**Figure 5**

*Z*

_{4}

^{0}(

*WF*) of the pseudophakic SyntEyes lies close to 0.20 μm, which is in the range of the values reported in the literature for spherical IOLs.

^{31}(37 patients; male/female: 29.7%/70.3%; age: 33.6 ± 9.4 years; postoperative date: 6.2 ± 4.5 months; no ocular comorbidities; no cylinder correction >1 D). These eyes were treated with an InPro GAUSS laser system (InPro GmbH, Norderstedt, Germany), which delivers a broad Gaussian laser beam with a diameter of 6.5 mm. This process can be simulated by removing a similar profile from the anterior cornea for all 365 SyntEyes with a myopia higher than −1 D (taken from the same set of SyntEyes as before). The amount of ablation required for emmetropia was determined using the Munnerlyn algorithm for myopia.

^{32}Assuming perfect centration and a perfect beam profile, the postrefractive

*SE*of the SyntEyes is similar to what is found in the real post-LASEK eyes (Fig. 6a), albeit with a sharper distribution. The Zernike coefficients of both the pre- and postoperative wavefronts are similar for both the real eyes and the SyntEyes (Figs. 6b, 6c), apart from the preoperative astigmatism (

*Z*

_{2}

^{+2}), which is higher in the real eyes and two postoperative third-order coefficients (

*Z*

_{3}

^{−3},

*Z*

_{3}

^{−1}) that are more pronounced in the real eyes as well. As before, the variability of the Zernike coefficients is higher in the real eyes that in the SyntEyes (not shown). The SyntEyes data for this example are included in Supplementary Material C.

**Figure 6**

**Figure 6**

^{33}Similarly, the effects of laser refractive surgery can be simulated for a wide variety of biometry values, which could potentially improve the postoperative outcome. Alternatively, one could simulate how the IOL power should be calculated in postrefractive patients to identify combinations of biometric parameters for which a surgical procedure could lead to previously unexpected outcomes. Finally, the SyntEyes model may be useful for governments or regulatory authorities (e.g., Food and Drug Administration and European Medicines Agency) to refine the indications for reimbursement of certain therapies or improve laser safety thresholds.

^{34,35}

^{36–39}wide angle refraction,

^{2,3,40–44}a gradient index lens,

^{4,36,44,45}accommodation,

^{1,2,36,46,47}or corneal biomechanics.

^{48}

^{16}Given the stochastic nature of the SyntEyes model, which only requires the mean and covariance values in order to produce plausible estimates for the lens curvatures values, these parameters are probably well integrated into the model. Nevertheless, once phakometry equipment becomes available clinically, for example, based on the crystalline lens topography,

^{49}we will consider updating the current model with these parameters to keep the current model relevant.

**J.J. Rozema**, None;

**P. Rodrigue**z, None;

**R. Navarro**, None;

**M.-J. Tassignon**, None

*. Hamburg, Germany: Voss ; 1909; 301–358, 382–415.*

*Handbuch der Physiologischen Optik*. 3rd ed*. 1985; 2: 1273 –1280.*

*J Opt Soc Am A**. 1999; 16: 1881 –1891.*

*J Opt Soc Am A**. 1997; 14: 1684 –1695.*

*J Opt Soc Am A**. 2014; 91: 713 –722.*

*Optom Vis Sci**. 2005; 116: 80 –85.*

*Optik**. 2006; 53: 1605 –1618.*

*J Modern Opt**. 2007; 15: 2204 –2218.*

*Opt Express**. 2006; 83: 371 –381.*

*Optom Vis Sci**. 2011; 2: 1649 –1662.*

*Biomed Opt Express**. 2011; 52: 4525 –4533.*

*Invest Ophthalmol Vis Sci**. 1981; 65: 805 –811.*

*Br J Ophthalmol**. 2002; 19: 2329 –2348.*

*JOSA A**. 2009; 34: 7–9.*

*Optics Lett**. 2014; 34: 667 –677.*

*Ophthalmic Physiol Opt**. 2012; 53: 2533 –2540.*

*Invest Ophthalmol Vis Sci**. 2011; 52: 7937 –7942.*

*Invest Ophthalmol Vis Sci**. 2010; 10 (14): 34.*

*J Vis**. 2007; 47: 71–84.*

*Vis Res**. 1988; 8: 53 –59.*

*Ophthal Physiol Optics**. 2015; 56: 3577 –3583.*

*Invest Ophthalmol Vis Sci**. 2002; 19: 657 –669.*

*JOSA A**. 1976; 9: 741 –744.*

*Biometrics**. 1987; 15: 657 –680.*

*J Pharmacokinet Biop**. 2002; 23: 151 –169.*

*Ann Rev Publ Health**. 2004; 4 (4): 328.*

*J Vis**. 2011; 37: 2120 –2129.*

*J Cataract Refract Surg**. 2007; 33: 1267 –1272.*

*J Cataract Refract Surg**. 2009; 35: 1385 –1388.*

*J Cataract Refract Surg**. 1990; 16: 333 –340. (Erratum:*

*J Cataract Refract Surg**J Cataract Refract Surg*. 1990;16:528).

*. 2010; 51: 2800 –2804.*

*Invest Ophthalmol Vis Sci**. 1988; 14: 46–52.*

*J Cataract Refract Surg**. 2015; 160: 403–405.*

*Am J Ophthalmol**. Wright-Patterson AFB, OH: 2012.*

*Probabilistic Model for Laser Damage to the Human Retina**. 2012: 822114 –822116.*

*SPIE BiOS**. 2014; 14 (13): 21.*

*J Vis**. 2008; 8 (4): 29.*

*J Vis**. 1992; 9: 2111 –2117.*

*JOSA A**. 1989; 29: 1685 –1692.*

*Vision Res**. 1983; 73: 1544 –1550.*

*JOSA A**. 1984; 61: 166 –176.*

*Am J Optometry Physiol Optics**. 1974; 58: 709 –714.*

*Br J Ophthalmol**. 2006; 46: 2236 –2250.*

*Vision Res**. 2007; 24: 2157 –2174.*

*JOSA A**. 2004; 29: 1197 –1199.*

*Optics Lett**. 2002; 22: 201 –208.*

*Ophthalmic Physiol Optics**. 1999; 76: 720 –727.*

*Optometry Vision Sci**. 2011; 133: 011002.*

*J Biomechan Engin**. 2014; 5: 3547 –3561.*

*Biomed Opt Express*