1.

Introduction

Knowing the light fluence distribution or the amount of light absorbed inside biological tissue is useful in many applications, such as photoacoustic (PA) imaging, diffused optical tomography, Raman scattering, and fluorescence imaging. Monte Carlo (MC) simulation for light propagation in biological medium solves the radiative transfer equation numerically. Therefore, it is considered to be the gold standard for predicting the fluence distribution in biological tissue for many optical imaging modalities. MC for light propagation in biological medium was introduced by Wilson.¹ It was modified by many others for usability.^2–4 MC modeling of light transport in multilayered tissues (MCML) coded in standard C, brought to public domain by Wang et al., has been extensively used for various studies.⁵ Error percentages are higher in MATLAB® ported MCML, but the graphical outputs are easier to obtain.⁶ Initially, the drawback of MCML was the long time taken for each simulation. As mentioned by Zhu and Liu, with high-end computation facilities available these days, computation time has reduced from many hours to a few minutes.⁷ Multicanonical MC⁸ has been introduced, which is a speeded-up form of classical MCML. Moreover, variations of MCML simulations for fluorescence propagation and Raman generation were also reported, which are computationally time demanding.^9,10

Another drawback of MCML was its inability to model geometries other than layers. Tissue models with embedded objects have been published for various applications.⁷ Currently, hybrid models are used for simple geometries, such as skin tumor, where the tumor is considered to be cuboidal.¹¹ Another hybrid, MC and diffusion theory, was developed by Golshan et al. for the study of light propagation in skin.¹² MC simulations for photodynamic therapy for tumor assume the embedded object to be of spherical shape.¹³ For modeling of illumination configuration for light focusing in tissue with blood vessels and capillaries, the embedded object is modeled as cylindrical in shape.^14,15 In the above cases, the refractive index of the embedded object is matched with that of the surrounding tissue. In our previous work, the refractive index mismatch was taken into consideration for an embedded sphere for simulation of light delivery configuration of PA imaging of sentinel lymph nodes.¹⁶ Modeling of light propagation in a cylinder with embedded cylinder depicting blood vessel and sphere within a sphere depicting head has also been reported, which are specific to applications.^17–20 MCML was also modified to specify optical properties for each grid instead of an entire layer to model skin.²¹

Mesh-based MC (MMC) and GPU-based MC by Fang et al. has been used to study light propagation in a rat model and adult/neonatal brain, where the inhomogeneity is of irregular geometry.^22,23 Input to MMC is meshes generated either from images (grayscale or binary) or by specifications of geometric dimensions and location. Mesh generation demands expertise in MATLAB® (or Octave) programming. Programs have to be written to decode the output of MMC, which is in the form of vertices, nodes, and faces of meshes. Due to the complexity involved in using MMC in comparison with MCML, tumor is modeled as layer.²⁴ Since the MCML code is readily available, customization, such as Raman scattering, Tyndall effect, and fluorescence imaging, needs lesser efforts and can be done faster. For modeling photon reflectors one has to stick to MCML, as reflecting the photons that are transmitted or reflected at the boundaries is easier in MCML than MMC. To model photon reflector in MMC, the tracked photons that reach the reflection surface have to be remitted once the run is over.^9,16 CONV MCML is also readily available if one wants to model broad laser beam rather than pencil beam, whereas MMC is specific to one pencil beam (point source).²⁵ MMC can model regular and irregular geometries, but nonconventional photon tracing, such as Raman fluorescence, is not possible. So there is a need for MCML with the capability of simulating multilayered tissues with embedded objects. The hybrid models are easier to model than MMC, but they reduce accuracy.

In this work, MCML with object of regular geometry (sphere, ellipsoid, cylinder, or cuboid) was embedded in the layered surrounding tissue; we refer to this as MMCL with embedded objects (MCML-EO). The refractive index mismatch between the embedded object and the surrounding tissue is also handled in this simulation. Tumors and lymph nodes can be modeled close to sphere or ellipsoid. Red blood cells are ellipsoidal in shape. Simulation of blood vessels and bone need cylindrical geometry embedded in tissue. For the completeness, modeling of cuboidal geometry is also discussed. The MMC is also run with embedded objects with similar optical properties and dimensions for comparisons.

