Suprathreshold laser injuries in excised porcine skin for millisecond exposures at 1070 nm

Abstract. Skin injury response to near-infrared (NIR) laser radiation between the minimum visible lesion threshold and ablation onset is not well understood. This study utilizes a 1070-nm diode-pumped Yb-fiber laser to explore the response of excised porcine skin to high-energy exposures in the suprathreshold injury region without inducing ablation. Concurrent high-speed videography is employed to determine a dichotomous response for three progressive damage categories: observable surface distortion, surface bubble formation due to contained intracutaneous water vaporization, and surface bubble rupture during exposure. Median effective dose (ED50) values are calculated in these categories for 3- and 100-ms pulses with beam diameters (1  /  e2) of 3 mm (28, 35, and 49  J  /  cm2) and 7 mm (96, 141, and 212  J  /  cm2), respectively. Double-pulse cases are secondarily investigated. Experimental data are compared with the maximum permissible exposure limits and ablation onset simulated by a one-dimensional multiphysics model. Logistic regression analysis predicted injury events with ∼90  %   of accuracy. The distinction of skin response into progressive damage categories expands the current understanding of high-energy laser safety while underlining the unique biophysical effects during induced water phase change in tissue. These results prove to be useful in the diagnosis and treatment of NIR laser injuries.


