2.1.

Optical Transport Model

Our simulation of light transport in skin utilizes the well-established weighted-photon MC technique.¹⁴ In this approach, a large number of energy packets (“photons”) propagate through the tissue and deposit fractions of their energy into specific volume elements, according to local absorption properties. Probabilistic relations are used to compute each photon trajectory, based on the macroscopic laws of light scattering, reflection, and refraction. These relations involve optical properties such as absorption coefficient (${\mu }_a$), scattering coefficient (${\mu }_s$), anisotropy coefficient ($g$), and index of refraction ($n$) of the involved tissues. With an increasing number of simulated photons and ever finer spatial discretization of the treated volume, the resulting energy deposition maps converge toward a realistic spatial distribution of light absorption in the modeled organ.

Our implementation of the MC code is an extrapolation of the popular multilayer MC model (MCML)¹⁴ to three spatial dimensions. This allows treatment of skin inclusions of arbitrary shape, position, and orientation, as well as irradiation with spatially nonuniform light beams.¹ We have implemented two versions of the 3-D MC code, differing only in the treatment of curved boundaries between neighboring tissues.

In the customary version of the code, tissue geometry is represented within the confines of the Cartesian discretization grid. Each voxel is assigned to a single tissue and thus obtains a uniform set of optical parameters (${\mu }_a$, ${\mu }_s$, $g$, and $n$). Consequently, boundaries of anatomical structures are approximated by flat or serrated surfaces pinned to the corners of the discretization grid. In the following, we refer to this approach as “model D.”

In order to avoid the potential deficiencies of the above approach, we have implemented an augmented MC model which utilizes analytically defined tissue borders (“model G”), e.g., a blood vessel is represented by a cylinder with a perfectly smooth surface (Fig. 1). The cuboid voxel structure is retained to store the 3-D energy deposition data for subsequent rendering and analysis. However, individual voxels do not necessarily contain only one tissue type. Voxels which include a boundary between neighboring tissues thus accumulate the laser energy absorbed on both sides of the boundary.

Fig. 1

Our main model geometry includes a superficial epidermal layer and semi-infinite dermis with a horizontal blood vessel (red). The green circle marks the perimeter of the incident laser beam, and black cube indicates the computational volume of interest (VOI).

In terms of basic light transport, our model G proceeds along the same MC simulation steps as model D (i.e., photon generation, relocation, absorption, and scattering). The only difference between the two versions is in the treatment of light reflection and refraction at curved boundaries between different tissues, as illustrated in Fig. 2.

Fig. 2

Reflection and refraction of light at the vessel boundary in MC model with (a) tissue boundaries pinned to the Cartesian grid and (b) analytically defined tissue boundaries.

In short, when a propagating photon packet crosses a tissue boundary, the normal unity vector at the intersection point, $\hat{n}$, is calculated. The cosine of the incidence angle ($\alpha$) is then computed as the dot product between $\hat{n}$ and the current photon direction, $\mu$. If $\alpha$ is smaller than the critical angle for total internal reflection, ${\alpha }_c=\mathrm{arc}\sin (n_2/n_1)$, and a random-number generator is applied to determine whether the photon will be reflected or refracted based on the reflectivity value calculated using the Fresnel’s equation. Detailed expressions and their derivations are presented in the Appendix.

In both implementations (models D and G), special attention is devoted to treatment of the VOI side boundaries. As pointed out before, neither an “escape” nor “mirror” (insulating) boundary condition is appropriate for high-resolution treatment of small volumes in human skin irradiated with significantly larger yet spatially confined laser beams.⁶ One prominent example is the cutaneous vascular lesions called port-wine-stains (PWS). These lesions contain blood vessels with diameters of 20 to $300\mu \mathrm{m}$ and are commonly treated with lasers beams of several millimeters in diameter. Under such conditions, the most common solution—i.e., to control edge-related artifacts by moving the “escape” or “mirror” VOI boundaries sufficiently far from all relevant structures—becomes prohibitively memory demanding and time consuming.⁶