2.

Simulation Setup

In MCML, photons are traced in the medium with known optical properties, such as absorption coefficient (${\mu }_a$), scattering coefficient (${\mu }_s$), scattering anisotropy ($g$), and refractive index ($n$) of the medium. Figure 1 depicts the flow chart of MCML modified for embedded object. Launched photon in MCML takes a random step-size and checks if the layer boundary is hit with this step-size based on the $z$ direction cosine ($u_z$) of the photon. If $u_z$ is negative, then the boundary check is for the layer above the photon’s current layer. If $u_z$ is positive, then the boundary check is for the layer below the photon’s current layer. If there is no boundary hit, then the photon hops to the next scattering site with the step-size and drops weight based on absorption coefficient of the medium. When the photon hits the boundary, the step-size is recomputed to find the distance between photon’s current location and the boundary. Photon after taking the recomputed step-size moves to the boundary and checks if it has to cross the boundary (transmit) or not (reflection). The decision is governed by Snell’s law based on the refractive indices of the layers across the boundary. The direction cosines of the photon either on reflection or on transmission are updated based on Fresnel formula. Photon’s direction cosines remain unchanged under matched boundary conditions. The photon moves the remaining part of the step-size either in the same layer or in the transmitted layer based on whether the photon is reflected or transmitted. The photons are tracked till they die. A photon’s tracking is terminated when its weight is $<{10}^{-6}$ (variable) using Russian roulette or when it escapes from the simulation geometry into the outer medium (launch surface or transmission surface, typically air or water).

Fig. 1

Modified Monte Carlo modeling of light transport in multilayered tissue (MCML) flow chart for the embedded object. $s$ is the dimensionless step-size, $\xi$ is a random number, $d_b$ is the distance between current location of the photon and the layer, $d_{ob}$ is the distance between current location of the photon and the object, and ${\mu }_t$ is the total interaction coefficient.

For embedded objects there are two challenges: one is computing the distance of the photon from the scattering site to the object boundary and the other is determining the direction cosine of the photon after it is reflected/transmitted at the object boundary. Both the problems arise due to the curved nature of the sphere, cylinder, and ellipsoid object.²⁶ In case of cuboid, the calculation of the distance to the boundary is much simpler, and it is explained in the original MCML for layers.⁵ We used the same method here as well for the embedded cuboid. The rest of the section explains how to find out the distance between the photon’s current position and the boundary of the object (sphere, ellipsoid, and cylinder), and whether the photon is going to hit the boundary with the step-size it needs to travel.

Equation of an origin-shifted sphere centered at $C(C_x,C_y,C_z)$ with radius $r$ is

Equation of an origin-shifted cylinder aligned along the axis vector $(1,0,0)$, centered at $C(C_x,C_y,C_z)$ with radius $r$ is

Equation of an origin-shifted ellipsoid centered at $C(C_x,C_y,C_z)$ with radius along $x$, $y$, and $z$ axis as $(r_x,r_y,r_z)$ is

A ray of origin $(O_x,O_y,O_z)$ with direction cosines $d(d_x,d_y,d_z)$ parameterized by distance $t$ is represented by the equation $(O+t\circ d)$, where “∘” implies multiplication of scalar $t$ with every element in vector $d$. Any point on the ray at distance $t$ from origin of ray is given by $[(O_x+t*d_x),(O_y+t*d_y),(O_z+t*d_z)]$. The intersection point between the ray and the curve object has to satisfy the equation of ray and the equation of curved object. Conversely, the ray would intersect with a curved object, if the solution to the quadratic equation in $t$ [Eq. (4)] is real. The quadratic equation is obtained by substituting the ray equation $(O+t\circ d)$ in the equation of sphere, cylinder, and ellipsoid, which are quadratic equations themselves [as seen in Eqs. (1) to (3)]. In Eq. (4), $t$ is the distance the ray has to travel to hit the curved object. Numerical values $I$, $J$, and $K$ are the coefficients of the equation, which vary for the three geometries. Equation (4) is given as

