Laser-induced corneal injury: validation of a computer model to predict thresholds

: The exposure and emission limits of ICNIRP, IEC 60825-1 and ANSI Z136.1 to protect the cornea are based on a limited number of in-vivo studies. To broaden the database, a computer model was developed to predict injury thresholds in the wavelength range from 1050 nm to 10.6 µm and was validated by comparison with all applicable experimental threshold data (ED 50 ) with exposure duration between 1.7 ns and 100 s. The model predictions compare favorably with the in-vivo data with an average ratio of computer prediction to ED 50 of 0.94 (standard deviation ± 30%) and a maximum deviation of less than 2. This computer model can be used to improve exposure limits or for a quantitative risk analysis of a given exposure of the cornea.


Introduction
Laser radiation in the infrared range with wavelengths above 1 µm is absorbed in the cornea to a sufficient degree so that thermally induced corneal injuries become relevant. Exposure limits to protect the eye are promulgated on the international level by ICNIRP [1] and in the United States of America in ANSI Z136.1 [2], where they are referred to as the maximum permissible exposure (MPE). Cornea-related MPEs are historically based on injury thresholds obtained with rabbit and non-human primate models. The injury thresholds, usually expressed as corneal radiant exposure or corneal irradiance, show a strong dependence on wavelength and pulse duration and, to a degree, on the diameter of the laser beam incident on the cornea. While the collection of experimental data covers the parameter range sufficiently to set MPEs, computer models can form the basis for a systematic comparison of injury thresholds with MPEs for the full range of wavelengths, pulse durations and beam diameters of interest. Particularly interesting for such a comparison is the interdependence of trends with wavelength and pulse duration. To our knowledge, studies comparing experimental corneal injury thresholds with predictions of computer models based on the Arrhenius integral have been limited to single wavelengths [3][4][5]. A complete validation of a computer model for the whole applicable wavelength and pulse duration range has not been available.
The current paper describes in detail a computer model and demonstrates its validity on the basis of a comparison with experimental data in the relevant range of pulse durations longer than 1 ns and wavelengths above 1 µm so that a thermal bulk model is applicable. An earlier version of the model was presented in 2011 [6]; since then the model has been improved and the validation extended into the nanosecond pulse duration regime. Recent experimental data and multiple pulse data is included for the validation of the model. A systematic comparison of multiple pulse thresholds against the respective MPE values is discussed in another publication [7]. For completeness it is mentioned that in the wavelength range up to approximately 1400 nm, the retina is also at risk. The potential hazard to the retina is an issue that is biophysically not related to corneal injury (other than that absorption in the pre-retinal ocular media results in a higher retinal injury threshold when expressed as power or energy incident on the surface of the eye) and does not need to be considered in the present work. Correspondingly, safety standards such as ANSI Z136.1 in the wavelength range between 1300 nm and 1400 nm also define separate and independent MPEs to protect the retina and to protect the cornea, respectively.

Anatomy and physiology
The cornea is the outermost transparent part of the eye, with a diameter of about 11 mm in humans. The space between the cornea and the lens of the eye is filled with the aqueous humor. The cornea functions as the primary refractive element in the human eye. The outermost layer of the cornea is the epithelium, which consists of about six cell layers resulting in a thickness of about 50 µm [8,9]. New cells are formed in the basal cell layer of the epithelium and these cells gradually move outwards, where they are shed into the tear film. One cycle from cell formation to cell shedding takes about 10 days in the normal healthy cornea. Upon injury, cell turnover becomes faster and the epithelial cell layer heals completely within a few days. The bulk of the cornea is made up by the stroma, separated from the epithelium by the Bowman's layer. The innermost part of the cornea consists of the Descemet's membrane and the endothelium, all 5 µm to 10 µm thick [10]. The cornea of non-human primates (NHP) from the genus macaque is approximately 500 µm thick [11], while that of rabbits is thinner at about 300 µm to 400 µm [12,13]. The corneal epithelium is covered by a tear film. Its thickness, including lipid and aqueous layers, depends on many physiological and environmental factors but for humans is reported to be between 3 µm and 8 µm [9,14,15].

