1.

Introduction

Significant work has been done in characterizing the bulk optical absorption, scattering, and refractive properties of biological tissues separated into reasonable classes, such as skin, brain, bone, blood, etc.¹ This work provides the foundational data for noninvasive optical analysis, which can, in theory, estimate the physiology and chemistry of an illuminated sample to a given degree of precision. That is, with the optical properties of the major constituents understood, it is possible to predict how light will interact with an anatomical part, such as a human finger, and extract from this interaction information about the blood and tissue. Of course, as has been summarized elsewhere,¹ this is a difficult challenge since any optical property will change from person to person, from site to site on one person, and from time to time on any given site. It is, therefore, essential that detailed time- and position-sensitive data be modeled to understand what information can and cannot be extracted if noninvasive techniques are to be extended to measuring biological chemistry and/or physiology at increasingly finer detail.

One approach to understand how light propagates through tissues is to solve the diffusion equation as Steinke and Shepherd² did with whole blood. However, the complex geometries that exist in tissue make this approach impractical without simplifying approximations. It has, therefore, been a common practice to use less constrained Monte Carlo simulations to model photon interactions in specific applications. We have adopted this approach, as have others,^3–12 because of its relative simplicity and accuracy.

Numerous studies involving Monte Carlo models have addressed the various geometries that are present within the human body, particularly with blood and tissue. Wang et al. created a multilayer grid model to explore diffuse reflectance and transmittance initiated by an infinitely narrow beam incident on the upper layer.⁸ Others have created more detailed models with more layers to explore skin reflectance^12–14 and blood flow with layers.^10,11 Some have explored the effects of complex boundaries.¹³ Still, others have modeled blood vessels within tissue,¹⁵ including using cylindrical models¹⁶ and by grouping the blood cells into discrete regions.¹⁷ This is a rich field, and we touch upon only a few of the many studies done within it.

As an exploration of the fundamental aspects of biological modeling itself, we have developed a Monte Carlo model to examine a homogeneous mix of basic biological material to understand better the complexities that may arise as structures are modeled in increasingly finer detail. We chose the finger for our own research interests, but the concepts and results presented in this article are general. We are interested in how near-infrared light transmits, particularly attenuates, through a semi-infinite slab of such a material. Others, such as Schmitt and Kumar,¹⁸ have used mathematical relationships and other physical principles applied to a mixed model to study how the scattering and back-scattering coefficients as well as the anisotropic factor vary based on the scattering contributions of the individual scatterers within a tissue as a whole. In the study reported here, we use Monte Carlo simulations to simulate light interaction with a mixture of two whole materials. Our aim is to find what effect mixing materials has on the attenuation of light within such a model as opposed to finding the spectral dependence of the scattering coefficients.

2.

Monte Carlo Model

Our model begins with launching individual photon packets normally down into a specified material, where they interact in accordance with the material’s scattering and absorbing coefficients, ${\mu }_s$ and ${\mu }_a$, respectively. The propagation distance, $t$, between interactions is given as follows in terms of the interaction coefficient ${\mu }_t$:^6,7

The scattering of the photon packet is determined using directional cosines that include an anisotropic factor, and intensity is removed from the photon packet according to the absorption and scattering coefficients, as was done by Prahl et al.⁶

We set up our model in a Cartesian coordinate system, where the source is at the origin, and photon packets are launched in a collimated beam in the $-z$ direction, normal to the interface. A circular detector 1 mm in diameter is placed at a location $(0,0,-d)$, where $d>0$. We measure the intensity (termed “weight”) of the photon packets that reach the detector region. As we increase the thickness, $d$, of the material, the transmitted intensity should follow a modified form of Beer’s Law for highly scattering materials given as

3.

Methods and Procedures

There are cases where layered models are preferred over mixed models and vice versa. Blood and plasma are best represented in the finger by a mixed model since they travel together through blood vessels and veins. For epidermis and blood, a layered model is more useful in many regions, whereas a mixed model may be more accurate in other places. Although epidermis is a bloodless tissue,²⁰ some areas of the epidermis layer are penetrated by capillaries, deviating from a layered model.²¹ A mixed model, although not completely accurate, approximates this penetration and is one step closer toward adding greater complexity to our models.

We chose to examine mixtures of epidermis, blood, and plasma because of their prevalence in the situations of interest to us^21,22 and for simplicity in understanding how light attenuates differently for layered and mixed models. Since water makes up $\sim 92\%$ of plasma,²³ we approximate plasma with water.