In case of a sphere

In case of a cylinder aligned along axis $v(v_x,v_y,v_z)=(1,0,0)$ (in our case)

In case of an ellipsoid

In all the equations ($*$) denotes the multiplication of two numbers. If $(J^2-4*I*K)>0$, then the ray does not intersect the curved object. In that case, there is no hit. Since there are two solutions to the quadratic equation, there are two points of intersection. The distance between the points of intersection $O+t\circ d$ and origin of ray $O$ is computed. If the distance is greater than step-size, there is a boundary hit.

Once the photon reaches the boundary, it will undergo either reflection or transmission. For layer interface it is easier to find the normal to the tangential plane. In case of layer, the normal to tangent plane coincides with the global $z$ axis. However, in case of curved geometry, the normal to tangent plane (normal to the incident plane in case of cuboid) does not coincide with the global coordinate $z$ axis; therefore, it needs to be transformed into a local coordinate system whose $z$ axis matches with the normal to tangent. Hence, reflection/transmission is done in local coordinate system and then converted back to the global coordinate just as it is done for spinning of photon in MCML. For a photon of polar and azimuthal angles $({\theta }_0,{\varphi }_0)$, the steps involved in conversion of global coordinate system $(u_x,u_y,u_z)$ to local system $(u_x^{\prime },u_y^{\prime },u_z^{\prime })$ are as follows:

The formulation for the above concept is

If $u_z\to 1$, then $u_x^{\prime }=\sin {\theta }_0*\cos {\varphi }_0$, $u_y^{\prime }=\sin {\theta }_0*\sin {\varphi }_0$, and $u_z^{\prime }=\mathrm{sgn}(u_z)*\cos {\theta }_0$, where

The normal to the tangential plane is needed for Fresnel computations. Angle of incidence $({\propto }_i)$ is ${\cos }^{-1}|u_z|$. When there is a refractive index match, angle of transmittance $({\propto }_t)$ is equal to $({\propto }_i)$. In case of refractive index mismatch, Snell’s law [Eq. (9)] is used to compute $({\propto }_t)$.

During transmission, the direction cosines of the photon are updated to $[u_x*(n_i/n_t),u_y*(n_i/n_t),\mathrm{sgn}(u_z)*\cos {\propto }_t]$. For reflection, only $u_z$ of the photon is negated.

Figures 2(a) and 2(b) give a few examples of ray interaction at layer boundaries; line MN is the normal to the plane of incidence PQ. MN is parallel to $z$ axis since plane of incidence is perpendicular to it.

Fig. 2

Examples of ray boundary interaction under matched and mismatched boundary conditions. (a) Reflection and transmission under mismatched boundary condition for rays AB, DE, and HI, which are incident at 30, 30, and 60 deg. (b) Matched boundary condition for ray AB. MN is the normal and PQ is the plane of incidence.

This section will explain how the normal to the tangential plane is computed in case of curved objects (sphere, ellipsoid, and cylinder) as this is the reference to compute critical angle and transmittance angle. In case of sphere, the line joining the center of the sphere to the point of tangency $P(P_x,P_y,P_z)$ on the surface $[(P_x-C_x),(P_y-C_y),(P_z-C_z)]$ is the line normal to the tangential plane. This is pictorially represented in Fig. 3(a). A local Cartesian coordinate system is created with the normal line (line perpendicular to the tangential plane) as $z$ axis. The direction cosines are calculated in the rotated coordinate system; angle of reflectance or transmittance is then computed and again converted back to global coordinate system. For an $x$-axis aligned cylinder, the normal line to a tangent at a point $P(P_x,P_y,P_z)$ on surface is $(0,P_y-C_y,P_z-C_z)$ [Fig. 3(b)]. For an origin-shifted ellipsoid [Fig. 3(c)], the normal to the tangent at a point $P$ is $\{[(P_x-C_x)/r_x^2],[(P_y-C_y)/r_y^2],[(P_z-C_z)/r_z^2])$. Equation of tangential plane for ellipsoid was derived from Eq. (3) to find its normal.

Fig. 3