Nature of the injury
This work is concerned with thermally induced corneal injury; photochemically induced corneal injury associated with ultraviolet radiation is not considered. Absorbed laser energy translates into heat and an injury can occur when critical temperatures are exceeded in the tissue. At the cellular level, thermal injury to corneal keratocytes and stromal collagen fibrils disturbs the periodic structure of the cornea that is the basis for its transparency. Consequently, a thermally induced injury of the cornea can be detected by ophthalmoscopic means as clouding. Superficial lesions that are confined to the epithelium layer show partial or total repair within a few days [16,17] while deep lesions affecting the stroma remain permanent [18,19]. Unhealed injuries leave scars, which lead to scattering and distorted wavefronts, thereby reducing visual acuity [20].

Endpoint for threshold studies
Studies conducted to support the setting of safe exposure limits need to be based on the determination of the exposure level that just results in a minimum visible lesion (MVL), referred to as threshold lesion. The common endpoint for such studies is the detection of a MVL with a slit lamp, appearing as greyish white spot as shown in Fig. 1. Lesions as small as 100 µm in diameter can be detected by ophthalmic means [18]. Near threshold lesions may take some time to become macroscopically visible, thus a time delay must be introduced between the exposure and the observation and determination of the lesion. A threshold lesion typically appears within 1 hour and lesions involving the stroma show no worsening at 24 h after exposure [18,19]. In the case of a threshold lesion, healing processes can be observed starting at 6 hours post exposure [21]. Exposure levels for which an injury is visible immediately after exposure constitute a "super-threshold" exposure when compared to thresholds determined by examination at some time after exposure. Figure 1 shows a corneal lesion obtained with a CO 2 laser for a radiant exposure of 6.8 J·cm −2 significantly above the injury threshold of 2.3 J·cm −2 and the MPE of 0.3 J·cm −2 . The exposure duration was 100 ms and the 1/e diameter of the Gaussian beam was 2.5 mm [10]. A threshold is usually defined as the dose having a 50% probability of resulting in a visible lesion and is referred to as the ED 50 . Discussions on the meaning of the ED 50 level and its application to corneal injuries can be found in Ref. [22,23], respectively. The ED 50 can be obtained by probit analysis of the yes/no injury data [24]. Alternatively, an injury threshold level can be obtained by bracketing responsive and non-responsive exposure levels [25]. There is no evidence that the results obtained with these two methods differ significantly when applied to corneal threshold injuries. However, probit analysis offers some advantages, such as the possibility to quantify the steepness of the tissue response expressed by the ratio of ED 84 to ED 50 , commonly referred to as slope of the dose-response curve. For the identified experimental data that were obtained by probit analysis, the average slope was 1.16 (36 samples, maximum slope 1.47), which corresponds to a standard deviation of ±16%, thus indicating a rather steep response as compared to other bioeffects or tissues. Noticeably, the average ratio of lowest exposure level leading to an observable injury to the ED 50 level was 0.85 (20 samples, minimum ratio 0.63). In the remainder of the paper, the symbol ED 50 is used even for the case that the threshold was determined by bracketing and not by probit analysis.

Optical properties
For wavelengths up to 2500 nm, absorption data are available for the NHP cornea, aqueous and lens in a comprehensive review of the International Commission of Illumination [26]. Between 2500 nm and 20 µm, absorption values for the anterior parts of the eye can be assumed to be equivalent to that of water, as shown in Fig. 2. The penetration depth at which the local irradiance drops to 1/e of the irradiance incident on the surface (see section 3.1) exhibits a particularly strong dependence on wavelength, varying from approximately 1 cm at a wavelength of 1000 nm to less than 1 µm (i.e. well within the tear film) at wavelengths around 3 µm. For wavelengths of 2.5 µm and longer, the incident radiation is completely absorbed within the cornea.