We ran semi-infinite models (finite in $z$, infinite in $x$ and $y$) of pure material, layered material, and homogeneous mixed material using Mathematica. Illumination was from a collimated beam of monochromatic light at $\lambda =810\mathrm{nm}$. We ran each simulation starting at a thickness $d$ of 0.010 mm and repeated the simulation at new thicknesses increased in increments of 0.100 mm up to a maximum thickness of 5.00 mm. We then plotted the transmitted intensity versus thickness $d$ and fit it with an exponential curve. The derived attenuation constant was then compared to $\alpha$ given in Eq. (3) for the respective tissues.

3.1.

Semi-Infinite Pure Tissue Model

To verify our approach, we first modeled semi-infinite pure tissue structures of epidermis, water, and whole blood, separately. These simulations successfully modeled the change in detected intensity versus the thickness of the sample for a modified form of Beer’s law in materials that exhibit high scattering properties [Eq. (3) and Fig. 1]. The derived $\alpha$ values agreed with published values from diffusion models after they are normalized to our approach by dividing out the $\sqrt{3}$ factor in Eq. (4) (see Table 1).

Fig. 1

Results of the Monte Carlo simulation for a semi-infinite slab of whole blood, represented by solid blue points. Error bars from Poisson statistics are shown. The modeled attenuation coefficient ${\alpha }_{\text{derived}}$ is derived from the solid black curve. Substituting the theoretical normalized value, ${\alpha }_{\text{normalized}}$, into modified Beer’s law [Eq. (3)] gives the dashed red line. One can see that the derived exponential decay curve of the data (solid black) is in good agreement with the modified form of Beer’s law and Eq. (4) (dashed red).

Table 1

Absorption (μa), scattering (μs), transmission (μt), and bulk attenuation (α) coefficients and anisotropy factor (g) of epidermis (Ref. 19), water (Ref. 24), and whole blood [Hct=0.45−0.46, oxygenation >98% (Ref. 19)] at a wavelength of 810 nm. The values for αnormalized are the values calculated using Eq. (4), normalized by dividing out the 3. The values for αderived are the attenuation coefficients we obtained from fitting an exponential decay to the results from the simulations (see Fig. 1). Comparing αnormalized to αderived, we see that the two values agree quite well.

Material	${\mu }_a$ (${\mathrm{mm}}^{-1}$)	${\mu }_s$ (${\mathrm{mm}}^{-1}$)	${\mu }_t$ (${\mathrm{mm}}^{-1}$)	$g$	${\alpha }_{\text{normalized}}$ (${\mathrm{mm}}^{-1}$)	${\alpha }_{\text{derived}}$ (${\mathrm{mm}}^{-1}$)
Epidermis	4.0	42.0	46.0	0.85	6.419	7.815
Water	$2.00\times {10}^{-3}$	$3.00\times {10}^{-7}$	$2.0003\times {10}^{-3}$	0.9	$2.00\times {10}^{-3}$	$1.14\times {10}^{-3}$
Whole blood	0.65	69.0	69.65	0.989	0.957	0.931

3.2.

Semi-Infinite Layered Tissue Model

We next modeled layers of two different types of materials to compare against previously published results. This was accomplished by dividing up a semi-infinite volume into two separate layers, each composed of a different material. The interface between the two layers was given a specific coordinate, say $d_1$, in the $xy$-plane (Fig. 2).

Fig. 2

Semi-infinite layered model for two tissue layers.

We modeled water and whole blood with water being the top layer and blood being the lower layer. We varied the percentage of each component from $10\%/90\%$ to $90\%/10\%$ and generated the same intensity versus thickness curves as earlier. We derived a value for $\alpha$ as earlier and compared them with theory (Fig. 3).

Fig. 3

Bulk attenuation coefficient $\alpha$ as a function of volume fraction of water in a layered model of water, as the layer near the source and blood as the layer near the detector. Statistical error bars are shown for each point but are too small to see for most points. The coefficient of determination, $R^2$, is the statistical value that represents how close the data points follow the linear relationship between $\alpha$ and the volume fraction of water (dashed line). The results for this model yield $R^2=1.0\pm {10}^{-4}$. A value of $R^2=1$ represents a perfect correlation.

One would expect there to be a linear relationship between the total attenuation and the volume fraction of each layer since the individual $\alpha$ values mathematically add together. That is, if $d_1$ is the thickness of the layer of water, $d_2$ is the thickness of the layer of blood, ${\alpha }_{\text{water}}$ is the attenuation coefficient of water, ${\alpha }_{\text{blood}}$ is the attenuation coefficient of blood, $x$ is the fraction of the total volume that is water, and $y$ is the fraction of the total volume that is blood, then multiplication of the two curves yields

Eq. (6)