Our approach to the above problem, called “extended layers,” is to also propagate the photons outside of the VOI in laterally infinite tissue layers, until they lose most of their energy or escape into the air. Optical properties of these layers are inherited from the border VOI elements at the respective depths, and the majority rules when they do not represent the same tissue type. Most importantly, the energy deposition outside of the VOI is not recorded. This keeps the memory requirements and/or the number of hard drive read/write operations at a manageable level and thus reduces the computation time by several orders of magnitude.⁶ The same principle is also applied at the bottom boundary of the VOI. Specifically, the entire semispace below the VOI bottom boundary inherits the optical properties of the dermis.

Both versions of the MC code were implemented in Visual C++ 2010 (Microsoft, Redmond, Washington) with utilization of multithreading to help speed up the computations.

2.2.

Model Geometry and Optical Parameters

Our main test geometry consists of two laterally infinite layers, representing the superficial epidermal layer (thickness $60\mu \mathrm{m}$) and semi-infinite dermis underneath. Energy deposition is recorded within the VOI with dimensions $x=y=z=0.6\mathrm{mm}$. A single horizontal blood vessel has its axis located at depth $z=195\mu \mathrm{m}$ and laterally at $x=0.3\mathrm{mm}$ (see Fig. 1). Vessel diameters of 30 and $150\mu \mathrm{m}$ will be considered, representing small and large PWS vessels, respectively.

The photon-launching pattern imitates a perpendicularly incident laser beam with a top-hat intensity profile, centered above the blood vessel. Irradiation with a narrow (diameter $d_b=1.2\mathrm{mm}$) and a wider laser beam ($d_b=5.0\mathrm{mm}$) will be discussed.

All tissue optical parameters correspond to the wavelength of 532 nm (Table 1).^6,7 Melanin is assumed to be homogenously distributed within the epidermis. The blood absorption coefficient corresponds to a hematocrit of 40% and an oxygen saturation level of 70%.

Table 1

Absorption coefficient (μa), scattering coefficient (μs), anisotropy factor (g), and refractive index (n) values used in our Monte Carlo model of optical transport (λ=532  nm). The μa value for blood corresponds to hematocrit Ht=40% and oxygenation level of 70%, and the scattering properties were determined for Ht=33%.

3.1.

Laser Energy Deposition in a 30$\text{-}\mathit{\mu }\mathrm{m}$ Vessel

Figure 3 presents the spatial distribution of the energy deposited in a $30\text{-}\mu \mathrm{m}$ vessel geometry upon irradiation with a 1.2-mm diameter laser beam, as obtained using the MC model D. The indicated energy deposition values, $a$, correspond to a laser radiant exposure of $10\mathrm{J}/{\mathrm{cm}}^2$. The number of incident “photons” was $1.5\times {10}^8$ and the spatial discretization step in the cross-sectional plane ($x,z$) equals $3\times 3{\mu \mathrm{m}}^2$.

Fig. 3

Cross-sectional maps of energy deposition in our model geometry with a $30\text{-}\mu \mathrm{m}$ blood vessel irradiated with a 1.2-mm diameter hat-top laser beam: (a) in the central plane perpendicular to the vessel axis ($y=0.3\mathrm{mm}$); and (b) in the vertical plane containing the vessel axis ($x=0.3\mathrm{mm}$).

As expected, the energy deposition values are largest within the blood vessel, significantly lower in the epidermis, and minimal in the dermis. Energy deposition within the vessel lumen is nonuniform, with the largest values indicated along its top boundary [Fig. 3(a)]. An energy deposition map in the vertical plane containing the vessel axis [Fig. 3(b)] demonstrates the absence of considerable edge artifacts.

In Fig. 4, we compare the spatial distributions of laser energy deposition within the blood vessel lumen and nearby dermis as obtained using our two simulation models. The differences are evident, especially along the vessel border. In the results of model D [Fig. 4(a)], the energy deposition region is sharply demarcated by the serrated vessel boundary. In Fig. 4(b), in contrast, the transition between the vessel and nearby dermis is blurred due to specific treatment of the tissue boundaries in model G. Here, the border voxels accumulate energy deposited in both blood and dermis on either side of the vessel boundary.

Fig. 4

Energy distribution within the vessel of $30\text{-}\mu \mathrm{m}$ diameter and in the nearby dermis obtained by the MC models D (a) and G (b).

