2.1.

Simulation Framework for Thermal Damage Estimation

To determine damage, tissue-photon interactions were first modeled during laser energy delivery using three-dimensional (3D) Monte Carlo simulations.¹² These simulations were performed with the optical properties summarized in Table 1, which are specific to swine liver tissue characteristics at a wavelength of 750 nm, for later comparisons with experimental results. The optical absorption coefficient dictates the depth of photons penetrating tissue prior to being absorbed. The scattering coefficient dictates how the tissue scatters photons outside of the laser beam, and the anisotropic factor determines the amount of forward direction propagation retained after scattering. These parameters were modeled within a $20\times 20\times 20{\mathrm{mm}}^3$ homogeneous porcine liver block. A 5 mm-diameter laser source touching one surface of this cubic volume irradiated energy for 1, 10, and 20 min with an optical wavelength of 750 nm. The output of the 3D Monte Carlo simulations is a spatial distribution of normalized energy density, which can be used to assist with defining the optical laser delivery as a heat source.

Table 1

Optical and thermodynamic parameters utilized in simulations.

To define a heat source, heat conduction was modeled using COMSOL Multiphysics 6.1.¹⁸ A homogeneous cubic tissue with the same size as the porcine liver block described above (i.e., $20\times 20\times 20{\mathrm{mm}}^3$) was also modeled considering three domains. In the first domain, which constitutes the surface directly exposed to the 5 mm-diameter laser beam, triangular elements with a maximum size of 2 mm were used. In the second domain, which covers the surface not exposed to the laser beam, triangular elements with a maximum size of 3 mm were used. In the third domain, which is the remainder of the modeled volume, quadrilateral elements with a maximum size of 3 mm were used. All boundaries were kept at constant temperature except for the liver surface. Heat exchange in the liver surface was based on free convection in air, resulting from the liver exposure to the environment. This convective effect influences the local temperature dynamics induced by the laser.

To monitor the tissue temperature over time and ultimately predict local tissue damage, the Time-Dependent Bioheat Transfer interface within COMSOL Multiphysics 6.1 was used to solve the Pennes Bioheat Transfer equation. This equation describes the time-dependent biological heat transfer to model hyperthermia processes in perfused tissue, as follows:

The remaining variables are based on properties of the liver (i.e., $\rho$, $T$, $c_p$, $k$), associated blood (i.e., ${\rho }_b$, $T_b$, $c_b$, ${\omega }_b$), and heat sources (i.e., $Q$, $Q_{\mathrm{met}}$). In particular, $\rho$ and ${\rho }_b$ are the liver and blood density, respectively; $T$ and $T_b$ are the temperatures of the liver and blood, respectively; $c_p$ and $c_b$ are the specific heat capacities of the liver and blood, respectively (representing the energies required to raise $T$ and $T_b$, respectively, by one unit of temperature per unit mass); $k$ is the thermal conductivity (which quantifies the ability of the liver to conduct and retain heat); ${\omega }_b$ is the blood perfusion rate (which accounts for the circulation of blood through an in vivo vascularized liver); $Q$ and $Q_{\mathrm{met}}$ are the laser and metabolic heat sources, respectively; and $t$ is the time. The values of ${\omega }_b$, $c_b$, ${\rho }_b$, $T_b$ employed in our study are reported in Table 1 for the porcine liver, whereas $\rho$, $c_p$, and $k$ were considered to be temperature-dependent, based on the details reported by Rossmann and Haemmerich.¹⁹ The impact of body heat production on temperature was not considered due to the lack of sufficient literature values for the porcine liver (i.e., $Q_{\mathrm{met}}=0$), and $Q$ in Eq. (1) was obtained as follows:

A domain point probe was placed at the irradiated surface in the COMSOL simulations to monitor $T(t)$ in Eq. (1), and the associated degree of tissue injury, $\alpha (t)$, was obtained by solving the following differential equation of the Arrhenius Kinetics model:

2.2.

Experimental Procedures to Assess Necrosis