Pictorial representation of incidence plane and normal plane considered in executing Fresnel formula for MCML with embedded objects (MCML-EO). Point $C$ is the center of the object and $P$ is the point on the surface of object where the photon reaches when it hits the object boundary. (a) Sphere with tangent (tangential plane in case of three dimension) and normal $N$, (b) cylinder with tangent and normal, (c) ellipsoid, and (d) cuboid depicting two cases—positive $x$ axis ($P_1$) and positive $z$ axis ($P_2$) being the planes on which photon has reached.

Now looking into cuboid, which is a collection of six planes, based on the direction cosines of the photon, the boundary check is done for the planes that lie in these directions. Cuboid center $(C_x,C_y,C_z)$ and dimensions [length (lt) along $x$ axis, breath (bt) along $y$ axis, and height (ht) along $z$ axis] are read from the input file. The six planes are defined as follows:

When $u_z$ is positive, the plane below the center of cuboid is checked. When $u_z$ is negative, distance computation is with respect to the plane above the center. Sign of $u_x$ and $u_y$ decide if the distance check should be for the plane on positive axis or negative axis. Reflectance angle and transmittance angle are calculated with respect to the $x$, $y$, and $z$ axis based on the plane that the ray intersects. This is shown in Fig. 3(d).

MCML-EO for sphere, ellipsoid, and cuboid was filled with methylene blue for demonstration, as sphere/ellipsoid filled with methylene blue can be used to model the sentinel lymph nodes in PA imaging. Embedded cylinder of radius 0.25 cm was assumed to have properties of blood, which can be used as the model for blood vessels. The dimensions and the optical properties for the embedded objects and the surrounding tissue are given in Tables 1 and 2, respectively. For all the layers and embedded objects, $g=0.9$. The outer medium (launch surface and transmit surface) was taken as air with refractive index 1.0. The depth of the surrounding tissue layer was taken as 6 cm. Simulation geometry for the four embedded objects are given in Fig. 4. A pencil beam light was launched at the origin $(0,0,0)$ along the $z$ axis with direction cosines $(0,0,1)$. Weight dropped was accumulated in $x$, $y$, and $z$ grids. For fluence map, weight dropped was divided by absorption coefficient. Number of grids in each axis is 301 and the size of the grid was 0.02 cm. The volume spanned is $-3$ to $+3\mathrm{cm}$ in $x$ and $y$ axis. Grid covers 0 to 6 cm in $z$ axis. For display (either absorbance map or fluence map) $y=0$ plane is presented. Contour plot of the fluence (averaged over five grids) is for 5-dB spacing. The geometries with same optical properties and dimensions (approximately) are modeled using MMC for comparison.

Table 1

Dimensions, location, and optical properties of the embedded objects for Monte Carlo modeling of light transport in multilayered tissue with embedded objects (MCML-EO) and mesh-based Monte Carlo (MMC).

Table 2

Optical properties of various layers used in the MCML-EO and MMC simulation model at 664-nm wavelength.

Fig. 4

Schematic diagram of the simulation geometry. Skin of semi-infinite depth with $n=1.4$, ${\mu }_a=0.2525{\mathrm{cm}}^{-1}$, and ${\mu }_s=254{\mathrm{cm}}^{-1}$. (a) 1-cm diameter spherical object placed 0.7 cm below the skin filled with methylene blue of $n=1.3$, ${\mu }_a=1.7049{\mathrm{cm}}^{-1}$, and ${\mu }_s=180{\mathrm{cm}}^{-1}$. (b) 0.5-cm diameter cylinder placed 0.95 cm below the skin filled with blood of $n=1.35$, ${\mu }_a=2.10{\mathrm{cm}}^{-1}$, and ${\mu }_s=773{\mathrm{cm}}^{-1}$. (c) 1-cm diameter along $z$ axis ellipsoidal object placed 0.7 cm below the skin filled with methylene blue. (d) 1-cm (height) cuboid placed 0.7 cm below the skin filled with methylene blue. $g=0.9$ for all the layers and objects. Inputs are in ${\mathrm{cm}}^{-1}$ for MCML-EO and in ${\mathrm{mm}}^{-1}$ for mesh-based Monte Carlo (MMC).