$$\begin{gathered} I(d_1,d_2)=I_0{\mathrm{e}}^{-{\alpha }_{\text{water}}{d}_1}{\mathrm{e}}^{-{\alpha }_{\text{blood}}{d}_2} \\ I(d)=I_0{\mathrm{e}}^{-{\alpha }_{\text{total}}d}, \end{gathered}$$

where ${\alpha }_{\text{total}}={\alpha }_{\text{water}}x+{\alpha }_{\text{blood}}y={\alpha }_{\text{water}}x+{\alpha }_{\text{blood}}(1-x)$ and $x+y=1$. Figure 3 shows that the expected linear relationship between ${\alpha }_{\text{total}}$ and the volume fraction of water is in excellent agreement with theory. We simulated a layered model of blood and epidermis as well and produced the expected linear relationship between volume fraction of either substance and the total attenuation coefficient in this case, too.

4.

Results and Discussion

After making the adjustments above, we ran the simulation for a mixture of water and blood and blood and epidermis. The results are in Figs. 4 and 5 and are plotted as ${\alpha }_{\text{total}}$ versus volume fraction. These figures show that attenuation in a mixed medium does not have the linear relation with material fraction that is found in layered models (see Table 2).

Fig. 4

Bulk attenuation coefficient $\alpha$ as a function of volume fraction of water in a mixed model of water and blood. Error bars from Poisson statistics are smaller than most data symbols. The coefficient of determination, $R^2$, is 0.889 and represents the value of a linear relationship between $\alpha$ and the volume fraction of water (red, dashed line) if one existed.

Fig. 5

Bulk attenuation coefficient $\alpha$ as a function of volume fraction of epidermis in a mixed model of epidermis and blood. Error bars from Poisson statistics are smaller than most data symbols. The coefficient of determination, $R^2$, is 0.997 and represents the value of a linear relationship between $\alpha$ and the volume fraction of blood (red, dashed line) if one existed.

The linear relation for ${\alpha }_{\text{total}}$ in Beer’s law [Eq. (6)] is a consequence of the materials being physically and statistically independent. In the layered model, the materials are separated into distinct sections, where the propagation distance between interactions depends solely on the material the photon packet is in at any given moment [Eq. (1)]. Statistically, we can say that the probabilities of interaction of the different materials are independent of one another, as seen in the modified form of Beer’s law [Eq. (6)]. However, as soon as the materials are mixed together, the propagation distance between interactions is now dependent on all of the materials within the total interaction volume [Eq. (9)]. Therefore, the mixed model no longer yields statistical independence between the differing materials, and the linear relationship in ${\alpha }_{\text{total}}$ no longer holds as it did in the layered model.

Table 2

Comparison of total attenuation coefficient of a semi-infinite slab of 50% blood and 50% epidermis by volume derived from three different models: Monte Carlo multilayered (MCML), diffusion theory homogeneous mixed, and Monte Carlo homogeneous mixed.

In terms of the bowling-ball–ping-pong-ball analogy discussed in Sec. 3.3, a material with the larger ${\mu }_t$, or larger interaction cross-section, would be represented by the bowling ball, whereas the one with the lower ${\mu }_t$, or smaller interaction cross-section, would be represented by the ping-pong ball. If we took a mix of balls that was 50% bowling balls and 50% ping-pong balls by number and then separated them into two volumes, such as in the layered model, we would find that the bowling balls would occupy more volume than the ping-pong balls. Thus, ${\alpha }_{\text{total}}$ would not prorate as 50% ${\alpha }_{\text{bowling}}$ and 50% ${\alpha }_{\text{ping-pong}}$ but rather have a higher weighting of ${\alpha }_{\text{bowling}}$ than ${\alpha }_{\text{ping-pong}}$. As a result, the $\alpha$ of the material with the higher ${\mu }_t$, in this case, the bowling ball, would dominate the curve for ${\alpha }_{\text{total}}$, an argument that is also mentioned, but not fully characterized, by Schmitt and Kumar.¹⁸ The model with a mixture of water and blood is comparable to this analogy (Fig. 4). In this case, if blood is subscript “1” and water is subscript “2,” ${\mu }_{t,1}>{\mu }_{t,2}$ and ${\alpha }_1>{\alpha }_2$.

