Modeling Retinal Damage for Threshold Prediction and Probabilistic Risk Assessment

Summary

A modeling approach for laser-induced thermal damage to the human eye is presented. It aims to improve the evaluation of laser-based hazards by providing a means of calculating the damage for a specific irradiation scenario.

Abstract

With the increasing number of laser applications in medicine, defense, and technology, both intentional and accidental exposure of the human eye to laser sources has become a major concern. Model-based prediction of retinal damage thresholds could enable more scenario-specific laser safety assessment, especially for laser parameters not covered by experimental data. Ideally, such models would allow the calculation of ED50 (effective dose at which the likelihood of damage is 50%) values based —among other factors — on wavelength, pulse duration, and spot shape. This requires a detailed understanding and modeling of all damage regimes to reflect the dependency between the key parameters and the dominant damage mechanism.

This work discusses the status of this approach (validated here for the thermal regime, or simply "in the thermal damage regime"); critical aspects that could block its success are highlighted, and potential benefits are outlined. These range from improved accuracy of laser exposure limits in eye safety standards to optimized dosimetry in retinal laser surgery to probabilistic risk assessment for the use of lasers in outdoor environments.

Introduction

This work describes the development and validation of a physiologically detailed thermal damage model of the human eye. Within the thermal injury regime, the model predicts retinal temperature evolution and injury thresholds using an Arrhenius damage formulation (compare explanation in the photothermal damage section). Representative applications cover all situations in which the prediction of ocular temperatures and damage is of interest. This includes, e.g., evaluating damage thresholds for scanned retinal irradiation, understanding the effect of pulse-train additivity on damage thresholds, and comparing computed thresholds with safety limits from the laser safety standard. Outside the thermal regime, current modeling approaches under consideration are presented, and a roadmap is provided toward extending the framework to additional damage mechanisms.

The work presented here concerns the modeling and, therefore, the prediction of retinal damage caused by laser irradiation. While a critical dose could, in theory, always be determined by experiments with animal retinas similar to human retinae, there is a strong need to predict the damage without performing experiments. The laser parameter space of variation (wavelengths, pulse durations, and repetition rates) is vast, implying a prohibitively large number of animal experiments for each new parameter set. Additionally, for long irradiation times, retinal blood flow must be considered as well, which would require in vivo experiments. Consequently, modeling the interaction of the laser with the eye seems to be the only realistic way forward.

The need for a detailed understanding of damage mechanisms and, therefore, damage thresholds (which could be used as ED50 surrogates) is also linked to the situation of the eye safety standard (IEC 60825 or ANSI Z136.1). Because the standard must address the full range of wavelengths, pulse durations, repetition patterns, and spot sizes, it necessarily incorporates simplifying assumptions, interpolations, and conservative safety factors to account for uncertainty. As only a limited number of ED50 values—primarily derived from non-human primate studies—are available, interpolation is required to establish comprehensive MPE (maximum permissible exposure) limits. Although this approach provides broad and practical applicability, a framework grounded in a detailed mechanistic understanding and modeling of damage processes offers clear advantages in terms of physical transparency, scenario-specific accuracy, and applicability without detailed knowledge of the laser safety standard.

For example, pulsed and scanning lasers are evaluated as pulsed sources, although retinal scanning introduces additional temporal and spatial effects. The appropriate treatment of scanning in the derivation of safety limits has been the subject of ongoing discussion within the community over the past decade. Even with regular updates to reflect technological advances, it is not feasible for the standard to cover every complex configuration of new laser systems with new sets of parameters without simplifications and conservative safety factors. Consequently, room for interpretation remains, which may lead to inconsistencies or errors in safety evaluation.

A physics-based modeling approach could substantially reduce reliance on interpolation and conservative margins and extend the standard's applicability to complex scenarios. Because the development and validation of such models are directly linked to a deeper understanding of the underlying damage mechanisms, the resulting insights could also support a more transparent and physically grounded derivation of MPE values from existing ED50 data.

In the long term, eye safety assessment could be streamlined through an integrated, plug-and-play modeling framework. Such a tool could be supplied either with the relevant system parameters or directly with an optical design file (e.g., a Zemax model), which is typically available during product development, thereby enabling consistent and scenario-specific safety evaluation¹.

Another field of application is the growing field of high-energy lasers (HEL), e.g., to counteract drones. Here, the difficulty mostly resides in the laser reflections from targets, especially metallic targets, that can change rapidly and randomly and constitute a hazard to military personnel and civilians^2,3. This random non-deterministic situation requires an adequate evaluation mechanism – usually a probabilistic approach is used. This approach creates statements about the occurrence of certain irradiation situations (intensity, exposure time, wavelength), which still must be translated into a likelihood of damage. Here, the damage modeling discussed in this work can close the gap by translating probabilistic scenarios into likelihoods of damage.

Understanding the principles of damage mechanisms in detail and mimicking these in a software model is a straightforward way to determine damage thresholds instead of performing experiments. Depending on the pulse duration, retinal damage occurs via different laser-tissue interaction mechanisms (Figure 1)^4,5,6,7:

Figure 1: Overview of the damage mechanisms. The type of damage mechanism depends on the exposure time and irradiance. It ranges from photomechanical damage following from very high irradiance at short pulse durations to photochemical damage due to low irradiance for long time durations. Please click here to view a larger version of this figure.