Introduction
Advances in diode-pumped, solid-state, and fiber lasers have brought high-energy applications in the near-infrared (NIR) band to the forefront. The potentially hazardous nature of laser exposure is well established, particularly with regard to possible injurious biological effects to the skin. Thresholds for the transition from healthy skin tissue function to optically induced damage are typically quantified in terms of a laser irradiation level for a given exposure duration and damage criterion. [1][2][3][4][5][6] These thresholds can vary greatly across different types of lasers, as skin and subcutaneous tissues feature different degrees of optical absorption and scattering dependent on the particular wavelength of light involved.
Due to increased deployment of high-powered NIR lasers in medical, industrial, and military applications, there exists significant motivation for accurate characterization and prediction of the suprathreshold response in skin and subcutaneous tissues. At the far, end of this response is ablation, the removal of material from the surface of the skin. In the biological context, there are several models for this action. Verdaasdonk et al. 7 described a three-phase process of tissue ablation using Nd: YAG (1064 nm) and argon (488∕514.5 nm) lasers, consisting of tissue denaturation, explosive vaporization ("popcorn" effect) in conjunction with an increase in light absorption, and then carbonization expanding in a cyclic fashion. LeCarpentier et al. 8,9 observed similar effects, though they only employ argon lasers at 488 and 514.5 nm.
Jacques 10 expanded upon the mechanisms of NIR ablation by highlighting the dynamic role of tissue optics, noting that tissue desiccation and carbonization will form a surface layer of char with drastically increased optical absorption at 1064 nm. This char layer will then heat to extreme temperatures and induce steam formation in lower tissue layers, which will then violently escape through the highly porous char material and take carbonized particles with it, resulting in ablative mass loss. Majaron et al. 11,12 proposed a different mechanism, where the exceptionally large deposition of thermal energy into tissue over a relatively short period of time will result in near-immediate vaporization of water and generation of confined micro-explosions that expel material from the irradiated zone.
Understanding the onset and progression of NIR laser ablation is critical for a variety of applications, such as directed tissue necrosis or removal for medical purposes. A unique characteristic of NIR laser wavelengths, such as 1064 nm in Nd:YAG lasers or 1070 nm in Yb-fiber lasers, is a low optical absorption coefficient for water and thus high depth of penetration in soft biological tissue. In the context of ablation, this feature means that a larger volume of tissue will heat and undergo coagulative damage before the onset of tissue removal. NIR lasers have been employed for soft tissue ablation in the pancreas, 13 thyroid, 14 and liver, 15 where this high-volume heating with accompanying coagulative damage is medically beneficial. These lasers have also been applied for extensive use in hard tissue ablation of dense cortical bone, [16][17][18][19] which features a low water fraction and a greater composition of highly absorbing NIR chromophores. There are currently multiple published studies on laser ablation of skin for surgical purposes, where the goal is to produce a clean cut without substantial adjacent thermal damage. [20][21][22][23][24] However, the NIR band is not optimal for this effect due to the low absorption of water, and thus other laser types are more suitable, such as those emitting wavelengths of 193 nm (argon-fluoride), 20 248 nm (krypton-fluoride), 20 2.94 μm (Er:YAG), 21,22,24 and 10.6 μm (CO 2 ). 23 Several studies have explored exposure of the skin to NIR lasers, but these are typically focused on lower energy therapeutic applications in dermatology, such as nonablative scar treatment, 25 melasma reduction, 26 other forms of laser skin resurfacing, [27][28][29] and laser hair removal. 30 These therapies typically employ pulsed Nd:YAG lasers, where the dose is dynamically adjusted until the desired effect is reached without causing substantial burn damage to the skin. Vincelette et al. 3 performed one of the few investigations involving porcine skin injury thresholds in the NIR band (1070 nm). This comprehensive study included a wide variety of continuous-wave exposure durations (0.01 to 10 s) and laser beam diameters (0.6 to 9.5 cm). As with many laser safety studies, the injurious event of interest was the generation of a minimum visible lesion (MVL), which occurs at much lower irradiance levels than the onset of ablation. Various computational modeling and simulation techniques have been developed specifically for application to laser irradiation of biological tissue, characterizing both the thermodynamic 31-33 and ablative 12 responses. However, severe skin and subcutaneous injuries prior to ablation onset due to suprathreshold, injurious NIR laser exposures have not been well-characterized with experimental data.
Given the aforementioned high depth of NIR laser penetration into soft tissue, incident high-energy laser exposure on the skin will initially heat a larger volume of tissue to a lower peak temperature than other more absorptive wavelengths. As such, visual indications of surface tissue damage must be interpreted differently depending on the type of laser used. Damage due to exposure to a high-absorbing wavelength will be more confined to the surface, whereas exposure to a low-absorbing wavelength that produces a lesion of similar appearance will involve substantially more subsurface and subcutaneous tissue heating.
These distinctions are important in the context of medical diagnosis, particularly in the suprathreshold exposure region prior to ablation onset, where knowledge of the depth of tissue necrosis is critical for proper treatment. 34 The purpose of this study is to investigate skin suprathreshold NIR laser exposure by developing dose thresholds for various damage criteria evident on the surface. These visible injuries may be indicative of more significant tissue damage than the evident surface skin burn, such as subsurface coagulation and deep cellular necrosis, which would require a more complex medical response. This set of experimental data is utilized in concert with the standard recommended skin exposure limits, existing MVL studies, and computational models of ablation onset to provide a comprehensive picture of skin dose-response to NIR lasers.
Traditionally, laser damage studies have been content to report exposure dose thresholds for a specific set of parameters, such as beam diameter and exposure duration, 1-6 and make extrapolations for untested parameters based on observed trends. The possibility of pulsed laser exposure complicates this parameter space, adding terms such as duty cycle and number of pulses. This study provides the expected parameter-specific threshold results but also explores the application of multivariate logistic regression analysis to the dichotomous injury response data. These statistical models, built from the independent variable values and measured outcomes for each experimental exposure, can predict an injury response based on an arbitrary set of input parameters. Though they have not been historically employed in laser damage threshold studies, logistic regression models offer a useful predictive capability for hazard assessment, particularly when multiple laser exposure parameters are involved.