A 5 mm-diameter fiber bundle was connected to a Phocus Mobile laser containing an internal power meter (Opotek, Carlsbad, California) to deliver output desired energies to the surface of exposed in vivo liver tissue. The laser delivered 5 ns pulse widths at a pulse repetition frequency of 10 Hz. To monitor and compensate for known fluctuations in energy and maintain a desired sliding-average pulse-to-pulse energy range (with a 50-pulse window) of 71 to 75 mJ, the internal power meter of the laser and a custom command-line interface were employed.⁵ In particular, a calibration between the internal power meter and an external power meter was performed prior to the start of the procedure, and the internal power meter was used to record and adjust real-time laser energies throughout the duration of the laser application. Figure 1 shows the pulse-to-pulse energy measurements during different laser application times. The interquartile ranges of delivered energy was 70.29 to 76.75 mJ, 70.28 to 73.50 mJ, and 71.37 to 74.59 mJ for the 1, 10, and 20 min time durations, respectively. The mean energy is reported in Table 2.

Fig. 1

Distribution of laser energy delivery per laser application duration, shown as box-and-whisker plots. The median energy is indicated by the red horizontal line, the interquartile range is indicated by the top and bottom bounds of each box, and the maximum and minimum values (excluding outliers, appearing as red datapoints and defined as values $>1.5$ times the interquartile range) are indicated by the lines extending from each box.

Table 2

Laser energy delivery positions and duration.

A laparotomy was performed on the abdomen of a female Yorkshire swine (36 to 40 kg) to access and expose the left lateral liver lobe. Laser energy was applied to three positions on the surface of the left lateral lobe for 1, 10, and 20 min, creating three samples for analysis, as detailed in Table 2. The three differences in irradiation time are expected to induce minimal, moderate, and severe necrosis, based on our previous results.⁵ The irradiated positions were marked by first placing a line of suture perpendicular to and 3 cm away from each intended laser application site. The three suture lines are observable in Fig. 2(a). After laser application, the irradiated regions were directly inked using a tissue marking dye. After euthanasia, the entire left lateral liver lobe was removed from the abdomen, irradiated regions were excised [Fig. 2(b)], and immediately fixed in a 10% formalin solution. This study was approved by the Johns Hopkins University Institutional Animal Care and Use Committee.

Fig. 2

The marking strategies to identify irradiated areas on the swine liver samples included (a) sutures before laser application and (b) tissue dye after liver lobe resection (and prior to sample fixation).

2.3.

Qualitative Histopathology Assessment

To provide baseline comparisons with previous qualitative reports, the excised samples described in Sec. 2.2 were embedded in paraffin, then placed in a microtome with the irradiated surface oriented approximately parallel to the microtome blade. Each sample was sectioned into 250 slides with a section thickness of $4\mu \mathrm{m}$. Out of the collection of 250 sections, one section was selected for staining with Hematoxylin and Eosin (H&E), which was extracted from a depth of $928\mu \mathrm{m}$, to guarantee the representation of the entire cross-section of the tissue sample. As H&E was previously used to qualitatively identify the presence of necrosis, hemorrhage, and inflammation,⁵ the samples described above were similarly assessed for these three pathological features (i.e., absent, minimal, mild, moderate, or severe).

2.4.

Qualitative Immunohistochemistry Assessment

Out of the collection of 250 sections per liver sample described in Sec. 2.3, a representative subset consisting of one per every 10 sections (i.e., $40\mu \mathrm{m}$ spacing) was selected for immunohistochemistry (IHC) staining with an antibody specific to cleaved Caspase-3, diluted at a ratio of 1:1000. In addition, three additional sections near the surface of each sample, separated by $8\mu \mathrm{m}$ per section, were stained with IHC to provide mean and standard measurements for comparison to surface simulation results. Cleaved Caspase-3 was the antibody chosen to identify necrosis biomarkers because of its well-established role in apoptosis, a programmed cell death mechanism. This choice was additionally motivated by our preliminary qualitative observations, which revealed a positive correlation between cleaved Caspase-3 and the identification of necrotic areas resulting from prolonged laser irradiation.