Experimental thresholds
A literature review was conducted in order to identify all relevant experimental studies that report laser induced injury thresholds in-vivo for the cornea. A total of 178 experimental injury thresholds were identified; 28 with the non-human primate model (NHP), 150 with the rabbit model. The range of relevant experimental data reached from 1064 nm to 10.6 µm in terms of wavelength and from 1.4 ns to 1800 seconds in terms of exposure duration. In order to use a consistent term when referring to a wide temporal regime, such as in plots, "exposure duration" is preferred to "pulse duration" even when that regime includes nanosecond pulses. Ultrashort-pulse data in the sub-ns time domain was excluded from the validation due to non-linear absorption effects as well as potentially non-bulk thermal damage mechanisms. No experimental data is available in the range between 250 ns and 100 µs. The irradiance profile on the cornea was either Gaussian or "top hat". For the specification of the "diameter" of a Gaussian beam profile, the 1/e irradiance level is consistently used throughout this paper, so that dividing the total power by the area defined by the 1/e diameter results in the peak irradiance of the Gaussian beam profile. The beam diameter incident on the cornea ranged from 100 µm to 10 mm. Multiple pulse data with up to 1024 pulses were also identified and included.
The relevant injury thresholds [3,13,[16][17][18][19][20][21]25, are listed in Table S1 along with the predictions of the computer model in the form of a ratio referred to as R ED50 , defined as the ratio of the computer model injury threshold to the experimental ED 50 . The predicted lesion depth is also given (zero being set as the exterior surface of the tear film; more details in section 3.3). From a total number of 178 thresholds, 9 were not included in the comparison with the predictions of the computer model, for the following reasons.
For the wavelength of 10.6 µm and a pulse duration of 1.4 ns, the threshold of 14 mJ·cm −2 reported by Mueller and Ham [52] was not included in the validation of the computer model since it is more than a factor of 10 lower than four other thresholds in the same pulse duration regime [46,48]. The inconsistency was discussed in Ref. [55,23]. The injury threshold obtained with a HF laser (2.9 µm) in that same study was also not included because the beam profile was not sufficiently defined to permit modelling. According to the authors, the beam profile was "poorly Gaussian with hot and cold spots".
For the wavelength of 10.6 µm and a pulse duration of 70 ms, the thresholds of 0.48 J·cm −2 and 0.41 J·cm −2 for the NHP and rabbit models reported by Borland et al. [30] were a factor of about 3 lower than comparable exposures obtained with a pulse duration of 55 ms [32] and 31 ms [13] (these thresholds can be scaled, based on typical pulse duration dependence, to compare against the Borland data more accurately). This discrepancy was discussed at a bio-effects research meeting [56] under the Technical Cooperation Program TTCP in Ottawa, Canada in 1972, with one of the authors of this paper (BES) present. The thresholds reported by Borland et al. were found to be affected by diffraction due to an aperture placed about 12 cm to 15 cm in front of the eye, leading to an annulus beam profile at the cornea instead of a top-hat. Compared to a top-hat and the reported average irradiance value, the actual irradiance was calculated to be a factor of 2 to 2.5 higher. Some of the lesions shown at the meeting also were annular in shape. Therefore the Borland data are not included in the validation of the computer model.
The injury threshold reported by Fine et al. [34] for a 30-minute exposure was not included in the model validation based on the rationale that the cornea was intentionally not irrigated but left to dry during the exposure, resulting in an injury threshold that is not representative for normal exposure conditions.
Two threshold data from a study by McCally et al. [37] with pulse durations of 36 ms and 56 ms were not included in the validation on the grounds that they are consistently lower than threshold data obtained with similar exposure parameters but shorter pulse durations (see Fig. 3(a)), which is inconsistent with biophysical principles. . The injury thresholds were computed for a 1 mm beam diameter, and for exposure durations longer than 100 ms, experimental data points were selected to have approximate beam diameters of 1 mm in order to avoid deviations resulting from a beam-diameter dependence.
One datum reported by Zuclich et al. [46] for a pulse duration of 25 ns (wavelength 10.6 µm) was not included in the model validation because it is significantly higher than other injury thresholds with comparable laser parameters (see Fig. 3(b)), including from the same study; the absence of a pulse duration dependence is also dictated by biophysical principles in the regime of thermal confinement, which is confirmed by the other available data.
Finally, one threshold reported [47] for a wavelength of 1410 nm by Archibald and Taboada was not included for reasons discussed in section 5.1.

Optical properties
The attenuation of irradiance E as a function of depth z due to absorption was modelled according to the Beer-Lambert's law: where R is a reflection coefficient and the depth z equals zero at the anterior surface of the absorbing medium (here, the tear film). The wavelength dependent absorption coefficient α is discussed in subsection 2.4. Specular reflection at the anterior surface of the eye is taken into account by the Fresnel equations using the refraction index for water and air as function of wavelength [27]. Reflections from other surfaces (such as from the back of cornea or from the lens) are neglected. The model assumes the laser is incident on a flat absorbing surface ("the cornea"). It is considered as justified that the curvature of the cornea does not play a significant role in the thermal problem at hand. Simple geometrical optics confirms that the optical power of the cornea does not significantly alter the beam size within the cornea. By applying Snell's law of refraction, it can be shown that the largest angle of refraction within the cornea does not exceed 11°for a collimated beam of 8 mm in diameter and the larger the beam, the less important the local variations at the edges become since the highest temperatures will be obtained in the center of the corneal irradiance profile. Therefore it can be assumed that the beam diameter remains constant throughout the corneal tissue.

Transient increase in temperature
A time-dependent thermal model was developed using a finite-element Multiphysics package (COMSOL 3.5a, Comsol AB, Stockholm, Sweden, 2008), with two contiguous flat slabs representing the cornea and aqueous chamber. Both tissues exhibit homogeneous and isotropic properties. In this model, the tear film was not differentiated from the other tissues. Axial symmetry of the laser beam and the cornea reduces the problem to a 2D case. Besides heat conduction, convection was considered in the aqueous chamber via a negative heat source (heat sink) representing the flow of the aqueous humor. The extent of the finite-element domains, both radially and in depth, was set to allow sufficient heat diffusion, i.e. varying with exposure duration and penetration depth. The depth of the aqueous chamber was considered sufficient so that it was not necessary to model the deeper tissues such as lens or iris. The boundary conditions were adiabatic, except at the air-cornea interface where heat losses were accounted for by non-linear boundary conditions. The equation to be solved can be written in the form of the Penne's bio-heat equation: where T(t,r,z) is the increase in temperature, Q(t,r,z) represents the heat source (laser radiation) and S aq the heat sink in the aqueous humor. The parameters k, ρ and C are conductivity, density and heat capacity, respectively. The time, radial and axial variables are noted t, r and z, respectively. Heat losses occurring in the aqueous chamber via humor flow can be inserted in the bio-heat equation as a perfusion term for this domain, ω being the perfusion rate and T aq the average aqueous temperature. The boundary conditions were as follows: where the heat flux at the air-cornea boundary depends on a differential between surface temperature T s and ambient temperature T a (both in Kelvin). The first term represents the contribution of radiative heat loss (σ and ε being the Stefan-Boltzmann constant and the emissivity of the cornea, respectively), the second term the contribution of free convection with h as heat convection coefficient. A third arbitrary term containing a scaling factor e and an exponent n to be optimized was introduced. The values used for all parameters discussed above are summarized in Table 1, including references, if applicable. Values annotated as "optimized" were adjusted so as to minimize the spread of R ED50 ratios (ratio of computer model injury thresholds to experimental ED 50 ).

Damage model
A thermally-induced threshold injury can be conceptualized as an accumulation of microscopic sub-lethal damage [59], which eventually leads to cell death by apoptosis or necrosis [60]. Such pathway can be modelled by using the Arrhenius equation which describes the temperaturedependent rate of reaction κ and when integrated over time results in a measure of macroscopic damage: where T(t) is the solution of the heat equation (unit: K) and R the ideal gas constant. The parameter E is related to an activation energy of the reaction (energy barrier to overcome so the reaction takes place) and the parameter A represents the rate of reaction. The output of the computer model is the lowest exposure level that leads to Ω = 1 within the cornea and for a given lesion diameter; this predicted injury threshold can be directly compared to experimental ED 50 values. The reference point for calculating the damage integral was set to 100 µm away from the beam axis, which is equivalent to a minimum visible lesion diameter of 200 µm. The depth at which the lesion first occurs within the cornea is not predefined, but is automatically detected by an algorithm. The tear film was excluded from injury determination. In order to improve the quality of the computer model predictions for highly absorbing wavelengths, the tear film thickness δ was varied between zero and 15 µm and an optimum thickness of 9 µm was found. The results reported in Table S1 show that the threshold lesion can occur at various depths depending on the exposure parameters. All parameters used in the damage model are listed in Table 2. For a valuable review of the concept of Arrhenius integral and its application to thermal insult see Ref. [61]. Values annotated as "adapted" were adjusted so as to reduce the spread of R ED50 s.

Model validation
The computer model described in section 3 was used to reproduce injury thresholds for those wavelengths, exposure durations and corneal laser beam diameters for which relevant experimental threshold data is available, according to section 2 (169 data). The main figure of merit used to evaluate the validity of the computer model was the ratio R ED50 , defined as the ratio of computer model injury threshold to experimental ED 50 (see Table S1). Thus for R ED50 > 1, the predictions of the computer model can be interpreted as potentially to err on the unsafe side. Further figures of merit such as the trend of R ED50 with different laser parameters, animal models or assessment delays were also investigated to confirm the general validity of the computer model. Figure 4 shows the distribution of ratios R ED50 with a mean value of 0.94 and the standard deviation of 30%, as listed in Table 3, together with other relevant descriptive statistics.  The collection of ratios R ED50 can also be illustrated as function of exposure duration, wavelength and beam diameter (Fig. 5). No trend indicative of a systematic deviation of the  5. Distribution of individual R ED50 ratios as a function of (a) exposure duration, (b) wavelength and (c) beam diameter; to provide further information, the data in each diagram is split in (a) wavelengths shorter or longer than 3 µm, (b) exposure durations shorter or longer than 100 µs and (c) assessment delays shorter or longer than 1 h after exposure. model predictions with a laser parameter was found. It is worth mentioning that the computer model was set up with parameters from the literature for the geometry, thermal and optical parameters, leaving only the heat losses and damage model to be adjusted in order to reduce the deviations of model results against experimental data. On the topic of corneal thickness, it was not necessary to differentiate between rabbit and NHP models. The computer model predictions were best in line with experimental ED 50 s for a corneal thickness of 400 µm regardless of the animal model. For penetrating wavelengths (and hence a shallow temperature gradient over depth), a thicker cornea -such as 500 µm -would result in a slightly higher threshold levels as the lesion depth increases. The difference in threshold values is not very significant, less than 10% for the exposures considered in the validation and approaching 0 for long exposures). It can be seen that the ratios R ED50 are approximately evenly distributed around a value of 1 regardless of the range of the laser parameter under consideration. Furthermore, for all relevant "binary" variables such as animal model (NHP vs. rabbit), assessment delay (less than 1 h or more than 1h and up to 24 h), beam profile (uniform or Gaussian), injury threshold determination technique (probit analysis or bracketing) or reported lesion depth (epithelium or stroma), the difference between R ED50 ratios were not statistically significant for the two groups (p-value consistently larger than 0.05 on both geometric mean and standard deviation). The validation also included repetitive pulses with up to 1024 pulses at different wavelengths and pulses durations including ns-duration pulses.