For quantitative analysis of the related effects, we assess the relative differences between the energy depositions as predicted by models D and G for each VOI voxel (marked with matrix indices $i,j$, and $k$):

Because the test geometry is uniform along the $y$-axis, we control the inevitable noise in the MC results by using significantly coarser discretization in this direction [$\mathrm{\Delta }y=20\mu \mathrm{m}$; see Fig. 3(b)]. In addition, we can average the results over a certain number of ($x,z$) slices near the VOI center, where the light fluence is practically constant along the $y$ direction. This leads to the noise-reduced version of Eq. (1):

Based on the analysis of the data in Fig. 3(b), in the following we apply averaging over the central 16 slices. The remaining part of the VOI (i.e., the last 0.14 mm on either side) features a small decrease in energy deposition, resulting from the relatively narrow laser beam used in this example (see Fig. 1).

The obtained map of ${\langle \delta \rangle }_{ik}$ values is presented in Fig. 5. In accordance with the above discussion of Fig. 4, a systematic and thus rather symmetrical pattern is evident along the vessel border. Individual values of ${\langle \delta \rangle }_{ik}$ are either positive (yellow to red colors) or negative (blue), depending on whether the respective voxel is assigned to the inside or outside of the vessel in model D. In some of these voxels, the amplitudes approach 100%.

Fig. 5

Relative difference of the energy deposition values in MC models D and G, ${\langle \delta \rangle }_{i,k}$ [Eq. (2)] for the $30\text{-}\mu \mathrm{m}$ diameter vessel. Red circle indicates the vessel boundary.

However, we are interested in systematic differences between the two implementation models beyond such inevitable effects in the border voxels. We can expect the latter to largely average out and lose relative importance when applying sufficiently fine spatial discretization.

To that end, let us consider the vertical profile of ${\langle \delta \rangle }_{ik}$ through the vessel center, as presented in Fig. 6 (solid line). The value is positive throughout the vessel lumen, increasing from $\sim 2\%$ at its top to 6% at the bottom. In contrast, ${\langle \delta \rangle }_{ik}$ in the epidermis is mostly negative and amounts to $-2\%$. The second effect is in good agreement with the first. Namely, as more energy is absorbed in the vessel in model D as compared with model G, fewer photons are scattered back toward the epidermis, thus reducing the fluence and energy deposition in the latter.

Fig. 6

Vertical profile of ${\langle \delta \rangle }_{i,k}$ [Eq. (2)] through the center of the $30\text{-}\mu \mathrm{m}$ diameter vessel ($x=y=0.3\mathrm{mm}$; dark blue line). Standard deviation and standard error of the mean are indicated with thin lines (red) and shaded band (light blue), respectively. Vertical lines indicate the boundaries of the epidermis and blood vessel.

In order to test the significance of these effects, we compute the standard deviation of ${\delta }_{ijk}$ values at each location ($i,k$):

Since numerical noise in different voxels is uncorrelated, we can also assess the standard error of ${\langle \delta \rangle }_{ik}$

For further analysis, we consider the energy uptake by the blood vessel per unit of length, $a_{\text{L}}$. For a selected cross-sectional plane (i.e., slice $j$), this quantity can be computed by summing up the energy deposition contributions, $a_{ijk}$, across the vessel lumen:

In Fig. 7, we compare such linear energy deposition values, $\langle a_{\text{L}}\rangle$, as obtained using our two MC models with the lateral discretization step ($\mathrm{\Delta }x=\mathrm{\Delta }z$) varied between 1 and $15\mu \mathrm{m}$. The values obtained with model D (open circles) evidently feature the significant and rather erratic variations with the discretization step. In contrast, no such unphysical behavior is observed in the results of our augmented model G (solid squares).

Fig. 7

Linear energy deposition, $\langle a_{\text{L}}\rangle$ [Eq. (5)] as a function of the spatial discretization step ($\mathrm{\Delta }x$) for model D (open circles) and model G (solid squares).

Our analysis shows that model D mostly overestimates absorption of laser light by the blood vessel. However, the opposite effect can be seen at $\mathrm{\Delta }x=12\mu \mathrm{m}$. The relative difference of linear energy depositions between the two models,