MMC mex files along with preprocessing and postprocessing MATLAB® functions were downloaded and compiled from website.²⁷ Ray tracing using Plücker co-ordinate system was implemented to trace photons in nonhomogeneous medium. In MCML-EO the surrounding medium was infinite along $x$ and $y$ axis. Simulating infinite medium in MMC is not possible; thus, in MMC, the objects were embedded into a relatively large box of dimension $20\times 20\times 6\mathrm{cm}$ (which can be assumed approximately infinite compared to the inclusion dimension). All inputs are in millimeters in MMC. The maximum volume for a single mesh was 0.020 cm³. The incident light beam is a point source (pencil beam) illuminating at $(10.1,10.1,0)$ with the incident direction cosine $(0,0,1)$. Optical properties of the surrounding tissue (box) is ${\mu }_a=0.2525{\mathrm{cm}}^{-1}$, ${\mu }_s=254{\mathrm{cm}}^{-1}$, and $n=1.3$ at 664-nm wavelength.²⁸ The optical properties of methylene blue (of concentration 10 μM) are ${\mu }_a=1.7049{\mathrm{cm}}^{-1}$ and ${\mu }_s=180{\mathrm{cm}}^{-1}$.^29,30 Optical properties of whole blood are ${\mu }_a=2.10{\mathrm{cm}}^{-1}$, ${\mu }_s=773{\mathrm{cm}}^{-1}$, and $n=1.4$.^31,32 All simulations were run for ${10}^6$ photons on a desktop with Intel i7 64-bit processor and 8 Gb RAM. Plane $y=100$ is imaged and contoured with 10-dB spacing.

3.

Results and Discussion

The major drawback of MCML of not being able to handle boundaries other than layers is addressed to increase the flexibility of the simulation tool. The contribution of the work is consolidation of distance calculation to the curved boundaries of sphere, cylinder, and ellipsoid and computation of tangential plane and its normal for execution of Fresnel formula. MCML modified for embedded objects is compared with MMC in terms of fluence map. First, MCML-EO simulations were run for all four types of embedded objects, considering both boundary refractive-index matched and mismatched conditions between the embedded object and the surrounding tissue. Table 3 shows the absorption inside the embedded object for four cases. As observed, the difference in absorbance (error) within the object ranges from 4 to 6% between refractive-index matched and mismatched cases. Absorbance within the embedded object ranges from 0.03 to 0.007. Note that the absorption is more when there is refractive-index matched boundary (which is not a practical scenario) compared with mismatched boundary condition. Absorption of light within embedded object is of interest in many applications; for example, in case of sentinel lymph node detection, one needs to optimize the total absorbance inside the node to achieve high signal-to-noise ratio.¹⁶

Table 3

Absorbance within the embedded object under matched and mismatched boundary conditions.

Figure 5 shows the absorbance map (in log scale) of the four MCML-EO simulations. The dotted black lines represent the boundary of the embedded object. Figures 5(a) to 5(d) show absorbance along $y=0$ plane in case of sphere, cylinder, ellipsoid, and cuboid, respectively. Number of pixels (grids) along $x$ and $z$ axis are 300 each. Absorption map is important to study effects of illumination during phototherapy.¹⁴ Once again, one can play with the light delivery configuration, optical properties, and various object sizes to see what kind of absorption map one wants to achieve. This makes MCML-EO a very flexible tool for various clinical uses.

Fig. 5

Absorption map from Monte Carlo for turbid medium with embedded object along $y=0$ plane. (a) to (d) Sphere, cylinder, ellipsoid, and cuboid. Black dotted lines are boundaries of the embedded object.