Photomechanical breakdown (10^-12–10^-6 s, 10¹⁰–10¹⁶ W/cm²):
Photodisruption: At ultrashort to nanosecond pulse durations and very high irradiances, nonlinear absorption initiates optical breakdown. A dense plasma of free electrons and ions forms, expanding explosively and generating strong shock waves. The process, called photodisruption, mechanically tears tissue apart with minimal bulk heating and creates sharply defined lesions even in weakly absorbing areas^4,8,9,10.

Plasma-induced ablation: In the ns–µs range at very high irradiance, plasma again dominates. Here, tissue removal is driven not just by shock waves but also by plasma expansion and explosive ablation. The mechanism is known as plasma-induced ablation and produces significant material ejection^4,11.

Photoablation (≈ 10^⁻9–10⁻⁶ s; ≈ 10⁷–10¹⁰ W/cm²):
For pulse durations above the nanosecond and below the microsecond range, photoablation occurs. In this process, molecular bonds are directly broken by the irradiation. This method is, for example, used to correct refractive errors of the eye by reshaping the cornea (the so-called LASIK method). Typically, photoablative procedures are applied in a power range where plasma formation has not yet occurred^4,11.

Thermomechanical damage (≈ 10^-9 – 10^-6 s; ≈ 10⁶ – 10⁸ W/cm²):
In pigmented ocular tissues, strong absorption by melanosomes in the ns–µs range can cause rapid superheating. When the melanosome surface reaches ≈ 150 °C, microbubbles nucleate^12,13. Their expansion and collapse generate mechanical stress waves that damage the retinal pigment epithelium (RPE). This thermomechanical mechanism bridges the gap between photoablation and thermal injury: it is not plasma-driven, but involves mechanical transients coupled to localized heating. Thresholds depend on melanosome size, shape, orientation, and local illumination^4,11.

Photothermal damage (≈ 1 µs– 60 s; ≈ 10– 10⁶ W/cm²):
From microseconds up to seconds, tissue heating dominates. Energy deposition raises temperature, leading first to protein denaturation and, with higher exposures, coagulative necrosis and carbonization. Literature bounds differ: Niemz⁴ cites 1 µs – 60 s, Zuclich¹⁴ 10 µs–60 s. Approximate reciprocity with radiant exposure (~1–1000 J/cm²) holds, with deviations for very short pulses (limited heat diffusion) and very long pulses (enhanced cooling by perfusion)^4,6,15.

The modeling of the thermal damage in the presented work is based on the Arrhenius integral, which is defined as:

   (1)

With τ denoting the exposure time, E_a the activation energy, R the universal gas constant, T the temperature during exposure, and A a pre-exponential scaling factor, the parameters applied in this study are A = 1.3 × 10⁹⁹ s⁻¹ and E = 628 kJ/mol¹⁵. A condition of Ω = 1 is taken to indicate the onset of tissue damage. For a comprehensive discussion of the underlying model, refer to a previous publication⁴.

Photochemical damage (≈ 10 s – 10⁴ s; ≈ 10⁻³ – 10² W/cm²):
For long exposures at low irradiances, cumulative photochemical reactions dominate - such as photopigment bleaching or reactive oxygen species (ROS)-mediated pathways. They occur at exposure durations above 10 s⁶ (or above 1 s⁴). Thus, there exists an intermediate range in which both thermal and photochemical damage can occur simultaneously¹⁶. Chronic low-level blue-light hazard is a typical example.

Protocol

Section 1 describes the general steps necessary for model building, as they can be performed with a variety of software options. Section 2 gives these instructions for the specific and exemplary case where Altair Hypermesh (finite element modeling and analysis [FEA] software) and Ansys Fluent (fluid simulation software) are used. In addition, a supplementary file (Supplementary File 1^{17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45}) is provided, which describes the modeling approach and gives theoretical background in a non-stepwise manner (Figure 1 - Figure 5,  Table 1 - Table 3).

1. Implementation of the modeling approach – General procedural steps

1. Create a 3D model according to the geometry as described above.
   1. Within this work, create the geometry in a 3D mechanical computer-aided design software (Version 16) according to the parameters described in Supplementary File 1.
2. Create the volume mesh according to the properties and extensions of all parts of the eye, according to Supplementary File 1.
   1. Carry out meshing and post-processing with FEA analysis software (Version 11) and computer-aided engineering (CAE) analysis software (Version 11). The model consists of 761,766 tetrahedral elements. The mesh cells have dimensions between 37 µm and 491 µm.
      NOTE: A coarser meshing would result in calculational errors/deviations. A higher number of cells is always possible, but it would increase computational time. This mesh proved to be size-independent, as further refinement of the element size resulted in negligible differences in the calculations.
   2. To maintain reasonable computation times, use this mesh and adopt only if specific regions require a finer resolution. Accordingly, remesh the retina with a minimum feature size of 5 µm to enable absorption within the RPE layer, as required and explained later in the protocol.
      NOTE: Figure 2 shows an exemplary image of the mesh with a hidden vitreous to make the retina visible.
3. Integrate boundary conditions.
   1. At the sclera, use a conduction-only boundary coefficient h_scl = 20 W/m²K and a corneal value h'_corn,1 = 14.21 W/m²K of as derived above in detail. Use these values to define the boundary conditions of the surface representing the sclera and cornea, respectively.