Laser and Optics Configuration
The optical radiation source used in these experiments was an IPG Photonics YLR-3000, which is a Yb-doped, diode-pumped, quasicontinuous-wave fiber laser with a maximum rated power of 3000 W and a central emission wavelength of 1070 nm. This laser has a circular beam profile with a Gaussian distribution. Figure 1 shows the optical configuration utilizing this laser for tissue exposures.
Optics 1 to 3 are high-reflectivity mirrors at 1070 nm, with the small amount of transmitted light through optic 2 being captured by a PM10 (Coherent, Inc.) reference detector. This detector was read by a PM5200 (Coherent, Inc.) power meter to provide an approximate proportional power and energy measurement of each exposure. Optics 4 to 7 were silver mirrors and optic 8 was a 50∕50 beam splitter. These optics were used to guide a helium-neon laser (633 nm) of negligible power onto the target from two directions, one coaxially with the 1070-nm beam and one at a small angle from the side. The intersection of these two beams defined the target plane and ensured proper positioning of the tissue sample.
A beam telescope consisting of two positive lenses and one negative lens (CVI Melles Griot) was used to set the laser spot size to ∼7 mm in diameter (1∕e 2 ). Lens identifiers and optics were as follows: L1 = PLCX-50. 8  path to further reshape this diameter from 7 to 3 mm (1∕e 2 ), when desired.
The profiles of the two laser beam configurations (3 and 7 mm) were analyzed using a Spiricon ThinCam TC-1122 camera, a 10 bit m∕n LBA-710PC frame grabber, and the Spiricon software package. The beams were imaged from the surface of a uniform diffusely reflective Macor ® plate (Corning Inc.), to enable the preservation of their characteristic profiles. Images of these approximately Gaussian profiles are shown in Fig. 2.
It was necessary to determine the relationship between the nominal power setting of the laser and the actual power incident upon the sample. Prior to each day of experiments, the system output at the sample (after all optics) was measured with a PM5K detector and EPM1000 meter over a wide range of input power settings. These measurements were used to generate transformation curves relating power setting and actual output. All equipment utilized for power measurements was up to date with factory-specified, NIST-traceable calibration.

Camera Systems and Exposure Monitoring
The metric of interest for this experiment concerned skin surface changes due to short-duration laser exposures. A high-speed digital camera is particularly suited to these observations, as dynamic changes could not be measured with much certainty through postexposure tissue analysis. In this experiment, a Photron FASTCAM SA1.1 was positioned at a 30-deg angle from the sample plane to observe progressive changes in tissue during laser exposure without being limited to the parallel or perpendicular planes of action. The high-speed camera frame rate was set to 10 kHz with typically 682 frames collected per exposure.
A FLIR SC6800 thermal camera with a spectral range of 3 to 5 μm was positioned in front of the tissue mount such that the image frames were in the plane of the sample surface. The camera was internally calibrated such that its sensor outputs corresponded to temperature values using a Mikron M345-X6-DW black-body radiation source. The emissivity of the excised porcine skin was assumed to be 0.98, similar to human skin. 35 The thermal camera monitored the samples in real time and allow the investigators to enforce a consistent initial sample temperature prior to laser exposure.
A custom LabVIEW program served as the control system for camera and laser activation. Upon initialization of a laser exposure routine, the program waited 500 ms and sent a triggering pulse to the high-speed camera and laser. The high-speed camera employed a rolling buffer, such that this pulse was used to define the start-point of the laser exposure while the first actual video frame was acquired retrospectively at −500 ms.

Sample Preparation and Ambient Temperature Control
The swine model has become a conventional investigatory tool for laser safety studies due to the inherent structural and pigment similarities between porcine and human skin. 36,37 Samples of full-thickness skin were collected from terminated Yucatan miniature pigs from a separate study protocol that allowed tissue sharing. Tissue was surgically excised from the flank, back, and leg regions immediately following euthanasia in a clean necropsy suite, placed in plastic bags, and stored in a −30°C freezer until testing could be conducted. Samples were removed from frozen storage the night prior to testing to defrost at room temperature. After defrosting, skin samples had hair removed by means of plucking to avoid local laser-induced hair-incineration events that would distort skin thermal response. Prior to testing, samples were placed in an incubator for approximately four hours and warmed to 30°C. This sample warming was done to simulate the normal tissue temperature of a live subject. Once a sample reached an adequate temperature, it was mounted in a custom frame within an x − y translation stage and placed in the laser beam path for exposures. Two floodlights positioned below the sample were used as heat lamps, a UtiliTech FU1201A (120 V, 60 Hz, 250 W) and a UtiliTech FU2202 (120 V, 60 Hz, 500 W). The lamps were used to maintain the sample temperature close to the levels expected in live tissue while also providing additional illumination for the high-speed camera.
The temperature of the target area was constantly monitored on the live feed of the thermal camera, with the laser only being activated if the tissue was between 30°C and 32°C, with 31°C being the ideal point. If the measured temperature was over 32°C , the heat lamps were turned off and the tissue was allowed time to cool down. Similarly, if the measured temperature was under 30°C, the heat lamps were reactivated and the tissue was allowed time to warm up.