At finer spatial discretization, $\langle {\delta }_{\text{L}}\rangle$ assumes exclusively positive and somewhat smaller values, below 8% at $\mathrm{\Delta }x\le 4\mu \mathrm{m}$. However, the values do not show any tendency to converge toward zero, but rather stabilize at 4% to 5%. This points to a fundamental difference in interaction of incident light with the model vessel, as predicted by our two MC models.

3.1.1.

Shifting of discretization grid

In model D, the approximated shape of the model vessel depends not only on the discretization step, but also on location of the origin of the computational grid. In Fig. 8, we illustrate the latter effect for the $30\text{-}\mu \mathrm{m}$ diameter vessel (at $\mathrm{\Delta }x=\mathrm{\Delta }z=6\mu \mathrm{m}$). In Fig. 8(a), the vessel axis coincides with the center of the border between two neighboring voxels. This was the case in the example discussed thus far (Figs. 3–7) and is referred to as the “original grid” in the following. In Fig. 8(b), the discretization grid is shifted in the $x$ direction by $\mathrm{\Delta }x/2$, such that the vessel axis passes through the voxel center (“shifted grid”). As is evident from Fig. 8, this results in a considerable change of the vessel’s cross section in model D.

Fig. 8

Transversal cross section of the energy deposition in the $30\text{-}\mu \mathrm{m}$ diameter vessel centered (a) at the boundary between two voxels, and (b) at the voxel center. Result of model D.

When using the augmented model G, we do not see any influence of the grid positioning (or $\mathrm{\Delta }x$) on the linear density of energy absorbed in the vessel, $a_{\text{L}}$ [Eq. (5)]. Therefore, we compare in Fig. 9(a) only the relative deviations $\langle {\delta }_{\text{L}}\rangle$ [Eq. (6)] as assessed for both grid positions and varying discretization steps.

Fig. 9

(a) Relative difference of the linear energy depositions, $\langle a_{\text{L}}\rangle$, and (b) relative error of the cross-sectional area in model D ($\delta A$) as a function of discretization step ($\mathrm{\Delta }x$) for the original (open circles) and shifted computational grids (solid, red).

For most small discretization steps (i.e., $\mathrm{\Delta }x\le 10\mu \mathrm{m}$), considerably smaller deviations are indicated when using the shifted grid (solid circles; red) as compared with the original one (open circles). Similarly, $\langle {\delta }_{\text{L}}\rangle$ assumes significantly lower values (below 2%) than with the original grid at the finest discretizations under test. This likely reflects the fact that with the shifted grid, the vessel cross section in model D always features a fourfold rotational symmetry and thus more closely resembles the actual shape of the vessel lumen (see Fig. 8).

At the larger discretization steps ($\mathrm{\Delta }x\ge 12\mu \mathrm{m}$), the shifted grid yields somewhat larger mismatches than the original grid. We consider this effect to be incidental and attribute it to the increasingly deteriorating representation of the vessel cross section at larger $\mathrm{\Delta }x$, which dominates over the influence of grid positioning.

For the present analysis, it is also important to note that the vessel’s cross-sectional area in model D ($A^D$) can differ from its actual value ($A=\pi d^2/4$) and can also vary with both $\mathrm{\Delta }x$ and grid positioning. In Fig. 9(b), we plot the relative mismatch, $\delta A=(A^D-A)/A$, for the involved discretization steps and both grid positions. A close similarity between the behavior of $\delta A$ and that of $\langle {\delta }_{\text{L}}\rangle$ [see Fig. 9(a)] is evident for the original grid (open circles). For the shifted grid (solid circles), the similarity is less apparent and limited to smaller discretization steps ($\mathrm{\Delta }x<5\mu \mathrm{m}$).

3.2.

Energy Deposition in a 150-μm Vessel

Our next example involves a larger blood vessel (diameter $d=150\mu \mathrm{m}$), irradiated with the same laser beam as above ($d_b=1.2\mathrm{mm}$).