Figure 6 shows the log of fluence map from MCML-EO and MMC. Figures 6(a) to 6(d) show the fluence map from MCML-EO, where each pixel is 0.2 mm and the number of pixels is $300\times 300$. Figures 6(e) to 6(h) are the fluence maps from MMC, which cover the meshes from 7 to 13 cm. Log of fluence distribution is represented in the color bar. The contours of the fluence distribution from MMC and MCML-EO are shown in Figs. 6(i) to 6(l). MCML-EO contours are smoother compared to that of MMC, probably due to the much finer grids in MCML-EO compared to the meshes in MMC. Again, the dotted black lines are the boundaries for the object embedded. We can see agreement between MMC and MCML-EO. The simulation parameters do not exactly match in both MMC and MCML-EO as mentioned earlier. MCML-EO is for semi-infinite medium, whereas MMC is done for bounded medium (large compared to the embedded object size). But for qualitative comparison, it is acceptable.

Fig. 6

Fluence map of MCML-EO ($y=0$ plane) and MMC ($y=100$ plane) along with contours. (a) to (d) Fluence map of sphere, cylinder, ellipsoid, and cuboid from MCML-EO. (e) to (h) Fluence map of sphere, cylinder, ellipsoid, and cuboid from MMC. (i) to (l) Contours of fluence distribution of sphere, cylinder, ellipsoid, and cuboid from MCML-EO (5 dB line spacing) and MMC (10 dB line spacing).

Average runtime for MCML-EO is 4 min per geometry. Average runtime for MMC per geometry is $\sim 15\min$. There are a few seconds (40 s) spent in mesh generation in case of MMC, but MCML-EO needs no preprocessing. In applications where the interest is to know total absorbance in the embedded object (as is the case in illumination configuration for PA imaging of sentinel lymph node¹⁶), there is no postprocessing involved for MCML-EO as the weight dropped in the embedded object is accumulated and printed in the output file. But MMC requires postprocessing to track the weights dropped in the meshes contributing to the embedded object. For the fluence or absorption maps seen in Figs. 5 and 6(a) to 6(d), the postprocessing time goes up to 15 min in MCML-EO, but it takes only few seconds (45 s) in the case of MMC. With the increase in the volume of interest, the mesh generation time also increases. With increased mesh counts, the computation requires higher memory and increased CPU speed. The flexibility of MMC is that the number of meshes can be more for the embedded object and sparse in the tissue.

In case of finer inhomogeneity, such as brain modeling, one has to proceed with MMC, where the simulation geometry is given as volumetric data in the form of images. MCML-EO is $\sim 200\text{times}$ faster than MMC with respect to run time; however, embedding complex geometries will be a challenge in MCML. Optical properties can be assigned on nodal basis, which increases accuracy of modeling in the case of MMC. But for simpler modeling, such as blood vessels (cylinder), lymph nodes (ellipsoid), tumor (sphere), bone, and capillaries, one can use the MCML-EO.

4.

Conclusion

MCML-EO of geometries like sphere, cylinder, ellipsoid, and cuboid increases the flexibility of MCML for simulating more realistic structures of biological systems. Raman modeling, effect of photon reflectors, and distributed illumination sources can be easily studied for geometries other than layers. Total computation times for both MCML-EO and MMC are approximately the same. However, larger geometries need more memory for generating meshes in the case of MMC. From the programming language point of view, MMC requires knowledge of MATLAB®/Octave (for three-dimensional pre- and postprocessing). In MCML-EO outputs are either matrices (diffused reflectance, transmittance, and absorbance) or variable (absorption within embedded object). Standard C version of MCML code is available online, so remodeling it is easier compared to MMC, where preprocessing and postprocessing files are in MATLAB® and only mex files of MC are available. In this work we have considered refractive index mismatch between the embedded object and the surrounding medium. Both MCML-EO and MMC are capable of producing similar fluence map and absorption map. It is also seen that if boundary refractive index mismatch is ignored, 4 to 6% error is recorded in the total absorption within the embedded object. Future work is to combine MCML for all regular geometries, which would make it a more user-friendly tool, such as Online Monte Carlo, GNEAT.^33,34 Having discussed MCML-EO and MMC, the users can choose the tool that best suits their requirement. If one is interested in irregular geometry, then MMC should be the choice. Embedded object of defined geometry within a bounded medium can be simulated using MMC and MCML-EO. If the application requires nonconventional photon tracking, such as reflector design, Raman propagation, and fluorescence, then MCML-EO would be an apt choice.