Experimental Protocol
An exposure consisted of a single-or double-laser pulse sequence with a predetermined set of parameters on a known tissue sample location. These parameters were as follows: sample number, dwell time (total exposure duration), power setting, laser spot size, pulse frequency, and pulse duty cycle. These values can also be used to calculate pulse duration, laser-on time, total energy deposition, and radiant exposure (J∕cm 2 ). Highspeed camera videos were taken for each exposure. A few high-resolution still photographs of the surface and cross sections of skin lesions were taken for select samples during the investigation, but this was not comprehensively performed for each exposure.
A total of 380 exposures were performed, across 13 tissue samples. A laser beam diameter (1∕e 2 ) of 3 mm was used for 94 single-pulse exposures of 3 ms and 119 double-pulse exposures with a pulse duration of 1.5 ms and a duty cycle of 50%. A beam diameter (1∕e 2 ) of 7 mm was used for 109 single-pulse exposures of 100 and 58 double-pulse exposures with a pulse duration of 50 ms and duty cycles of 50% and 80%, respectively. The pulse durations were chosen such that Journal of Biomedical Optics 125001-3 December 2018 • Vol. 23 (12) various degrees of injury could be inflicted and observed over a wide range of tested powers. All of the double-pulse cases consisted of the same total laser on-time as the single-pulse cases, given the same beam diameter. The majority of data collection in this experiment was performed with probit analysis and injury threshold determination in mind, with radiant exposure set as the primary variable. The secondary variables were spot size diameter (3 and 7 mm) and number of pulses (single or double), with pulse duration, duty cycle, and frequency (for the double-pulse sequences) being held constant. The exception was in the case of 7-mm beam diameter double-pulse sequences, where the duty cycle was initially set to 50% for 28 exposures, but there was not enough thermal additivity to observe many injury events. As such, it was changed to 80% duty cycle for a further 30 exposures, resulting in higher thermal additivity and more evident suprathreshold tissue effects.

Ablation Onset Modeling
The onset of skin tissue ablation was computed using an idealized one-dimensional (1-D) implementation of the model proposed by Majaron et al., 12 with the dimension being the thickness of the skin along the vector perpendicular to the surface. The model combines the thermodynamic response of tissue-encapsulated water with the mechanical response of the tissue. Specifically, it predicts the onset of ablation when the water vapor pressure within microscopic bubbles exceeds the ultimate tensile strength of the surrounding tissue. This condition occurs at a known temperature with respect to the material properties, theorized by Majaron et al. to be 314°C based on an ultimate tensile strength of human abdomen skin of 9.5 MPa reported by Duck 38 and experimentally observed by LeCarpentier et al. 8,9 and Verdaasdonk et al. 7 to be between 200°C and 225°C for porcine aorta.
The simulation consisted of a bisection search over an input range of radiant exposure values for a specified exposure duration, with a search tolerance of 0.1 J∕cm 2 . The 1-D model was implemented using a vector of tissue consisting of an 80-μmthick surface epidermis layer, a 2.2-mm-thick middle dermis layer, and an infinitely thick bottom fat layer. The skin thickness values were chosen based on studies of Yucatan miniature pig skin. 36,39 Radiative, convective, and evaporative boundary conditions were imposed on the surface of the epidermis, using an emissivity of 0.98, a convective heat transfer coefficient of 15 W∕m 2 · K, and a relative humidity of 40%. The initial temperature of the tissue was set to 37°C, with the ambient room temperature held at 20°C.
Thermal and mechanical properties of the epidermis and dermis were taken from average porcine skin values compiled by Duck. 38 Optical properties for Yucatan miniature pig skin were measured and documented by Cain et al. 40 over a wavelength range of 1.0 to 1.6 μm. For 1070 nm, the optical absorption coefficients (μ a ) for the epidermis and dermis were 0.6 and 0.5 cm −1 , respectively, and the optical scattering coefficient (μ s ) was 105 cm −1 for both layers. The absorption and scattering coefficients for fat were 0.75 and 100 cm −1 , taken from Salomatina et al. 41 An anisotropy factor of 0.9 was assumed for the three tissue types. A temperature of ablation of 314°C was employed in accordance with the Majaron et al. estimations. 12,38 3 Results