A cross-sectional map of the relative differences between the energy depositions obtained using the two MC models, ${\langle \delta \rangle }_{ik}$ [Eq. (2)], is presented in Fig. 10. Apart from the trivial discrepancies evident in the border voxels, energy deposition in most of the vessel lumen is indicated to be somewhat larger in model D as compared with model G. Similar to the $30\text{-}\mu \mathrm{m}$ vessel (Figs. 5 and 6), the largest deviations are observed in the bottom part of the vessel lumen. The values of $\langle \delta \rangle$ here approach 10%, which is considerably more than in the previous example. This is somewhat counter-intuitive because, at the same discretization step of $3\mu \mathrm{m}$, the boundary of the $150\text{-}\mu \mathrm{m}$ vessel appears much better approximated than that of the five times smaller vessel [see Fig. 4(a)].

Fig. 10

Transversal cross section of the skin model involving the $150\text{-}\mu \mathrm{m}$ vessel, indicating the average relative difference of energy deposition between the MC models D and G, ${\langle \delta \rangle }_{ik}$.

A large area of positive $\langle \delta \rangle$ values is also evident in the dermis below the vessel. The discrepancies here amount to 8% to 18%, which is more than that within the vessel. We also find it interesting that the affected region does not shrink with subsurface depth. Moreover, as evidenced by the vertical profile in Fig. 11, the $\langle \delta \rangle$ values there do not decrease with depth but rather tend to increase toward the bottom of the VOI.

Fig. 11

Average relative difference of deposited energy profile through the axis of a vessel with a diameter of $150\mu \mathrm{m}$ (blue line). Standard deviation and standard error of the average value are presented with thin lines (red) and shaded band (light blue), respectively. Boundaries of the epidermis and vessel are indicated by vertical lines.

In contrast, small negative values of $\langle \delta \rangle$ are apparent throughout the epidermis and in the dermis laterally to the vessel border. As will be discussed further below, these discrepancies are mutually consistent and can be linked to the inaccurate treatment of light reflection and refraction at the vessel boundary in MC model D.

The influences of discretization step and grid positioning for this test geometry are presented in Fig. 12. In the original grid (open circles), the values of $\langle {\delta }_{\text{L}}\rangle$ are nearly independent of $\mathrm{\Delta }x$ and stabilize around 6% at sufficiently fine discretization ($\mathrm{\Delta }x<5\mu \mathrm{m}$). With the shifted grid, in contrast, $\langle {\delta }_{\mathrm{L}}\rangle$ varies much more (solid circles; red) and increases nearly monotonically from approximately –5% at $\mathrm{\Delta }x\ge 12\mu \mathrm{m}$ to $+5\%$ at the smallest discretization step under test.

Fig. 12

(a) Relative difference of the linear energy depositions, $\langle {\delta }_{\text{L}}\rangle$, and (b) relative error of surface area ($\delta A$) as a function of discretization step ($\mathrm{\Delta }x$) for the $150\text{-}\mu \mathrm{m}$ diameter vessel in the original (crosses) and shifted grids (squares).

Note that the shape of the $150\text{-}\mu \mathrm{m}$ diameter vessel in model D is represented rather accurately at the involved discretization steps. For both grid positions under test, the relative error of its cross-sectional area [$\delta A$; Fig. 12(b)] practically vanishes for $\mathrm{\Delta }x<3\mu \mathrm{m}$ and thus has a negligible influence on the observed behavior of $\langle {\delta }_{\text{L}}\rangle$ at the smallest $\mathrm{\Delta }x$. This strongly suggests that the latter results from an inherent deficiency in the treatment of light reflection and refraction in model D.

3.3.

Laser Beam Diameter of 5 mm

The same analysis as presented above is performed next for irradiation with a 5.0-mm diameter hat-top laser beam.

For these examples, we apply the “mirror” boundary condition at the VOI side boundaries. This is based on our earlier observation that the effect of a 5-mm laser beam on shallow blood vessels practically matches that of an infinitely wide incident wave.⁶ Because the photons are now launched only within the VOI top surface, a smaller number of them (i.e., $5\times {10}^7$) provide the same incident photon flux and thus a comparable level of numerical noise to that seen in the former results. This approach thus speeds up the modeling process considerably while yielding sufficiently accurate results, as was verified by direct comparisons with the more versatile extended-layer approach.