4. Implement the blood flow.
   NOTE: This field supplies flow direction and magnitude in every mesh cell and adds advection to the energy equation.
   1. Follow the construction steps as mentioned below:
      1. Specify physiological inlets/outlets—anterior inflow from the major arterial circle of the iris and posterior inflow from the short ciliary arteries, with venous drainage via the vortex veins.
      2. Compute the resulting flow field.
      3. Rescale all local speeds to a uniform 5 mm  (per Peyman⁴⁵) while preserving directions, implemented via a fluid simulation software user-defined function (UDF). The outcome is a directionally correct velocity field of constant magnitude across the choroid.
5. Calculate the temperature distribution.
   1. Begin the temperature calculation with the definition of the time-dependent positions and shapes of the laser beam. To obtain these time-varying positions on the retina (or within the eye more generally), use a raytracing software model of the eye in parallel with the thermodynamical model.
      NOTE: Based on the raytracing software modeling of the laser system under evaluation, this approach provides the temporal evolution of the laser spot positions and shapes. This is not explained in detail in this work for the sake of accessibility, but is described previously^1,18.
   2. Map the calculated spot-positions and shapes onto the corresponding sections of mesh cells over time, which are then assigned the respective laser power.
   3. Import these data into the fluid simulation software via a UDF, which processes them cyclically. The UDF identifies the relevant mesh cells and introduces an energy source term according to the specified laser power for the corresponding irradiation time. Consider both continuous-wave (cw) and pulsed irradiation, as well as static or time-dependent beam movements.
   4. For each position, calculate the absorbed energy using the absorption coefficient at the given laser wavelength and assign it to the affected mesh cells, which are heated in parallel.
      NOTE: As an example, for irradiation at a wavelength of 532 nm, 51% absorption occurs in the RPE with a thickness of 5 µm⁷, as shown in Figure 3. No absorption occurs within Bruch's membrane^7,46. The remaining laser energy decreases in the choroid, assumed to be 100 µm thick with an absorption coefficient⁷ of α_Ch = 270 cm^-1, according to the Lambert-Beer law. The UDF calculates the absorbed energy across the thickness of each mesh cell and assigns this value to the cell. For the transmission from cornea to retina, a value of 20 % is applied⁷. A change in laser wavelength is easily accommodated by simply adapting the absorption coefficient and the transmission from cornea to retina.
   5. In addition, define the duration of a single time step, since the spot-shape data (from the raytracing model or defined directly in the UDF if it is known without the raytracing model) only describes the irradiation sequence without a temporal scale.
   6. To accelerate the calculation, group all intervals without irradiation together and assign a correspondingly larger time step, as no fine temporal resolution is required there.
      NOTE: After these settings are applied, the solver can compute the temporal thermal response of the eye. Since the present model only requires heat conduction and blood flow to be considered, the solver is restricted to solving the corresponding differential equations. Refer to Supplementary File 2 for exemplary user-defined functions.
6. Calculate the damages.
   1. To derive a damage threshold based on the thermal response, extract the calculated temperature values from the Fluent result file using a C++-based tool (see^1,18 for details and code listings).
   2. Following the Arrhenius integral approach, calculate the Arrhenius values of the individual mesh cells from these temperatures and store them in a format compatible with the post-processing software. This enables convenient switching between temperature values and damage probabilities during analysis.
   3. Finally, obtain a binary damage outcome (yes/no) for each mesh cell. Determine the yes/no classification by a threshold value of 1 for the Arrhenius integral⁴. Identify this damage threshold iteratively by testing different input powers for a fixed irradiation duration.
7. Set solver-timestep size.
   1. Use the fluid simulation software (Version 14.5) as the solver, employing the finite volume method. The step size shall satisfy the Courant-Friedrichs-Lewy condition, which describes the ratio to the size of the mesh cells in order to create stable simulation conditions.
   2. Also, reduce the step size step-by-step until the result becomes invariant to the step size. E.g., for calculations with up to 10 s time steps of 1 ms and 100 ms, the temperature changes were found to be <1%.

Figure 2: Section cut through mesh (vitreous humor masked out). The figure shows the sclera (white), the choroid (red), the retina (yellow), the lens (white), the iris (green), and the aqueous humor (blue)¹⁸. Please click here to view a larger version of this figure.

Figure 3: Absorption characteristics at the ocular fundus. Within the RPE, 51% of the 532 nm laser radiation is absorbed; no absorption occurs in Bruch's membrane, and the choroid is modeled according to Lambert-Beer behavior¹⁸. Please click here to view a larger version of this figure.

2. Illustrative example – Specific steps

1. Create a 3D model according to the geometry as described above.
   NOTE: The 3D model is built by creating objects for all elements listed in Supplementary Table 1. As an example, these steps will lead to the sclera.
   1. Create a new part file (.ipt).
   2. Start a 2D sketch on the XY-plane.
   3. Draw a horizontal construction line (this will be the rotation axis).
   4. Draw an arc with a radius of 12 mm for the outer surface.
   5. At the center (posterior pole), create a point 0.99 mm inward from the outer arc.
   6. From this point, sketch the inner surface curve.
   7. At the limbus, set the distance between the inner and outer surface to 0.75 mm.
   8. Connect inner and outer curves at the edge to form a closed profile.
   9. Select 3D Model > Revolve.
   10. Choose Surface and confirm.