To systematically characterize our qualitative observations, individual IHC sections were digitized at 40× magnification using a Hamamatsu NanoZoomer S210 to generate NDPI files. The NDPITools Plugin Bundle of ImageJ²⁰ was used to extract the content of the NDPI file by dividing each digitized IHC section into a mosaic of adjacent JPEG images, encoding color information using Red-Blue-Green (RGB) channels. NDPITools automatically selected mosaic dimensions that were powers of two in each image dimension, resulting in individual JPEG image sizes ranging 700 to $1300\times 700$ to 1300 pixels, based on a 4MB storage restriction. When these multiple JPEG images were spatially arranged, aligned, and stitched together, the IHC section was accurately reconstructed. As the area of the digitized tissue sections generally increased with depth (see Fig. 3), due to the curved surface of the samples, the number of JPEG images that composed each mosaic varied. For example, the sizes of the IHC sections in Figs. 3(a)–3(c) are $\mathrm{30,720}\times \mathrm{30,976}$, $76,800\times 81,664$, and $\mathrm{103,680}\times \mathrm{101,376}$ pixels, respectively, which were decomposed into $32\times 32$, $128\times 64$, and $128\times 128$ mosaics, respectively, each comprising individual JPEG images sized $960\times 968$, $1200\times 638$, and $810\times 792\text{pixels}$, respectively (i.e., 1024, 8192, and 16,384 individual JPEG images, respectively). Overall, the number of individual JPEG images from the IHC sections associated with the 1, 10, and 20-min laser irradiation samples was 16,384 (constant for each section), 1024 to 16,384, and 2048 to 16,384, respectively.

Fig. 3

Digitized immunohistochemistry sections extracted at depths of (a) $0\mu \mathrm{m}$, (b) $480\mu \mathrm{m}$, and (c) $960\mu \mathrm{m}$ from a porcine liver sample irradiated for 10 min. The tissue section size increases with depth due to the tissue curvature.

To qualitatively identify the presence of biomarkers at the cellular level, we conducted a visual inspection of the digitized IHC sections, revealing two distinct types of cells: (1) cleaved Caspase-3-positive cells stained in intense brown and (2) cleaved Caspase-3-negative cells stained in intense blue, which will be referred to as brown and blue cells, respectively. The distribution of these cells presented four notable characteristics, which are shown in Fig. 4. First, a decrease in cell density was observed inside the irradiated region compared to non-irradiated areas. Second, blue cells (examples denoted with arrows in Fig. 4) were present both inside and outside visibly irradiated areas. Third, brown cells were predominantly observed surrounding the irradiated region. Fourth, the number of brown cells detected within the irradiated area was inversely proportional to the irradiation time.

Fig. 4

Visual inspection of a representative digitized IHC section. Blue cells are pointed out by arrows. The overall cell density was higher outside than inside the irradiated area. The boundary of the irradiated area presented a high concentration of brown cells.

To segment blue and brown cells from the digitized IHC sections, an Attention U-Net²¹ was employed. The architecture of the Attention U-Net consisted of four encoder layers, four decoder layers, and four attention gates. This network was trained using the Adam optimizer, a batch size of four samples, and standard data augmentation techniques (i.e., flip, elastic transformation, grid distortion, and optical distortion). Early stopping with a patience of 10 was employed to avoid overfitting. The Attention U-Net was trained with images from the current dataset. The network dataset contained 50 images along with their manual segmentations, with an 80%–20% training–testing set split, and 20% of the training set forming the validation set. Performance was quantified using the Dice similarity coefficient (DSC) as the main metric, following the guidelines of evaluation metrics for medical image segmentation.²² In addition, the associated intersection over union (IoU), recall, and precision were calculated. The weights obtained from this training process were stored and utilized for cell segmentation in each JPEG image. This segmentation was performed using Python 3.9.15 in a Jupyter Notebook.

Following cell segmentation, individual cells were analyzed for classification into blue or brown. A cell was categorized as brown if its RGB components satisfied the condition $60\le R\le 210$, $G\le 151$, and $B\le 130$. A cell was categorized as blue if its RGB components satisfied the condition $150\le R\le 186$, $G\ge 155$, and $B\ge 160$. These specific thresholds were empirically chosen based on the analysis of multiple images from our dataset. To reconstruct a section segmentation of blue and brown cells, the segmentation masks obtained from individual JPEG images were spatially arranged, aligned, and stitched together to achieve the same size and relative orientations that existed prior to creating the mosaic. This image processing was performed using MATLAB R2023a (Mathworks, Natick, Massachusetts) software.

2.5.

Quantitative Necrosis Mapping

To quantify necrosis as a percentage based on the digitized IHC sections, the observed spatial distributions and qualitative characteristics of blue and brown cells noted in Sec. 2.4 were modeled as an exponential decay, as follows:

2.6.

Registration of Digitized IHC Sections and Necrosis Maps

To place the independent two-dimensional (2D) serial digitized IHC sections and necrosis maps in a common coordinate frame for volumetric reconstruction, a rigid registration was first performed between digitized IHC sections. This registration enabled the computation of deformation maps for each section. The deformation maps were then used to warp the corresponding necrosis maps.