High-Speed Camera Feature Extraction
The high-speed camera, utilized with a frame rate of 10 kHz, enabled observation of dynamic epidermal change on the surface of the tissue, with many frames being recorded during sample laser exposure. The initial expectation was to use these video frames to make a determination of damage with a "yes" or "no" response, corresponding to similar MVL studies for in vivo experiments. 3 However, the investigators judged this method of evaluation to be too simplistic given the range of suprathreshold features evident on the high-speed camera videos. Past investigators have also noted the value of grading laser-induced skin injuries. 42 The large majority of laser-induced changes observed on the samples could be separated into three categories, representing progressive levels of tissue damage. These categories are described below, with accompanying high-speed camera frames. In Figs. 3-5, the left-most image is of a frame prior to laser exposure, with images progressing with time to the right. Each individual frame within these images displays a tissue sample region of ∼5 cm in height and 7.6 cm in width.
(1) Permanent damage. This feature was confirmed if there was a visible difference when comparing the preexposure and postexposure region. In some of these cases, permanent damage according to this definition was clearly observed in the high-speed video frames despite no lesion being evident on the target area upon naked-eye inspection. Figure 3 shows the infliction of a skin surface lesion during a laser exposure that remains evident after the laser has been turned off.
(2) Bubble formation. At higher-energy exposures, distinct spherical bubbles would form underneath the surface of the skin, as a result of water changing state from liquid to gas and collecting between the layers of the epidermis. These bubbles, essentially a steam blister, 43 would deflate once the laser was turned off, as the air and water vapor inside of them cooled down. A bubble formation is shown in Fig. 4.
(3) Bubble rupture. If energy continued to be deposited on a bubble, the outer tissue layer would eventual break down due to the expansion of hot air and water vapor within the bubble or degradation of the tissue from ablation. Hot air and vapor could often be observed exiting the bubble structure at a leakage point. Figure 5 shows images that demonstrate this escape of gas in the circled region; note the slight dark cloud not seen in the previous frames. Bubble rupture could also be distinguished from a deflation event if the bubble size reached a maximum before the laser was turned off, with subsequent size reduction due to escaping gas and not cooling.
In some instances of bubble rupture, the outer tissue layer failed due to rapid expansion of the bubble past a tensile limit, as opposed to small sections being burned away by the laser and causing a gas leak. These "tensile limit" rupture cases were distinguished on the high-speed camera videos by more violent burst characteristics with tissue often failing in a central location, whereas the ruptures due to degradation and ablation of tissue would occur slowly and appear as one or more leaks. High-speed camera frames for a tensile limit bubble rupture can be seen in Fig. 6, with the moment of rupture shown in Fig. 6(b). Although tensile limit ruptures tended to occur at higher powers than leakage-based ruptures, they were more subjective to quantify and did not feature often enough to merit a separate injury category.
High-speed video frames for each exposure were examined by three investigators, with an assessment of tissue damage being made based on these three categories of permanent damage, bubble formation, and bubble rupture. The majority consensus for the presence or absence of the damage category was recorded as a dichotomous response. It should be noted that due to the progressive nature of these damage bins, a positive response for any given category will by definition indicate positive responses for the preceding categories. Likewise, if the features of a given category are not observed, then a negative response for all subsequent categories can be assumed.