Overview of deviations
Over the entire range of laser parameters for which thermally-induced threshold injuries of the cornea are available in animal models, the largest deviation of our computer model predictions to the "unsafe" side, i.e. where model predictions are higher than experimental results, was R ED50 = 1.75. This deviation was ranked as relatively small considering the number of samples as well as the variation in absorption coefficient (four orders of magnitude) together with the eleven orders of magnitude in exposure duration. Also the deviation to the other side, i.e. where the model threshold predictions are lower than the experimental data, is not significantly larger with a maximum deviation of R ED50 = 0.51. This computer model even predicts the injury threshold for the rather extreme case of 1060 nm, where optical radiation is absorbed very weakly in the cornea (95% of the radiation is transmitted); the model prediction of 6.2 MJ·m −2 compares favorably with the experimental threshold of 8.2 MJ·m −2 for a 5-second exposure and beam diameter of 1.8 mm [18]. Also the case of a very long exposure duration of 100 seconds [38] can be reproduced well with R ED50 values of 1.3 and 1.4 for 2.2 mm and 7.0 mm beam diameters, respectively.
The results also permit the conclusion that a homogenous bulk model of the cornea is appropriate and it is not necessary to consider the layered structure of the stroma nor the birefringence exhibited by the cornea. We argue that it is not necessary to consider the birefringence in the model, as it can be assumed that over the thickness of the cornea any potential lateral displacement of a part of the beam due to birefringence is small compared to the laser beam diameter. This is also supported by the finding that even for penetrating wavelengths and very small beam diameters, the deviation of the model with experimental data is usually small, such as R ED50 = 0.9 for 1318 nm and 0.07 mm 1/e beam diameter [42], or R ED50 = 1.01 for 1338 nm and 0.28 mm [17].
Overall, the lack of a discernible bias of the ratio R ED50 as function of wavelength, exposure duration and corneal beam diameter (see Fig. 5) and other variables such as the end point supports the view that a significant part of the discrepancy between model predictions and experimental thresholds relates to experimental uncertainties.