To perform digitized IHC sections registration, a section in the middle of the stack was selected as the reference section based on overall appearance, optimal contrast, and depiction of the entire irradiated area. This reference section selection is critical to avoid error propagation and ensure accurate reconstruction results. We purposely avoided selecting the first section of the stack as the reference, which can introduce reconstruction artifacts such as skewed or helical volumes.²³ The carefully selected reference section was then registered with its two direct neighboring sections, and the registration process was performed in both forward and backward directions (i.e., pairwise registration until the last and first sections, respectively). During this automatic registration process, geometric features were extracted from the IHC sections with the salient feature being the brown cell region surrounding the irradiated area. After registration, the IHC sections were transformed from the RGB to the $L^{*}{a}^{*}b$ color scale to highlight the irradiated area and assist with automatic feature extraction. Each pair of adjacent IHC sections evaluated produced a deformation map, which was then used to warp the corresponding necrosis map.

To remove artifacts from the registered necrosis maps [i.e., necrotic areas appearing in addition to the irradiated region, which are artifacts introduced by Eq. (6)], we implemented a two-stage filtering process consisting of identifying the centroid of necrosis from all maps in a volume, followed by filtering individual maps in the volume, as illustrated in Fig. 5. In Stage 1, the 25 registered necrosis maps were first averaged to yield an average necrosis map, with the maximum value consistently located within the irradiated area. The histogram of the average necrosis map exhibited a quadrimodal distribution, where four peaks represented the background, low, intermediate, and high necrosis percentage values in increasing order of pixel intensity. Using Otsu’s method,²⁴ three threshold values based on pixel intensities of the average necrosis map were computed. As the irradiated region consistently presented regions of intermediate and high necrosis values, the second threshold value was used, converting the average necrosis map into a binary image. Each binary object within the binary image was assigned a non-negative integer as a label. The binary object encompassing the maximum value of the average necrosis map was preserved. The centroid of the selected binary object was computed and utilized to filter individual necrosis maps in Stage 2, which was initiated by converting each individual map to binary objects using the Otsu’s method described in Stage 1. Non-integer numbers were assigned as labels to individual binary objects for easy reference during image processing. Next, the solidity (i.e., a measure of compactness and convexity) of each individual object was computed. Objects with a solidity lower than 0.6 were discarded, as they were unlikely to represent the circular necrosis area. The distances between the centroid of each remaining binary object and the centroid computed in Stage 1 were calculated, and the lowest distance was identified as corresponding to the necrotic area, resulting in a filtered binary mask. The initial individual necrosis map was multiplied by the filtered binary mask to remove artifacts. The 25 filtered binary masks were used to generate a necrosis volume. These masks have a spatial separation of 40 um, resulting from staining one out of every 10 sections using immunohistochemistry (as described in Sec. 2.4). Therefore, linear interpolation was performed to recover the masks from the missing sections and complete the volume. This image processing was performed using MATLAB R2023a (Mathworks, Natick, Massachusetts) software.

Fig. 5

Examples of intermediate and final outcomes of the two-stage artifact filtering process.

2.7.

Comparisons between Simulated and Experimental Results

To compare simulation predictions and experimental results, we first consider that a direct comparison is challenging due to the liver samples exhibiting a curved surface, whereas the simulated tissue model was designed with a flat surface. To address this geometric disparity, we selected the most superficial IHC section that fully displayed the irradiated region and its adjacent non-irradiated area. Following this protocol, the sections selected for comparison were extracted at depths of 0, 280, and $160\mu \mathrm{m}$ for the 1, 10, and 20-min irradiation times, respectively. We then calculated the absolute error in percent necrotic tissue between simulated predictions (Sec. 2.1) and experimental results for the selected sections (Sec. 2.5), as a function of irradiation time. To strengthen the validity of our findings, five total sections at depths 280 to $320\mu \mathrm{m}$ and 160 to $200\mu \mathrm{m}$ from the 10- and 20-min samples, respectively, were each combined to report a mean and standard deviation per sample for additional comparison to the corresponding simulation results.

3.1.

Simulated Necrosis Predictions