For the $30\text{-}\mu \mathrm{m}$ blood vessel, the energy uptake per unit of length, $\langle a_{\text{L}}\rangle$, is 2.9 times higher than with the 1.2-mm beam if the same radiant exposure of $10\mathrm{J}/{\mathrm{cm}}^2$ is considered (see Table 2). However, the differences between the results from our two MC models are qualitatively and quantitatively very similar to those obtained with the 1.2-mm laser beam (Figs. 4–9) and are, therefore, not presented. For example, the relative difference between the linear energy deposition values at the finest discretization equals $\langle {\delta }_{\text{L}}\rangle =2\%$, which is the same as that obtained for the narrow laser beam [see Fig. 9(a)].

Table 2

Linear density of energy deposition, 〈aL〉 (in J/mm), in blood vessels with diameters of 30 and 150  μm irradiated with hat-top laser beams with diameters of 1.2 and 5.0 mm. Results from MC model G for radiant exposure of 10  J/cm2.

For the larger blood vessel ($d=150\mu \mathrm{m}$), the increase in $\langle a_{\text{L}}\rangle$ associated with switching to the larger beam diameter is 2.5-fold (Table 2). In this case, the MC model D predicts a somewhat higher energy deposition as compared with model G in most of the blood vessel (Fig. 13). In contrast to the result obtained for the 1.2-mm laser beam, however, the effect is now most prominent in the lateral parts of the vessel lumen (yellow/orange hues), while no significant difference is indicated near the vertical symmetry plane. The $\langle {\delta }_{\text{L}}\rangle$ values do not exceed 2%, which is considerably less than in the former example.

Fig. 13

Transversal cross section of the skin model involving a $150\text{-}\mu \mathrm{m}$ vessel irradiated by a 5.0-mm diameter beam, showing average relative difference of the energy deposition between the models D and G.

Negative values of $\langle \delta \rangle$ (blue) reaching up to $-6\%$ can be observed in the dermis lateral to the vessel boundary. The underestimation of energy deposition in this area by model D is significantly more pronounced as compared to irradiation with the narrower laser beam (see Fig. 10). Another evident difference with the former result is the absence of positive values underneath the blood vessel. In the epidermis, $\langle \delta \rangle$ amounts to approximately $-1\%$.

In Fig. 14, we compare the linear energy deposition values, $\langle a_{\mathrm{L}}\rangle$, obtained with the two MC models and by varying the discretization step. Just as seen above for the narrow-beam irradiation (Fig. 7), the values from our augmented model G (solid squares) are practically constant. Meanwhile, the results from model D (open circles) feature an apparently erratic variation with $\mathrm{\Delta }x$, which turns out to be closely correlated with the error of the vessel’s cross-sectional area, $\delta A$ [see Fig. 12(b)].

Fig. 14

Linear energy deposition, $\langle {\delta }_{\text{L}}\rangle$ [Eq. (5)] in the $150\text{-}\mu \mathrm{m}$ diameter vessel irradiated with 5-mm laser beam, as a function of the spatial discretization step ($\mathrm{\Delta }x$) for models D (open circles) and G (solid squares).

The relative difference between the two linear energy deposition values, $\langle {\delta }_L\rangle$ [Eq. (6)], is in the range of 0% to 2% and stabilizes around $+1.2\%$ with the finest discretizations. This value is markedly smaller than that obtained for the narrow laser beam [Fig. 12(a)]. Nevertheless, it is undoubtedly different from zero and positive, which is the same qualitative effect as seen in all the above examples.

Shifting of the computational grid by $\mathrm{\Delta }x/2$ did not have any appreciable effect on the predictions of integral physical quantities such as $\langle a_{\text{L}}\rangle$ at sufficiently fine spatial discretization and is, therefore, not documented.

3.4.

Undulating Epidermal–Dermal Boundary

In the test geometries analyzed above, the boundary between the epidermis and dermis was modeled as a flat surface (Fig. 1). In this way, we could isolate the effects resulting from different treatments of light reflection and refraction at the curved boundary between the blood vessel and dermis in two implementations of the 3-D MC modeling approach (i.e., models D and G). Interaction of model photons with the flat epidermal–dermal (ED) boundary was identical (and physically correct) in both models.