Nanosecond regime
In the nanosecond regime, no relevant heat flow occurs during the pulse and therefore the injury threshold, when expressed as radiant exposure, is directly determined by the optical absorption depth. For a given radiant exposure, a large absorption depth (low absorption coefficient) implies a large absorption volume and therefore a relative low temperature rise as compared to a shallow absorption depth (high absorption coefficient) and a small absorbing volume. The condition of "no relevant heat flow" applies to pulse durations less than roughly 1 ms even for highly absorbing wavelengths, and for this condition the predicted injury threshold, for a given wavelength, is a constant radiant exposure value without dependence on pulse duration or on beam diameter incident on the cornea. Thus the only dependence is on the wavelength dependent absorption coefficient. As shown in Fig. 6, the computer model predictions -as dictated by the extent of the absorption volume -and the experimental ED 50 data follow the water absorption trend for wavelengths up to about 2.5 µm, i.e. in the range where the absorption depth is larger than 100 µm (absorption coefficient less than 100 cm −1 ). For the wavelength range above 2.5 µm, where the absorption depth is very shallow, the model predictions do not follow the water absorption trend as closely. This is caused by the absorption of radiation within the tear film, which "protects" the cornea and results in a higher injury threshold. The energy absorbed in the tear film exceeds 10% in this wavelength range. For the extreme case at 2.9 µm, nearly 100% of the radiation is absorbed in the tear layer and no optical radiation penetrates through to the cornea. If the tear film were a biological tissue, a lesion would form there, at a lower threshold than in the cornea, which is damaged following conduction of heat from the tear layer. In order to study the effects of the tear layer, the computer model was also run for the assumption that the tear film is absent, and the predictions are "on top" of the scaled water absorption curve. By modelling the tear film with a depth of 9 µm, the model predictions are best in line with most of the available data (optimum in that sense), especially for the CO 2 laser data [46,48] where we see that the experimental data are consistently higher than the scaled water absorption curve. The deviation is largest for the wavelength bands around 2.7 µm and 3.7 µm obtained with chemical lasers [29].
The experimental ED 50 for 1.4 µm and 25 ns-duration exposure [47] reported by Archibald and Taboada is significantly lower than both the scaled water absorption curve and the computed model prediction by a factor of more than 3. The wavelength of 1.4 µm is in a range of very weak absorption, so that the tear film does not play a role. The model can predict all other experimental thresholds in this wavelength range, and these thresholds are fully consistent with the scaled water absorption curve. There is no explanation for the obvious deviation of this threshold by a factor of 3. Based on the evidence provided here we argue that this 1.4 µm threshold does not need to be considered for the comparison with the predictions of the computer model.
No experimental data was found in the range of exposure durations between 250 ns and 100 µs. Since the computer model, being based on thermally induced bulk injury, can predict the injury threshold for exposure durations between 1 ns and 250 ns well, it is possible to conclude that the injury mechanism in that exposure duration regime is of thermal nature and not, for instance, associated to thermomechanical mechanisms as it is in the pigmented epithelium of the retina. Thus, there is no reason to assume that the injury mechanism differs in the range between 250 ns and 100 µs, as compared to the two adjacent ranges for shorter and longer exposure durations, which can be modelled well with the thermal computer model. Therefore, we consider the range of exposure durations between 250 ns to 100 µs to be covered by the validation and the maximum deviations of R ED50 found for other exposure durations should apply.