Probit Analysis
As the aforementioned tissue damage categories were populated with binary assessments, the probit procedure 44 can be applied to generate a statistical distribution of expected response to a range of doses. The variable of interest to analyze with respect to tissue damage was radiant exposure (J∕cm 2 ), as the likelihood of damage tended to increase with delivered energy for a given beam diameter.   Probit analysis for the permanent damage, bubble formation, and bubble rupture metrics was performed with respect to radiant exposure for spot sizes diameters of 3 and 7 mm, in the cases of single-and double-pulse sequences. Calculated values for the radiant exposure at the effective dose at 50% occurrence (ED 50 ), lower and upper fiducial limits (95% confidence intervals), and probit slope are included in Tables 1-3. These tables also include the number of exposures performed per parameter set with the number of each dichotomous response, beam diameter, pulse duration, number of pulses, and duty cycle. The probit dose and fiducial curves for permanent damage, bubble formation, and bubble rupture are shown in Figs. 7 and 8. All three injury conditions are plotted on the same graph for a given set of laser parameters to compare the ED 50 levels.
Attempted probit analysis for the two-pulse exposures with a 7-mm beam diameter did not converge well due to limited available data. The permanent damage categories for these exposures had too few negative response observations to perform a probit analysis. The bubble formation and rupture categories for the 50% duty cycle cases converged but featured very poor fits to the probit distribution, so only the plots for 80% duty cycle cases are included for this analysis.

Logistic Regression Analysis
The relationship between the laser parameters and damage causality as specified by the three injury assessments was analyzed with a generalized linear regression for a binomial distribution (logistic regression). Logistic regression is appropriate as the measured outcomes are dichotomous decisions ("yes" or "no" evaluations of injury). In this case, the independent variables were laser spot size diameter (mm), duty cycle (as a fraction), number of pulses (unitless), pulse duration (ms), and power (W). The dependent variables of interest were the injury assessments for permanent damage, bubble formation, and bubble rupture. The generalized linear regression model for the binomial distribution fits a logit function, as follows:   Journal of Biomedical Optics 125001-6 December 2018 • Vol. 23 (12) where b is the coefficient estimates, x is the independent variable values, and n is the number of independent variables. The resulting model can be used to calculate the probability of observing a specific injury outcome based on the particular values of the independent variables. Each regression featured 380 observations and five independent variables, resulting in 374 degrees of freedom. The results of logistic regression analysis for permanent damage, bubble formation, and bubble rupture are given in Table 4.
The individual t-statistic for each parameter is the modeled parameter coefficient estimate divided by its standard error and is a metric of the contribution of that parameter to the model variance. The individual p-values test the null hypothesis that the parameter coefficient is zero; a lower value indicates a   Journal of Biomedical Optics 125001-7 December 2018 • Vol. 23 (12) decreased probability of this hypothesis, and thus a greater significance of the parameter coefficient to the model. It can be seen that the logistic regression models for permanent damage, bubble formation, and bubble rupture are all fairly accurate, with total predictive power of 93.4%, 90.5%, and 86.8%, respectively. Prediction accuracy was calculated by running the experimental data through the logistical regression formula. In the calculations of prediction, a positive response was defined as a model probability output ≥0.5, whereas a negative response was defined as <0.5.
For each model, all independent variables are significant parameters at the 95% confidence interval with p < 0.05. In many cases, including the entire case of bubble formation, the independent variables satisfy even more rigid significance criteria with p < 0.0001. As expected, laser power is consistently the parameter with the greatest contribution to the modeled variance, with the absolute value of its associated t-statistic always being the largest among the independent variables.