Figure 6 shows simulated temperature as a function of time for the delivery of 30 and 73 mJ laser energies. The initial temperature of each sample was 293.15 K (because laparotomy typically exposes the liver to the environment). The temperature profile for the 73 mJ laser energy exhibited three phases. During the first phase (0 to 0.34 min), there was a rapid linear temperature increase in the liver, up to 315.15 K where the inactivation of vital enzymes occurs.²⁵ In the second phase (0.34 to 3.88 min), the temperature continued escalating at a lower rate, eventually reaching a maximum temperature of 333.15 K, which is the threshold temperature for protein denaturation.²⁵ In the third phase, the temperature stabilized and was maintained at 333.44 K. In contrast, the temperature profile associated with the 30 mJ laser energy has two phases, with the first stage (0 to 4.94 min) characterized by a gradual temperature increase (up to 319.56 K, which surpasses the temperature known to inactivate vital enzymes²⁵) and the second phase characterized by temperature stabilizing at 319.66 K.

Fig. 6

Time-resolved temperature prediction at the tissue-laser interface for a simulated liver exposed to 73 and 30 mJ laser energy.

Figure 7 shows the time-resolved predictions for necrosis percentage at the tissue-laser interface for different exposure times and laser energies. The induced damage slowly progressed during the first minute of irradiation, resulting in necrotic tissue percentages of 0.03% and 1.23% for 30 and 73 mJ, respectively. Extending the laser exposure to 10 min led to necrotic tissue percentages of 7.49% and 66.23% for 30 and 73 mJ, respectively. Following 20 min of laser irradiation, the necrotic tissue percentage escalated to 15.05% and 76.84% for 30 and 73 mJ, respectively. Table 3 summarizes these simulation results.

Fig. 7

Time-resolved necrosis percentage prediction at the tissue-laser interface for simulated liver samples irradiated for (a) 1 min, (b) 10 min, and (c) 20 min.

Table 3

Comparison of simulated predictions and qualitative histopathological grading of necrosis in swine liver samples irradiated with ∼30 and 73 mJ laser energy for 1, 10, and 20 min.

3.2.

Empirical Necrosis Assessment

Figure 8 displays the characterization results of the H&E-stained sections. These sections were extracted at depths of $4\mu \mathrm{m}$ below the most superficial section that fully displayed the irradiated region and its adjacent non-irradiated area (i.e., 4, 284, and $164\mu \mathrm{m}$ for the 1, 10, and 20-min irradiation times, respectively) to closely match the quantitative assessment depths described in Sec. 2.7. Necrosis was defined by the loss of both hepatocellular and structural integrity, characterized by the loss of cell nuclei, cellular fragmentation, and content leakage. Inflammation manifested as aggregates of leukocytes within tissue. Hemorrhage was defined as the extravasation of red blood cells into interstitial spaces between cells. Based on these criteria, the porcine liver sample irradiated for 20 min exhibited severe necrosis and severe hemorrhage, whereas the sample irradiated for 10 min displayed mild necrosis and severe hemorrhage. In contrast, the porcine liver sample exposed to 1-min irradiation showed no pathological conditions. Minimal to no inflammation was observed in the three samples, as no leukocyte clusters were seen.

Fig. 8

H&E sections from porcine liver samples irradiated for 1, 10, and 20 min. Regions of interest within the irradiated region were selected to perform necrosis and hemorrhage grading.

Figure 9 shows example results from the digitized IHC section processing pipeline. A representative IHC section from a swine liver sample irradiated for 10 min [Fig. 9(a)] was digitized and subsequently divided into a mosaic composed of 16384 JPEG images [Fig. 9(b)], and individual JPEG images were segmented [Fig. 9(c)]. These JPEG images were segmented with acceptable performance by the Attention U-Net (i.e., 0.97 DSC, 0.94 IoU, 0.99 recall, and 0.95 precision). The arrangement of the segmentation masks [Fig. 9(d)] provided a clearer image of the spatial distribution of blue and brown cells within the section. The overall cell density inside the irradiated area was notably lower compared to external areas. In addition, blue and brown cells were present throughout the section, whereas the presence of brown cells was reduced within the irradiated area. This segmentation example confirms the features observed during visual inspection of the IHC sections. The local necrosis map [Fig. 9(e)] computed using Eq. (6) exhibited a rounded central region with the highest necrosis percentage (i.e., mean percentage of 48.45%), surrounded by two scattered regions of lower and variable necrosis percentage (e.g., 28.11% and 35.44% mean percentage per independent region, resulting in an unrealistic combined necrosis percentage exceeding 100%, which leads to their classification as artifacts). The necrosis values corresponding to the irradiated region were preserved after filtering the artifacts [Fig. 9(f)].