In reality, the ED boundary in human skin is strongly undulating. In addition, the respective refractive indices differ even more than those of blood and dermis, $n_{\mathrm{epi}}/n_{\mathrm{der}}=1.06$ versus $n_{\mathrm{der}}/n_{\text{bl}}=1.03$ (see Table 1). It is, therefore, of interest to check to what extent the above simplification might have affected the modeling results, and in particular our provisional conclusions regarding the deficiency of the approach implemented in model D.

To this effect, we model the ED boundary with a periodically undulating function:

For its parameters, we choose the mean epidermal thickness ($z_{\mathrm{epi}}=60\mu \mathrm{m}$), undulation amplitude ($z_0=30\mu \mathrm{m}$), and characteristic lateral dimension ($l_x=l_y=150\mu \mathrm{m}$) typical for human skin.

In MC model D, the ED boundary is approximated by a serrated surface, similar to the boundary between the vessel and dermis. Apart from this change of sample geometry, the modeling process follows identical steps as in previous examples.

In the augmented model G, however, each photon relocation now requires additional testing of whether it has crossed the ED boundary defined by Eq. (7). If so, the first crossing point (in the direction of photon travel) is determined by solving an implicit equation. After computing the unity normal vector ($\hat{n}$) at the intersection point, modeling proceeds in the same manner as at the vessel–dermis boundary (see Appendix).

We will present here only one example, the physically and clinically most relevant case involving the $150\text{-}\mu \mathrm{m}$ blood vessel irradiated with the wide laser beam ($d_b=5\mathrm{mm}$). Because the test geometry does not possess the same symmetry as before, the “mirror” boundary condition at the side VOI boundaries is no longer appropriate. However, the commensurability between the periodicity of the ED boundary ($150\mu \mathrm{m}$) and the VOI dimension ($600\mu \mathrm{m}$) allows us to apply the so-called periodic boundary condition. Whenever a photon crosses the VOI side boundary, it is translated to the opposite side of the VOI with the same direction vector, $\mu$.

3.4.1.

Differences between the MC Models D and G

Figure 15 presents the relative differences between the energy deposition values as obtained using our two MC models, ${\langle \delta \rangle }_{ijk}$ [Eq. (1)], in cross-sectional planes perpendicular to the vessel axis (left panel) and along the same (right). The results were computed with spatial discretization of $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z=3\mu \mathrm{m}$ and $1.5\times {10}^8$ incident photons.

Fig. 15

Relative difference of the linear energy deposition as a function of discretization step size for the $150\text{-}\mu \mathrm{m}$ diameter vessel and 5.0-mm diameter laser beam.

The relatively high noise level in the image reflects the fact that, due to lower symmetry of the model geometry, we could not average the results from multiple central slices ($x,z$). Nevertheless, the periodicity of $z_{\mathrm{ED}}$ [Eq. (7)] allows us to average the data from four identical blocks of thickness $l_y=150\mu \mathrm{m}$. Consequently, Fig. 15 presents the values along the $y$ direction for only one such segment.

The result clearly indicates that the energy deposition in MC model D exceeds that predicted by the augmented model G in the dermis under the “peaks” of the ED boundary. At the same time, it underestimates energy deposition in epidermal “troughs” (pockets) and the dermis just below them. The relative deviations in these areas reach values of $+6\%$ and $-6\%$, respectively.

The differences within the dermis decrease with depth due to the increasing influence of light scattering. Nevertheless, positive and negative deviations of energy deposition, spatially correlated with undulations of the ED boundary in both the $x$ and $y$ directions, are also evident inside the blood vessel.

However, the predicted energy uptake by the vessel, ${\langle a\rangle }_{\text{L}}=4.15\mathrm{J}/\mathrm{mm}$ (from model G), is practically identical to the value obtained with the flat ED boundary (Table 2). The deviation in model D is also the same as seen in the former case at similar spatial discretization steps, ${\langle \delta \rangle }_{\text{L}}=1\%$. Similarly, the assessed values for total reflectance of the modeled geometry amount to 24.5% and 24.6% from models D and G, respectively, and do not differ significantly from the values obtained when considering the flat ED boundary.