2. Create the volume mesh.
   NOTE: The mesh is created based on the 3D object. This is an example for meshing:
   1. Import the 3D model created in the 3D mechanical computer-aided design software.
   2. Click on the surface of the object to mesh, e.g., the sclera.
   3. Mesh the surfaces first – click Mesh in this interface, which appears once the surface is clicked upon. Select surfs in the interface and select automatic adaptive meshing.
   4. The 3D mesh (defining the object) is based on the surface meshes. To create, select the surfaces and click from the top left menu: Mesh > Create > Solid Map Mesh.
   5. Now, an interface appears. Click mesh.
3. Integrate boundary conditions.
   NOTE: Boundary conditions have to be defined for the sclera and cornea, since these surfaces are in contact with the surrounding world. This example is for the sclera:
   1. Go to Analysis from the top menu > BCs > Create > Constraints (or Loads).
   2. Select the entity type surfaces and pick the sclera from the menu.
   3. In the same panel, locate the value input fields and enter the value directly in the magnitude field.
   4. Click create to assign the boundary condition.
4. Implement the blood flow.
   1. Define surfaces for - anterior inflow (major arterial circle of iris), posterior inflow (short ciliary arteries), and outflow (vortex veins).
   2. Assign boundary conditions for these surfaces.
   3. Click on Analysis > BCs > Create > Loads.
   4. Select inlet surfaces > assign velocity direction vectors (approximate physiological directions).
   5. Select outlet surfaces > assign pressure outlet (or zero constraint depending on solver setup).
   6. Export model to the fluid simulation software – this is where the length of the flow vector will be adapted.
   7. Use the ScaleVelocity UDF from the annex within the fluid simulation software to adapt the length of the vector (UDF is a standard format to be loaded into the fluid simulation software).
5. Calculate the temperature distribution (Set Solver).
   NOTE: Use the fluid simulation software to calculate the temperatures like the following and save it in a file format readable by the CAE analysis software.
   1. Click File > Read > Mesh > Select mesh > Mesh > Check.
   2. Click General > Solver: Pressure-Based > Time: Transient.
   3. Click Models > Energy > Enable > OK.
   4. Click Materials > define ρ, p, kρ > Cell Zone Conditions > assign materials to all regions.
   5. Click Boundary Conditions > define thermal BCs (cornea, sclera, environment).
   6. Click Solution Initialization > Initialize.
   7. Click Run Calculation > set Time Step Size + Number of Time Steps + Iterations per step.
   8. Click Run Calculation > Calculate.
   9. Click File > Export > Solution Data… > File Type: Ensight Gold > select Temperature > select zones > Write.
6. Calculate the damage.
   1. Execute the AddingArrhenius.exe (Supplementary File 3 and Supplementary File 4) from the same folder as the FEA analysis software files containing the temperatures. It contains a simple calculation of the Arrhenius integral based on the temperatures stored in the FEA analysis software files.

Results

This section shows the verification and validation of the presented model first. Afterward, three exemplary applications are demonstrated.

Verification and validation
In this section, the presented model is first verified by comparison with established thermo-physiological eye models^{19,27,28,30,32} to ensure consistency of the predicted temperature distributions (steady-state situation). In a second step, the model is validated by comparing simulated temperature fields and Arrhenius damage predictions with experimental data obtained under laser irradiation (transient temperature distribution).

Validation of the blood-flow implementation is more challenging, as quantitative in vivo human retinal temperature data under long-duration irradiation are not available. For short irradiation times (several hundred milliseconds), temperature data from retinal surgery can be used for validation. These measurements show good agreement with the predicted temperature rise. As mentioned before, the impact of the blood flow becomes only relevant after a few seconds.

For this long-irradiation regime, only animal data are available. Comparison with retinal temperature measurements in rabbits shows qualitative agreement, while quantitative differences are observed. These deviations are expected due to species-dependent differences in blood flow.

Thus, although the vectorial implementation of blood flow provides a more physiologically realistic representation and a mean velocity of 5 mm/s constitutes a reasonable initial estimate, further validation against human or non-human primate (NHP) data is recommended if such data becomes available.

Model verification – Comparison with established models (without blood flow) (steady state)
The steady-state analyses performed to validate the model used a time step of 10 ms. It can be seen that the temperature distribution along the model is consistent with other models. The other models do not have a blood flow; accordingly, the blood flow in the model described here was deactivated for this verification. In the next step, the blood flow was then activated and compared successfully to in vivo measurements, compare results section (Supplementary Figure 3).

Supplementary Figure 3: Comparison of temperature distribution with other models¹⁷. The temperature within the human eye along the optical axis, following from the presented model, is in good agreement with existing models. Please click here to download this file.

Model validation – Retinae in water (mimicking vitreous body)
Since infrared measurements are impossible under water, a fiber-optic temperature sensor (GaAs probe) was inserted into a porcine eye sample of sclera, choroid, and RPE. The sensor was placed between the sclera and choroid and fixed in place using a holder (Figure 4). The measurement method using a GaAs probe was chosen because the use of an optical fiber to transmit the signal ensures very low thermal coupling, thereby minimizing the influence of the measurement method on the measurement itself.