In another situation, we could have ${\mu }_{t,1}>{\mu }_{t,2}$ while ${\alpha }_1<{\alpha }_2$. Since $\alpha$ of a given material is dependent on the anisotropic factor, $g$, and ${\mu }_t$ is independent of $g$, a situtation like this could exist if the anisotropic factor of material 2 is lower than that of material 1 [Eqs. (4) and (5)]. This scenario can be represented by large, light-weight plastic balls (material 1) and heavy shot-put balls (material 2). Since the plastic balls are larger, they have a higher interaction cross-section (${\mu }_{t,1}>{\mu }_{t,2}$), meaning it is much more likely for a photon packet to interact with them. On the other hand, since they are lightweight, they do not attenuate the intensity of the photon packet as much as the heavier shot-put ball would, which corresponds to a lower attenuation coefficient (${\alpha }_1<{\alpha }_2$). Since ${\mu }_{t,1}>{\mu }_{t,2}$, ${\alpha }_1$ will dominate the curve for ${\alpha }_{\text{total}}$. The model with a mixture of blood and epidermis is comparable to this analogy (Fig. 5).

The direction in which the curve for ${\alpha }_{\text{total}}$ bends is determined by comparing ${\mu }_t$ of the two materials. The material with the higher ${\mu }_t$ will have a greater weighting in the total attenuation coefficient, ${\alpha }_{\text{total}}$. Thus, if the material with the greater ${\mu }_t$ also has a greater $\alpha$, then the curve for ${\alpha }_{\text{total}}$ will bend above the linear relationship (Fig. 4). If the material with the greater ${\mu }_t$ has a smaller $\alpha$, then the curve will bend below the linear relationship (Fig. 5). How much the curve for ${\alpha }_{\text{total}}$ deviates from the linear relationship depends on how much the values of ${\mu }_t$ for the two materials differ. In the case of the water–blood mix, ${\mu }_{t,\text{blood}}\gg {\mu }_{t,\text{water}}$ (see Table 1), thus, the curve for ${\alpha }_{\text{total}}$ bends far from the linear relationship (Fig. 4). On the other hand, in the blood–epidermis mixed model, ${\mu }_{t,\text{blood}}>{\mu }_{t,\text{epidermis}}$, but they are of the same order of magnitude (Table 1). As a result, although the curve is more weighted to ${\alpha }_{\text{blood}}$, the deviation from the linear relationship is small.

The approach of Verkruysse et al. for their analysis of the attenuation by blood in a port wine stain is similar to ours although there are some important differences to note.¹⁷ They use absorption and scattering coefficients and anisotropic factors that are linearly scaled according to the volume fraction of each material in the total mixture. The nonlinear aspects of increasing attenuation from greater vesicular blood volume is reflected in their correction factor $C$. They then use these new mix coefficients to directly calculate the attenuation coefficient of a mixed model according to Eq. (4) (see Table 2). We linearly scaled the scattering and absorption coefficients according to the volume fraction of each material in the total mixture, which results in the transmission coefficient being linearly scaled according to volume fraction. However, we only use the rescaled transmission coefficient for the propagation distance and the selection window for which material to interact with since both of these parameters are dependent on all materials within the entire model. We did not rescale either ${\mu }_a$ or ${\mu }_t$ in the absorption factor (${\mu }_{\mathrm{a}}/{\mu }_t$), and we do not rescale the anisotropic factor, $g$. Our model is based on discrete events where a photon packet interacts with one material at each event. Both the absorption and anisotropic factors depend solely on the given event and, therefore, can only be based on the coefficients and factors of that one material.

Returning to the issue of the packing factor and its effects on ${\mu }_s$ introduced in Sec. 3.3, depending on whether a pure material is tightly packed or not, the nonlinear results could be mitigated or intensified. The results shown in Fig. 4 from the blood and water mixed model should not be affected as neither of the materials are tightly packed. However, in the case of the epidermis and blood mixed model, there should be some variations in the results once the packing factor is taken into account because tissue, such as epidermis, is considered to be tightly packed. By adding blood to a model of pure epidermis, we are essentially making the epidermis cells less tightly packed until we reach the model of pure blood. Reducing the packing factor, $W_p$, of epidermis lowers the ${\mu }_s$ of epidermis, and in turn, lowers the ${\mu }_t$ of epidermis as well [Eqs. (2) and (7)]. Since ${\mu }_t$ of epidermis is already lower than that of blood (Table 1), the difference between the ${\mu }_{t}s$ of the two materials becomes even greater, and thus, we should see the nonlinearities become more pronounced based on previous arguments.

5.

Conclusions

We have developed a Monte Carlo homogeneous mixed model and have studied how light should attenuate in such a material. We have found that the linear scaling of attenuation coefficients found in layered models does not hold when materials are homogeneously mixed together. In cases, such as those involving blood or plasma, where there is no clear boundary between constituents or the constituents are mixed together in a fluid, it may be necessary to use a mixed model to derive the total attenuation coefficient with the required accuracy. This is especially true if the ${\mu }_t$ values of the materials differ by an order of magnitude or more; otherwise, a layered model may be sufficient.