Lesion depth
Most experimenters have reported the appearance and approximate depth or thickness of the lesions observed by ophthalmoscopic means, whether it was confined in the epithelium or involved the anterior layers of the stroma or the full cornea up to the endothelium. This information is of great importance in the sense that corneal injuries differ in the long-term impact they can have on vision since the epithelium heals within a couple of days [3,46] while lesions involving the stroma and the endothelium persist [18,28,49]. Besides the loss of transparency, deep lesions can imply collagen shrinkage and a decrease in corneal curvature [20].
The diversity of absorption coefficients and exposure durations implies very different temperature profiles through the corneal tissue over time and thermal modelling can efficiently help identifying what emission parameters are likely to induce damage to the stromal layers. Figure 7(a) shows the location of damage onset according to the computer model as a function of exposure duration for a selection of wavelengths and beam diameters. The onset of the injury is the depth at which the damage integral reaches unity over an infinitesimally thin slice with a diameter of 200 µm. For wavelengths below 1900 nm (penetration depth > 0.3 mm), the injury appears unconditionally within the stroma, and in the deep stromal layers towards the corneal endothelium for long exposures. In the short wavelength range (1000 nm to approximately 1200 nm), it is noted that the lesion depth can be smaller than, for instance, at 1320 nm or 1560 nm although the absorption coefficient of the cornea is lower. This is because the absorption coefficient of the aqueous (water-like, see Fig. 2) is significantly lower than that of the cornea. Therefore the heat transfer from the cornea to the aqueous plays a significant role and the posterior temperature of the cornea drops faster than within the stromal layer. Consequently the lesion depth in the wavelength range of 1000 nm to approximately 1200 nm is smaller than at 1320 nm or 1560 nm. For wavelengths between approximately 1900 nm and 2500 nm, a threshold lesion can be confined only in the epithelium for short exposures. For wavelengths above 2650 nm (penetration depth < 100 µm) the stroma could be spared at threshold level even for long exposures. Heat losses at the interface with ambient environment are dependent on a temperature difference between the corneal surface and the ambient environment, thus the higher the temperature increase in the tissue, the higher the losses. It follows -in combination with the Arrhenius model in which the peak temperature to reach the desired effect decreases with exposure duration (see Fig. 7(b)) -that the lesion depth reaches a maximum for exposures lasting approximately 1s to 30s. For longer exposures, where the peak temperature needed to reach the injury threshold becomes smaller, heat losses also decrease and the Beer Lambert's law of absorption leads to a shorter lesion depth. For very short exposure durations the peak temperature is indeed much higher but heat losses do not play a major role (thermal confinement or nearly thermal confinement) and the Beer Lambert's law of absorption leads in the same way to a shorter lesion depth. This effect is mostly noticeable for small beam diameters as shown in Fig. 7(a). In general, the larger the beam, the deeper the lesion as radial losses becomes less important.
This dependence on the absorption coefficient is also reflected in the highest temperature predicted at threshold depth. The higher the absorption depth, the lower the axial conduction (i.e. the longer the thermal confinement for a given beam diameter) and the lower the temperature rise needed to reach the threshold injury level as shown in Fig. 7(b) for short exposure durations. For exposure durations longer than approximately 1 s, the peak temperature no longer depends on the wavelength. The reason for the continued wavelength-dependent lesion depth for long exposures lies in the difference of absorption coefficient and the impact of heat losses at the outer surface of the cornea.