The setup consisted of two optical paths, schematically shown in Figure 5. On the one hand, the laser passed through a beam splitter and an anti-reflection-coated entrance window to irradiate the tissue; on the other hand, light backscattered from the tissue was directed via the beam splitter and a 20 mm lens to a camera. This camera was primarily used to monitor the spot position during measurements, allowing the laser spot's position on the tissue and relative to the GaAs tip to be checked in real time. To improve the visibility of scattered light from the tissue, a polarizer suppressed reflections from optical surfaces. For optimal imaging, the tissue was additionally illuminated from behind, since the front part is strongly absorbing RPE; at the same time, this transmission approach allowed the sensor's position to be identified. Illumination was provided by a light-emitting diode.

For the measurements, the initial position was set such that the laser spot was exactly centered on the sensor tip. From this starting point, the tissue was moved laterally out of the laser spot while the sensor tip moved with it. As a result, the temperature response at different distances from the spot center could be recorded. Later, these measurements at various distances were compared with corresponding simulations. Therefore, the tissue holder was shifted along the vertical axis while the laser, beam splitter, and imaging optics remained fixed. The tissue was immersed in a water tank maintained at a constant 22 °C by a heating coil. A reference probe controlled this value and simultaneously served as a comparison point for the sensor reading in the tissue before irradiation was applied.

Figure 4: Measurement setup. The tissue sample is mounted inside the holder, including the GaAs measurement tip¹⁸. Please click here to view a larger version of this figure.

Figure 5: Experimental setup for underwater measurements. The GaAs measurement tip was inserted into tissue fixed in a holder. The assembly was immersed in deionized water at 22 °C and irradiated with the laser. To monitor the laser spot position relative to the measurement tip, the irradiated area was imaged onto a camera via a beam splitter placed in the laser beam path¹⁸. Please click here to view a larger version of this figure.

Irradiated retinae over time – Model vs. measurements in vivo (active blood flow)
This section presents the results for long-term retinal irradiation at 810 nm up to 10 s where it was compared against measurements in vivo on rabbits. As this is the first time such simulations are performed and no human measurements are available for such long time spans, the results are compared with rabbit data reported by Herrmann et al.⁴⁷. Figure 6 shows simulations for powers of 62 mW and a 2 mm spot-diameter using a time step of 10 ms (more results can be found in previous publications^1,17). While the overall trends are similar, the absolute temperature levels differ. This discrepancy is attributed to the higher blood flow in the human eye and the anatomical differences between rabbit and human eyes, as well as the higher absorption in the rabbit eye⁴⁶. Importantly, the point at which the curves for active and inactive blood flow diverge agrees well with the present model. Moreover, the relative differences between retinal temperatures in dead rabbits (no blood flow) and living rabbits (normal blood flow) are consistent with the relative differences predicted by the simulations.

Also, in vivo experiments on humans are available for shorter time spans and were used to test the model against. Figure 6 shows on the right-hand side the retinal temperature profiles under laser irradiation with 532 nm at 200 mW and a spot diameter of 300 µm as measured by Brinkmann et al.⁷. These measurements are compared with the temperatures calculated using the toolchain. Although the experimental data exhibit considerable fluctuations, the model is seen to predict the temperature rise with good accuracy. The comparison with measurements from rabbits' eyes is done here, because there are no in vivo experiments measuring the temperature of the retina under irradiation available for humans or NHPs outside the ms-regime.

Figure 6: Comparison model vs. measurements (in vivo). Left: Calculated temperature increase for the human eye ('model') versus measured temperature increase of a rabbit's eye - 62 mW / spot-diameter = 2 mm. Right: Temperature predictions compared to measurements (spot-diameter = 300 µm). This figure was reproduced with permission from Heussner et al.¹⁷. Please click here to view a larger version of this figure.

For the irradiation scenarios described in the left image here, a step size of 10 µs was applied, whereas for the long-duration irradiations in the right image a time step of 10 ms (Using time steps of 1 ms and 100 ms instead resulted in temperature variations of less than 1% after 10 s) was chosen.

Inferring human retinal injury thresholds from rabbit and ex vivo models obviously introduces uncertainties. Choroidal blood flow differs across species, and ex vivo tissue lacks perfusion entirely—both factors influence thermal load dissipation and can over- or underestimate damage for long-duration exposures. Absorption coefficients vary with melanin in the RPE/choroid. As a result, this validation approach has certain limits. For the field of thermal damage, these differences become less critical, though. While the blood flow cannot be transferred from, e.g., a rabbit to a human (the comparison below only considers the qualitative behavior), the matter of layer thicknesses and detail sizes has little effect, since the thermal properties of the different tissue layers are similar.

Irradiated retinae over time – Model vs. measurements ex vivo retinae (blood flow deactivated)
Figure 7 compares the measurements obtained under water with the calculated temperature distributions at the choroid–sclera interface. Results for 50 mW are depicted, but also results for 100 mW and 200 mW were used to ensure that the agreement between measurements and simulations is independent of power, see previous publications^18,36 for more details. For each power level, the temporal temperature profiles were compared for fiber-optic sensor positions at 0 mm, 1 mm, and 2 mm from the laser spot center. The observation period of 60 s corresponds to the timescale of thermal damage. The time step size for 60 s calculations was 100 ms, test runs with shorter time step sizes showed differences in the range of 1%. In general, the mesh size and step size can be adjusted to the specific needs of the model. If shorter calculation times are required and reduced accuracy of the eye geometry or of the calculated temperature is acceptable, an increase in mesh cell size and/or time step might be advisable.