Discussion
A principal motivation for this study was the determination of robust millisecond-duration laser-induced skin injury thresholds by generating well-defined ED 50 values in terms of radiant exposure with the probit procedure. This study attempted to address this issue by confining the laser parameters to minimally variable beam diameters and numbers of pulses with other factors held constant while sweeping the primary independent variable of power over an effective damage range. A large total number of exposures ensured that the resulting ED 50 s and fiducial limits for each injury category were statistically valid.
Probit analysis has become a standard practice for quantifying the damaging effect of laser exposures on biological tissue due to its ability to accept dichotomous (yes or no) assessments of injury as a dependent variable. However, the technique is limited in that only one independent variable can be used. This variable typically assumes the form of radiant exposure (J∕cm 2 ), which takes into account the exposure time, duty cycle, pulse duration, power, and beam diameter. These individual parameters contribute to tissue damage in different ways, and as such, it is necessary to separate probit calculations into groups according to these other categories.
The maximum permissible exposure (MPE) limits for laser safety, defined as the level of laser radiation to which an unprotected person may be exposed without adverse biological changes to the skin, can be determined by a set of rules and equations that are codified by the Laser Institute of America in the ANSI Z136.1 standard. 45 For Gaussian beams, the MPE calculations are expressed in terms of radiant exposure (J∕cm 2 ) over the 1∕e beam diameter, also known as a peak radiant exposure. In the case of a 1070-nm laser, the MPE is 1.29 J∕cm 2 for a single pulse of 3-ms duration and 3.09 J∕cm 2 for a single pulse of 100-ms duration. To compare the experimental data with the MPEs, the ED 50 values for the permanent damage (Table 1), bubble formation (Table 2), and bubble rupture (Table 3) conditions must be converted to represent the peak radiant exposure thresholds (as opposed to 1∕e 2 diameter, or average, radiant exposure thresholds). The 1∕e diameter is less than the 1∕e 2 diameter by a factor of p 2 for Gaussian beams, so the peak radiant exposure thresholds (J∕cm 2 ) are calculated by multiplying the average radiant exposure thresholds by 2. Ablation-onset thresholds in terms of peak radiant exposure were also calculated for 3-and 100-ms single pulses using the ideal 1-D ablation model described in Sec. 2.5. The values for these exposure durations are provided in Table 5.
For comparison of the experimentally determined injury ED 50 s with other relevant dose thresholds, the single-pulse ED 50 s at the 3-and 100-ms exposure times explored in this experiment were used to determine a best fit for peak irradiance (E) versus exposure time (T) according to a power function of the form E ¼ aT −b . Note that while only two points were used to fit each curve due to the limited number of exposure times  (12) tested, this method follows the same distribution of the MPE specification. 45 The MPEs as calculated by the ANSI Z136.1-2014 standard, 45 the ablation onset thresholds determined by the ideal 1-D implementation of the Majaron ablation model, the MVLs in live porcine subjects determined by Vincelette et al., 3 and the three suprathreshold damage categories explored in this study are plotted together in Fig. 9. In accordance with the standard, the MPE is set to 1 W∕cm 2 for exposure durations >10 s. The values shown in Fig. 9 form a progressive characterization of threshold and suprathreshold exposure response. It can be seen that the MPE dose is very conservative, with the in vivo MVLs being more than an order of magnitude greater than the MPE. The permanent damage injury threshold from this study is slightly above the in vivo MVLs, signifying that there is a physiological injury response in living tissue that cannot be visually observed in excised tissue, such as initial erythema and wheal formation as described by Vincelette et al., 3 or a possible mediating effect of highly-absorbing oxygenated hemoglobin 46 present in vivo but not ex vivo. Bubble formation and rupture events feature at progressively higher thresholds, with fairly proportional separation such that their characteristic slopes are approximately parallel over the experimental range of T ¼ 3 to 100 ms.
The MPE, in vivo MVLs, and excised tissue injury curves all feature similar slopes while the simulated ablation onset threshold curve is noticeably steeper and begins to intersect the excised tissue injury curves at long exposure durations. It is important to note that the simulation did not incorporate the inevitable changes in tissue properties (thermal, mechanical, and optical) and those that occur during the heating period leading up to ablation onset, as these are not currently well documented. The Majaron model itself is also optimized for midinfrared laser ablation, where the absorption coefficient of water is greater. Furthermore, the ablation onset values are generated using a 1-D environment, which does not account for radial heat diffusion. Although a three-dimensional (3-D) or axisymmetric two-dimensional (2-D) environment would be more appropriate for a rigorous study of ablation modeling, particularly at longer exposure durations, the 1-D implementation provided fast results for comparison with the experimental data.
The measured permanent damage, bubble formation, and bubble rupture ED 50 values provide useful insight into suprathreshold tissue response when plotted in Fig. 9, but it must be remembered that this data were generated using two different beam diameters. Vincelette et al. 3 observed that MVL radiant exposure thresholds at equivalent exposure durations tended to decrease as beam diameter increased but for durations of 100 ms or greater and 1∕e 2 beam diameters <2 cm. Further investigations into suprathreshold skin response may prefer to maintain a constant beam diameter and determine progressive injury ED 50 values across variable exposure durations. However, the maximum power of the laser system will limit this range and require a small beam diameter for all exposures. For example, with the optical setup used in these experiments, there was insufficient power to observe the suprathreshold effects at 3 ms for the 7-mm beam diameter. While suprathreshold effects at longer durations could be explored, these may be of less practical interest from a safety standpoint due to expected subject movement in the case of accidental exposure. 45 Although probit analysis will remain an integral part of laser safety studies to characterize tissue damage and generate thresholds, multivariate statistical analysis may also prove to be a useful tool for injury predictions based on specific parameters. A continuous laser exposure can be described by a variety of factors including its beam diameter, duration, power, and wavelength, all of which may contribute to different degrees toward its damaging effect on tissue. These conditions alone provide motivation for a multivariate approach, as the probit technique must group these effects into a bulk radiant exposure term to describe injury causality. However, pulsed lasers occupy an even larger parameter space, introducing independent variables, such as individual pulse duration, duty cycle, number of pulses, and pulse repetition rate.
Logistic regression is capable of generating a binomially-distributed model for laser-induced tissue injury from multiple independent variables and one dichotomous dependent variable. For a given set of exposure parameters and a specific injury condition, the regression coefficients and independent variable values can be used to calculate the probability of observing the injury. This method has been explored in this experiment to make determinations of individual laser parameter significance and contribution, though the regression models themselves were not used to make any prospective determinations. Future work can incorporate validation studies for logistical regression damage models in concert with more complex thermal damage simulations to predict degrees of tissue injury.
Cross sectional and surface photographs of laser-induced skin bubbles were inspected by a veterinary pathologist. An example of these photographs can be seen in Fig. 10. It was determined that the bubbles were formed between the layers of the epidermis, though this could not be definitively confirmed Journal of Biomedical Optics 125001-9 December 2018 • Vol. 23 (12) without proper histology. Future investigations into suprathreshold, preablative exposures should include tissue biopsies with hemotoxylin and eosin staining to isolate the relevant skin layers while also examining the degree of subsurface tissue damage. It is relevant to note that these studies were performed in excised porcine tissue samples that were harvested under tissue sharing protocols. The age and genders of the source swine were dissociated from the tissue samples following excision. These conditions are significant factors with respect to skin layer thickness, which in turn is an important parameter in skin response to optical energy deposition. Furthermore, live tissue will exhibit multiple immediate and latent physiological features including a thermoregulatory response, as it both generates heat through metabolic processes and draws excessive heat away through blood perfusion. Finally, the optical properties of excised tissue are different from in vivo tissue, which has not been damaged by the freezing process and features a greater chromophore fraction of highly absorbing hemoglobin. 46 If it is necessary to make stronger correlations with expected human skin response, future experiments could utilize live swine of consistent species, age, and gender.

Conclusion
This study investigated surface-level tissue responses to highenergy 1070-nm laser exposures under a variety of conditions, including single-and double-pulsed sequences. Progressive injury categories were utilized for determining suprathreshold tissue damage, as opposed to the traditional minimally visible lesion assessment. Subsequent probit analysis generated ED 50 values and fiducial ranges with a high degree of statistical integrity, over multiple levels of injury criteria. A simple 1-D tissue ablation model provided a dose boundary for ablation onset that could be referred to when comparing the suprathreshold ED 50 values with previous MVL thresholds and the standard MPE limits. Multivariate logistic regression analysis of the data predicted outcomes with ∼90% accuracy and could prove to be a useful tool in laser exposure hazard research moving forward, as a generation of experimental results may not be able to keep pace with increasingly complex laser system capabilities. These experiments and data gathered have substantially increased understanding of suprathreshold skin response to high-energy 1070-nm lasers. The results can serve as a basis for dose-correlated prediction of NIR laser injury and support the goal of informing postexposure treatment.

Disclosures
There are no conflicts of interest to disclose.