Fig. 9

Image processing workflow for necrosis quantification of IHC sections: (a) IHC section, (b) mosaic of 16,384 RGB images, (c) segmentation, (d) blue and brown cells segmentation, (e) necrosis map, and (f) filtered map.

Figure 10 shows the necrosis maps from the irradiated liver samples, each extracted from a depth of $440\mu \mathrm{m}$ below the tissue surface. Laser application for 20 min caused complete disruption of cells in the illuminated area, yielding a mean necrosis percentage of 85.12%. Laser application for 10 min caused moderate disruption of cells in the illuminated area, resulting in a mean necrosis percentage of 65.49%. There are no signs of cell disruption when the laser is applied for 1 min, yielding 0% necrosis percentage in the irradiated region.

Fig. 10

IHC section, blue and brown cells segmentation, initial necrosis map, and filtered necrosis map for three swine liver samples.

Figure 11 shows the results obtained when quantifying the volume and area of necrotic cells as a function of depth from the irradiated surface. In Figs. 11(a)–11(c), the necrosis volumes measured 0, 3.27, and $24.68{\mathrm{mm}}^3$ after 1, 10, and 20 min of irradiation, respectively. In Fig. 11(d), the cross-sectional necrotic area remained generally constant and large (0.21 to $0.24{\mathrm{cm}}^2$) within the analyzed depth range (1 mm of swine liver tissue) for a 20-min exposure time. When the exposure time was shortened to 10 min, the cross-sectional necrotic area decreased at shallow depths (i.e., $0.03{\mathrm{cm}}^2$), then further decreased to $0.01{\mathrm{cm}}^2$ at a depth of $960\mu \mathrm{m}$. No irradiation effects were observed with the 1-min laser irradiation time.

Fig. 11

Volumetric necrosis reconstruction for swine liver samples irradiated for (a) 1 min, (b) 10 min, (c) 20 min, and (d) cross-sectional necrosis area as a function of tissue depth and laser duration.

Figure 12 shows necrosis percentage as a function of depth for varying experimental laser exposure times. With a 20-min exposure time, the measured necrosis percentage ranged from 77.69% to 87.06%. When the exposure time was decreased to 10 min, the necrosis percentage ranged 38.95% to 66.22% for a similar mean energy of 73 mJ (see Fig. 1 and Table 2). There was no necrosis for the 1 min exposure time. When compared to the simulated results in Table 3 (obtained at a depth of $0\mu \mathrm{m}$), the experimental results obtained at depths of 0, 280, and $160\mu \mathrm{m}$ for 1, 10, and 20-min irradiation times, respectively (due to the curvature of the experimental tissue samples, as noted in Sec. 2.7) revealed 0%, 66.22%, and 84.94% tissue necrosis, respectively, whereas the simulation framework predicted 1.23%, 66.23%, and 76.84%, respectively, which corresponds to deviations of 1.23%, 0.01%, and 8.1%, respectively. Therefore, the overall experimental deviation from predicted values range 0.01% to 8.1%. The experimental results obtained at depths of 280 to $320\mu \mathrm{m}$ and 160 to $200\mu \mathrm{m}$ for 10- and 20-min irradiation times, respectively, yielded mean ± one standard deviation tissue necrosis percentages of $61.24\pm 4.49\%$ and $86.38\pm 1.15\%$, respectively. Therefore, the 0.01% to 8.1% deviations predicted with the simulation framework are within 1 to 2 standard deviations of the variations obtained over a minimal depth of $40\mu \mathrm{m}$.

Fig. 12

Necrosis percentage as a function of tissue depth and laser duration.

Overall, our empirical results demonstrate that liver imaging with no necrosis can be achieved with 73 mJ laser energy ($371.79\mathrm{mJ}/{\mathrm{cm}}^2$ fluence) applied for 1 min, with a laser wavelength of 750 nm, 5 ns pulse duration, and 10 Hz pulse repetition frequency. However, measurable necrosis was observed for exposure durations of 10 and 20 min, which is not considered safe under the same laser conditions (73 mJ energy, 750 nm, 5 ns pulse duration, and 10 Hz pulse repetition frequency). The 0.01% to 8.1% agreement between quantitative experimental results and simulation results, combined with $\le 4.49\%$ standard deviations on the quantitative measurements and the consistency between quantitative experimental results and qualitative outcomes from H&E staining collectively demonstrate the accuracy of our approach.