Figure 7B shows a comparison between measurements and simulations of the final temperatures after 60 s as a function of distance. This demonstrates that the spatial distribution of the calculated temperatures also matches the experimental data.

Figure 7: Comparison model vs. measurements (ex vivo). Left: Temperature measurements in water taken with the GaAs tip: The solid curves show the measurement results for a distance of 0 mm, 1 mm, and 2 mm between the laser spot and the GaAs tip. The dashed curves represent the respective simulations with (upper curve) and without (lower curve) blood in the choroid. Laser power: 50 mW, spot-size: 1.9 mm ( n = 14 ). Right: Lateral temperature profile of the measurements in water (solid lines) compared to the simulations (dashed lines) with (upper curve) and without (lower curve) blood in the choroid. The average is shown in orange, the error bars represent the standard deviation. Laser power: 50 mW, spot-size: 1.9 mm ( n = 14 ). This figure was reproduced with permission from Heussner et al.¹⁷. Please click here to view a larger version of this figure.

Thermal damage threshold prediction
The temperature over time for simple laser irradiation is shown in Figure 8:

Figure 8: Temporal development of retinal temperature. The temperature distribution on the retina after 5 ms (left), 7.5 ms (middle), and 10 ms (right) (grid size = 5 µm). This figure was reproduced with permission from Heussner et al.¹⁷. Please click here to view a larger version of this figure.

The Arrhenius value following from this is depicted in Figure 9:

Figure 9: Temporal development of retinal damage Arrhenius values distribution on the retina after 5 ms (left), 7.5 ms (middle), and 10 ms (right) (grid size = 5 µm). This figure was reproduced with permission from Heussner et al.¹⁷. Please click here to view a larger version of this figure.

Table 1 shows the predicted damage thresholds in comparison with measured values. The reference measurements¹⁵ defined damage as the appearance of a visible 20 µm lesion, which was adopted as the criterion here. Accordingly, the required power to reach an Arrhenius value of at least one was calculated for a 20 µm diameter at the specified irradiation time. For validation, the ambient temperature was adapted¹⁵ from 25 °C to 23 °C, and a top-hat laser profile was assumed. The resulting thresholds show good agreement with the experimental data.

Beam diameter [µm]	Irradiation time [ms]	Measured damage threshold [µJ]12		Calculated damage threshold [µJ]
120	1	85		86
	10	241		241
	100	1362		1294
288	1	456		497
	10	1212		1139
	100	4062		3697

Table 1: Comparison between calculated and measured damage thresholds¹⁷.

Thermal damage threshold prediction – dynamic irradiation (retinal scan)
When a scanning laser enters the human eye, the retinal irradiation pattern over time depends on both the distance between the scanning mirror and the eye, as well as on the accommodative state of the eye^48,49. A detailed analysis of these dependencies based on an additional optical eye model (Zemax) can be found in a previous publication¹. Here, we illustrate the principle by means of an exemplary retinal irradiation scenario. In the example shown in Figure 10 (left), the red point marks the center of the laser spot (blue circle). The scan proceeds sequentially from the upper left to the lower right, with a step size of 10 µs. This was again realized by a user-defined function writing energy depositions into the relevant mesh cells over time and evaluating the resulting temperature over time for all cells using the Arrhenius approach. In this example case, this means that each retinal location is exposed for 10 µs at a radiant power of 1 W. The corresponding thermal response at the point of maximum temperature on the retina, together with the derived Arrhenius integral, is depicted in Figure 10 (right). Under these conditions, the Arrhenius evaluation indicates the onset of retinal damage after approximately 150 ms (total scan duration) of exposure to the 1 W scanning laser.

Figure 10: Exemplary damage assessment. Left: Exemplary irradiation scenario on the retina. Right: Respective temperatures and Arrhenius values over time¹⁸. Please click here to view a larger version of this figure.

Evaluation of pulse additivity and its influence on the damage threshold
The presented thermal damage modeling approach was also applied to investigate the influence of different pulse patterns on retinal damage induction. In particular, the model allows a detailed analysis of pulse additivity effects, which are of central importance for eye safety evaluations involving repetitive or modulated laser emissions.

The simulations performed for pulse durations governed by thermal damage mechanisms support the concept of a "partial energy" or "partial N." In this framework, the effective number of pulses is no longer determined solely by counting discrete physical pulses. Instead, pulses are weighted according to their relative energy contribution within the emission sequence and their temporal overlap with the thermal relaxation time of the tissue. Consequently, the damage threshold depends not only on pulse count but also on pulse energy distribution, repetition rate, and thermal accumulation effects between successive pulses.

This interpretation provides a more physically consistent description of pulse additivity in the thermal regime. Rather than applying a purely geometric pulse-counting rule, the model evaluates cumulative temperature rise and the resulting Arrhenius integral directly. The concept of "partial N" therefore generalizes the current pulse-addition definition by linking it to the underlying thermal accumulation processes within the tissue.

This approach enables systematic investigation of how variations in pulse duration, inter-pulse spacing, modulation depth, and total radiant exposure influence the effective damage threshold. Since pulse additivity is a complex topic with significant implications for the laser safety standard, the reader is referred to previous publications^15,50,51 for a more detailed discussion.