Reduction factors to achieve negligible risk
It was not possible to explain the spread of R ED50 s from the direct comparison of calculated injury thresholds with their experimental counterpart (noting a wide range of exposure parameters with respect to wavelength, exposure duration, etc.) nor from an investigation of secondary parameters (such as assessment delay or animal model). The spread seems appears to follow a normal distribution (r 2 = 0.93, see Fig. 4) regardless of the variable under investigation. Besides the simplification of the damage model compared to the actual response of living tissues, the spread could be related for the most part to uncertainties including experimental uncertainties and to some degree inter-subject variability.
On the basis of this observation of computer model predictions being potentially "too high", the authors recommend a reduction factor when using this computer model for a simplified risk analysis (in a sophisticated quantitative probabilistic risk analysis, the uncertainty of the ED 50 and of the probit slope could be modelled, as was for instance done for the retina [63]). Since model predictions can be up to 1.75 times higher than the experimental ED 50 , reducing model predictions by that factor seems appropriate. This would then scale the predicted ED 50 to the lower edge of the assumed uncertainty band. Additionally, for a simplified risk analysis (avoiding dose-response curves), it is advisable to use a lower probability of injury, such as ED 1 (i.e. 1% probability of injury). Figure 8 shows the lowest determined thresholds from 12 studies [3,28,31,32,[49][50][51] for which both the probit slope and the lowest injury threshold have been reported. In these studies, corneal injuries occurred at levels below the ED 50 but not at levels below the ED 1 . According to the probit distribution model and assuming an ED 84 /ED 50 slope of 1.14 (the geometric average of all slopes reported by experimenters in the literature), the ED 1 level is 1.37 times lower than the ED 50 . The combination of uncertainty (reduction factor of 1.75) and reduced probability of injury (reduction factor of 1.37) leads to an overall reduction of the model-generated injury threshold by a factor of 2.4. Fig. 8. Per experimental study, the lowest injury threshold is plotted (as a relative value to the ED 50 , circles) as a function of the probit slope that was reported for the respective study. Additionally, the reduction relative to the ED 50 necessary to achieve a certain probability for injury is shown as lines; for the example of a probit slope equal to 1.1, the ED 1 (1% probability of injury) is found at an exposure level 20% below the ED 50 .
The value obtained with this reduction factor could be seen as being associated with a very small, and possibly negligible risk for minimally visible corneal injury. A 1% probability for injury might not be seen as sufficiently small to be referred to as negligible risk for injury. However, it can be argued that the log-normal dose-response curve no longer properly reflects the variability of individual thresholds, for exposure levels significantly below the ED 50 (see discussion in [22]). An ED 1 , after all, means that for 100 exposure cases at the given exposure level, only one will exceed the individual's injury threshold and result in a minimally visible lesion. While such a variability of the thresholds can be expected for retinal injury, due to the variability in pigmentation of the retina, it is not likely that such a variability has to be expected for the cornea in the thermal damage regime. Therefore, for an exposure level of ED 1 (which is never determined experimentally due to the high number of laboratory animals it would require) the actual probability for injury can be assumed to be lower as predicted by a log-normal distribution. At some exposure level, it can be assumed that no individual will exhibit a lower threshold. A factor in the derivation of a reduction factor is also the assumed slope of the dose-response curve. It could be argued that due to the reduction to a level of 1.75 below the computer model predictions, a good portion of the variability within the population is already accounted for. Consequently, for the ED 50 reduced by 1.75, a steep dose-response curve could be assumed, such as 1.05. An ED 1 level would in this case already be achieved with a reduction factor of 1/0.9 = 1.1. A reduction of 1.37 applied to the ED 50 is for a slope of 1.05 associated to a theoretical probability of inducing a lesion of less than 10 −10 . Such a line of argument can be used to support the interpretation of the combined reduction factor of 1.75 × 1.37 = 2.4 to result in a negligible risk for injury, which is probably even on the conservative side.
One could also argue that using the overall observed maximum deviation of 1.75 as a reduction factor is too conservative, particularly when the observed deviation for a specific wavelength and exposure duration regime is smaller. It could be possible to support a smaller uncertainty reduction factor, such as of 1.4 for the case of CO 2 laser radiation and exposure durations longer than 1 ms, as plotted in Fig. 3(b) (only exposure durations in the nanosecond regime exhibit deviations of the order of 1.75). However, it is then somewhat questionable if at the same time a steep slope can also be assumed, particularly for a "cautious" approach to a risk analysis. When the slope of the dose-response curve is somewhat larger, then a larger reduction factor is needed to scale the ED 50 value in order to result, with good reliability, in a negligible level of risk. Therefore the overall reduction factor (considering the uncertainty of the model as well as potential variability in the population), for a cautious approach, would be roughly a factor of 2, i.e. not much below the above value of 2.4 (where the more cautious uncertainty reduction factor permitted the assumption of also having accounted for a degree of variability of the thresholds in the population and therefore using a steeper dose-response curve). For a somewhat less conservative overall approach to risk analysis and applicable overall reduction factor, a value of 2 could for instance be justified. This overall reduction factor would result from a combination of an uncertainty reduction of 1.6 and a probability reduction, accounting for variability of thresholds in the population, of 1.25.
It is suggested to take the overall reduction factor as the main parameter, since applying a more conservative uncertainty reduction factor provides the basis to argue for a smaller remaining variability to be accounted for, i.e. a smaller probit slope. Overall the reduction factor would then remain approximately constant. For instance, an overall reduction factor of 2 would also result with an uncertainty reduction of 1.75 and a probability reduction of 1.14, which for a very steep (close to a step function) dose-response curve would still result in a negligible risk for injury.
The ultimate goal of threshold experiments being the development of exposure limits for humans, the computer model presented in this paper can be a valuable tool to further investigate the suitability of the current emission and exposure limits and eventual modifications thereof in the whole range of wavelengths where corneal tissues are at risk. An exhaustive comparison of injury thresholds with the current MPEs is the subject of another publication [64] by the authors.

Conclusion
The limited availability of corneal injury thresholds obtained by experimental means dictates the need for a reliable alternative in order to cover all combinations of laser parameters relevant to thermal laser-induced corneal injuries. The proposed model does meet that need since predicted thresholds are in good agreement with their experimental counterparts, with an overall deviation of 6% and a maximum deviation of less than a factor of 2 over a wide range of wavelengths (1050 nm to 10.6 µm) and exposure durations (1.7 ns to 100 s). For safety purposes, a reduction of the predicted ED 50 by a factor of up 2.4 is suggested to account for model and injury data uncertainties.
Overall, the lack of a discernible bias of the ratio R ED50 as function of wavelength, exposure duration and beam diameter (see Fig. 5) and other variables such as the end point supports the view that a significant part of the discrepancy between model predictions and experimental thresholds relates to experimental uncertainties.

Disclosures
The computer model described in this article is used -besides supporting the improvement of international safety limits -to predict corneal injury threshold for risk analysis studies for laser applications and products, offered by Seibersdorf Labor GmbH as commercial service. The co-authors Bruce E. Stuck and David J. Lund have no relevant financial interests in this article and no potential conflicts of interest to disclose. See Supplement 1 for supporting content.