Comparing the laser safety threshold with the calculated damage threshold
The presented model was used to derive and compare damage thresholds with class 1 limits⁵² following from IEC 60825-1. The calculated thresholds were found to be consistent with the implicit safety margin incorporated into the standard.

A key advantage of the modeling approach is that it is directly based on predicted damage thresholds (combined with a safety factor). In contrast, the laser safety standard relies on generalized exposure limits derived from experimental data and interpolations across wavelength and time regimes (Table 2). While this ensures conservative applicability across a broad range of systems, it may not always reflect the specific optical and temporal characteristics of a given device and has the potential of being too conservative.

The modeling-based evaluation, therefore, provides a physics-driven alternative that can complement standard-based classification. In addition, the formal application of the laser safety standard can be complex and often requires detailed expertise, particularly for systems involving scanning, pulsed, or modulated emissions. An automated modeling framework integrated into optical design workflows could facilitate early-stage safety assessment and enable systematic optimization of laser systems.

	Laser safety standard		Proposed model
Spot diameter (retina)	C6	Allowed emission	Damage threshold	Reduction factor	Allowed emission
50 µm	1.96	2.0 mW	5.89 mW	3	1.96 mW
100 µm	3.92	4.0 mW	10.24 mW	3	3.41 mW
250 µm	9.8	9.8 mW	28.86 mW	3	9.62 mW

Table 2: Calculated damage thresholds versus class 1 limits from IEC 60825-1. This table was reproduced with permission from Heussner et al.⁵².

Supplementary Figure 1: Geometry used. This model, its extensions, and tissue parameters were derived from the actual biological parameters and is the blueprint for creating the three-dimensional model of the human eye. Please click here to download this file.

Supplementary Figure 2: Boundary conditions used for eye models. The left-hand side shows state-of-the-art boundary conditions with a circular definition of   while on the right-hand side  is defined along the correct eyelid extensions. This figure was reproduced with permission from Heussner et al.³⁶. Please click here to download this file.

Supplementary Figure 4: Scheme of the physiological blood flow. The blood flow enters via the circulus arteriousus iridis major and the arteriae ciliares posteriores longae and exits via the venae vorticosae. This figure was reproduced with permission from Heussner et al.¹⁷. Please click here to download this file.

Supplementary Figure 5: Bloodstream implementation. Based on the actual blood flow, the bloodstream was modeled with the presented inlets and outlets, as shown in Supplementary Table 3. This figure was reproduced with permission from Heussner et al.³⁶. Please click here to download this file.

Supplementary Table 1: Geometrical values used in this work. Please click here to download this file.

Supplementary Table 2: Tissue parameters used within this work. Please click here to download this file.

Supplementary Table 3: Definition of inflow and outflow for the vectorial bloodstream. Please click here to download this file.

Supplementary File 1: Derivation of the modeling approach. Please click here to download this file.

Supplementary File 2: Exemplary user-defined functions. Please click here to download this file.

Supplementary File 3: AddingArrhenius.exe. Please click here to download this file.

Supplementary File 4: AddingArrhenius.cpp. Please click here to download this file.

Discussion

Within the protocol, there are no critical steps with respect to the creation of the model. Critical is the selection of the mesh size and type in combination with the time step size. A modification of the modeling technique is relevant with respect to employing other scanning patterns on the retina. The model itself shall not be modified. The technique is limited to the damage regimes as explained in the respective section. The significance lies in the possibility of damage prediction without animal experiments. One of the most relevant future applications is the probability analysis for retinal damage in scenarios involving HEL.

Once the model is validated for a certain damage regime (as it is now for thermal damage), it can be used to predict damage thresholds, as an experiment would provide an ED50 value. This value can then be used to either derive eye safety calculations instead of using the eye safety standard. For this, the ED50 values must be combined with a certain safety factor, as it is already assumed in the IEC 60825 standard. The difference is that since the ED50 value is now specifically known for all sets of parameters, the safety factor can be constant over all cases and does not need to reflect uncertainties with respect to the ED50 value itself; it only needs to reflect the likelihood of damage resulting from scaling the ED50 by a safety factor. The specific definition of this safety factor should be discussed within the laser safety community.

For other cases, where the standard is not applicable due to the unpredictable nature of the scenario, the combination of a probabilistic hit model with the presented damage model seems to be the best solution. A good example here is the evaluation of high-energy lasers for outdoor applications.

While it is obvious to adhere to legally mandated occupational safety and health regulations to safeguard operators and uninvolved third parties from unintended harmful exposure, it can be quite challenging to implement these measures in outdoor environments. This, in turn, imposes restrictions that limit experimental possibilities, especially when it comes to high-energy lasers (HEL) in the defense sector. These types of lasers can have 100 kW output power, typically at 1 µm wavelength and with excellent beam quality and low divergence, which means that not only the direct irradiation of humans must be considered, but also non-deterministic situations resulting from scattering, reflection by any kind of object in the path of the beam. Hazards can have three origins, as indicated in Figure 11. The direct beam, the atmospheric scattering, and reflections from the target. Another scenario, not described by the figure, would be the reflection of the laser on the sea surface, which could happen in a maritime scenario⁵³.

Not every parameter of these scenarios can be deterministically predicted. For instance, atmospheric turbulence is by nature a stochastic phenomenon and impacts the beam shape and its position on the target. The target itself and the interaction of the laser with the target are also never completely deterministic, especially when a metallic target is melting². Consequently, the reflected laser power, direction, divergence, and shape can vary rapidly. Different groups are performing experimental investigations on these complex irradiation scenarios^{2,54,55,56,57} and use different approaches for their analyses. However, no consensus exists yet on a proper methodology to translate the results of these experiments into a laser safety assessment tool. For these scenarios, the input parameters used to perform a laser safety assessment will have to be probabilistic and might be described using a probability distribution function. From there, a Monte-Carlo simulation could be used to calculate all possible resulting ocular hazard distances (OHD). The worst-case set of input parameters, i.e., the set resulting in the highest NOHD, can then be assumed to define a laser hazard area during a trial.

However, in a dynamic scenario with a moving target in altitude, it is not always the reflected beam with the largest OHD that exhibits the highest risk. The risk is usually defined as the combination of the probability of exposure and the extent of damage when exposed. A large OHD usually relates to a beam with a high power and a low divergence, for the reflected beam, which would imply a small beam diameter with large intensities on the ground. However, the probability of a person being exposed to a reflection exhibiting a small beam diameter might be lower than for a beam with a large diameter. Moreover, the exposure duration might be shorter for a small beam than a larger one. Consequently, it might result in a smaller overall risk despite a larger OHD. If exposed by the beam, a validated thermal model of the eye as presented in this paper is a crucial part of evaluating the potential extent of the damage to an exposed person and, therefore, evaluating the risk. Additionally, because a person does not necessarily directly look at the target during the engagement, reflected laser irradiation enters the eye at an angle and could be imaged far from the fovea. With a thermal model, one can therefore distinguish between severe and low damage depending on the size and position of the retinal damage. These considerations would greatly improve a risk assessment for uninvolved third parties in a laser engagement, similar to what already exists for conventional weapon systems.

Figure 11: Potential hazard zones during high-power laser application. The use of high-power lasers creates potentially hazardous situations due to direct beam irradiation, (partial) reflection from the target, and atmospheric scattering. Please click here to view a larger version of this figure.

Targeting a model that covers all ocular damage in the future, the approach has to be extended to other damage regimes. The first steps towards modeling of thermomechanical damage is described here.

For the modeling of thermomechanical damage, we suggest using the surface temperature at the melanosome as the relevant quantity. A nucleation temperature of 150 °C could be used as the threshold for bubble formation, as suggested previously^12,13. Based on the working hypothesis that bubble nucleation coincides with retinal injury, a melanosomes surface temperature of 150 °C is therefore taken to indicate the onset of thermomechanical damage (Figure 12).

Figure 12: Exemplary mesh of a melanosome. This mesh shows one option to model melanosomes (left) and an exemplary temperature distribution under irradiation. This figure was reproduced with permission from Heussner et al.⁵². Please click here to view a larger version of this figure.

In this approach, a homogeneous laser energy distribution (top-hat profile) is assumed whose size exceeded the dimensions of the melanosome, ensuring complete irradiation of the particle. The melanosome dimensions were set to 2.5 µm by 1 µm, and the absorption coefficient⁵⁸ to 13.000 cm⁻¹. Under the assumption of Lambert-Beer absorption throughout the melanosome. Calculation of the resulting surface temperature over pulse durations from 1 ns to 10 µs yields the damage threshold presented in the Results section.

The calculated surface temperature of melanosomes for pulse durations between 1 ns and 10 µs yields the damage threshold shown in Figure 13. Comparison with available experimental data shows good agreement for most data sets^12,59,60,61. However, accurate modeling requires a more precise definition of the damage thresholds, which necessitates the collection of new experimental data. It remains essential to determine whether the observed lesions are caused by thermomechanical or purely thermal mechanisms, as these pathways require different modeling approaches that must be validated by appropriate measurements. One possible strategy is the detection of shock waves, which are likely to accompany thermomechanical damage processes.

Figure 13: Computed and experimental thresholds for thermomechanical damage using the 150 °C criterion. First results from the melanosome modeling are compared with damage thresholds from the literature and show good agreement. This figure was reproduced with permission from Heussner et al.⁵². Please click here to view a larger version of this figure.

To summarize, the next steps towards modeling of thermomechanical damage will be to evaluate whether the 150 °C approach is, in general, sufficient to model the damage. Further analysis via comparison against robust experimental data is required. This analysis might confirm the approach or conclude that a deeper understanding of the damage mechanism is required, which would translate into the need for a model covering these aspects. Another aspect relevant to this is the creation of shock waves due to laser absorption and the comparison of their pressure with the pressure stability of cell membranes.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Altair Hypermesh 11.0	Altair Engineering Inc.	Version 11	Mesh creation and display; finite element modeling and analysis [FEA] software
Ansys Fluent 14.5	Ansys Inc.	Version 14.5	Thermodynamical solver; fluid simulation software
Autodesk Inventor	Autodesk	Version 16	3D mechanical computer-aided design software
Hyperview	Altair Engineering Inc.	Version 11	computer-aided engineering (CAE) analysis software
Optic Studio 13 (Zemax)	Zemax Development Corporation, today it is Ansys Inc.	Version 13	Raytracing Software
Optocon FOTEMP2	Optocon	https://comem.com/en/optocon/	Spectrometer; Temperature measurement of ocular tissue
Optocon TS2	Optocon	https://comem.com/en/optocon/	Measurement probe; Temperature measurement of ocular tissue