Monte Carlo simulation of near-infrared light propagation through homogeneous mixed media

Abstract. Noninvasive blood analysis devices that can measure levels of small constituents of blood are of interest in the medical community. An important step in creating these devices is to understand the interaction of photons with human tissue in increasingly greater physiological detail. Models based on layered biological materials give excellent results for many applications but may not be as accurate as needed when those materials are finely intertwined to the point of resembling a homogeneous mixture. To explore the ramifications of treating materials as layers versus a mixture, we have modeled, using a Monte Carlo technique, the interaction of photons through epidermis, blood, and water arranged both in layers and in a homogeneous blend. We confirm the expected linear relation between photon attenuation and material volumetric percentage in two-layer models. However, when the materials are homogeneously mixed together and volumetric percentage is replaced with interaction volume percentage, this relationship becomes nonlinear. These nonlinearities become significant when the values of the interaction coefficient, μt, differ by an order of magnitude or more.


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. 1 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, 1 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 2 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][4][5][6][7][8][9][10][11][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. 8 Others have created more detailed models with more layers to explore skin reflectance [12][13][14] and blood flow with layers. 10,11 Some have explored the effects of complex boundaries. 13 Still, others have modeled blood vessels within tissue, 15 including using cylindrical models 16 and by grouping the blood cells into discrete regions. 17 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, 18 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.

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, μ s and μ a , respectively. The propagation distance, t, between interactions is given as follows in terms of the interaction coefficient μ t : 6,7 where ξ is random number between 0 and 1 sampled over a uniform distribution and where 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. 6 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 where I 0 is the intensity of light emitted from the source and accounts for the size of the detector, d is the physical separation of the light source and detector, and is the bulk attenuation of the tissue sample. The form of μ 0 s , the reduced scattering coefficient, is where g is the anisotropic factor. 19 We nominally run 50 simulations of varying thickness and from them, we derive the attenuation constant, α. We keep the detector size and emitter intensity the same in all simulations so that the contribution to the scaling factor from the size of the detector is the same for all simulations and thus has no effect on the derived value of α.

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, 20 some areas of the epidermis layer are penetrated by capillaries, deviating from a layered model. 21 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 ∼92% of plasma, 23 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 λ ¼ 810 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 α given in Eq. (3) for the respective tissues.

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 α values agreed with published values from diffusion models after they are normalized to our approach by dividing out the ffiffi ffi 3 p factor in Eq. (4) (see Table 1).

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). 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 α as earlier and compared them with theory (Fig. 3).
One would expect there to be a linear relationship between the total attenuation and the volume fraction of each layer since the individual α 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, α water is the attenuation coefficient of water, α 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 where α total ¼ α water x þ α blood y ¼ α water x þ α blood ð1 − xÞ and x þ y ¼ 1. Figure 3 shows that the expected linear relationship between α 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.

Semi-Infinite Homogeneous Mixed Tissue Model
We next modeled the light propagation through a homogeneous mix of two materials distributed randomly throughout the volume and packed loosely enough so that there was no overlap of interaction cross-sections. For both blood and water, we can accurately make this assumption. For epidermis, a tissue of tightly packed cells, this is not physiologically correct. 20 Even so, we still consider mixed epidermis and blood. Those simulations can be thought of as gradually displacing the tightly packed epidermis cells with an increasing volume fraction of blood, thus reducing the packing factor that Schmitt and Kumar use in defining the μ s of constituents of a particular tissue: 18 where W p is the packing factor, x i is the volume fraction, v i is the volume, and σ i is the cross-section, each of a particular material i. Since the packing factor and scattering coefficient, μ s , are directly related, as the material of interest becomes more tightly packed or less tightly packed, we will see a change in μ s , and thus a change in μ t [Eq. (2)]. For now, we neglect this packing factor in the epidermis and blood mixed model since we are simply interested in the effects on attenuation that a mixed model exhibits.
In deciding which material of a mixed medium was likely to interact with the photon packet, we formed a new μ t;total , which was adjusted to account for the different interaction cross-sections of each type of material. Since each material has different scattering and absorption coefficients that are independent of volume, considering volume fraction alone is inaccurate. To illustrate why, compare the Monte Carlo photons to a bouncy ball moving around in a volume of bowling balls and pingpong balls with enough space between the balls for the bouncy ball to move freely. The bowling balls represent particles of the material with higher scattering and absorption properties. The Table 1 Absorption (μ a ), scattering (μ s ), transmission (μ t ), and bulk attenuation (α) coefficients and anisotropy factor (g) of epidermis (Ref. 19 . 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.  ping-pong balls represent particles of the material with lower scattering and absorption properties. Intuitively, it makes sense that as the bouncy ball moves around the volume, not only are the bowling balls going to absorb more energy than the ping-pong balls but it is also more probable that the bouncy ball will hit a bowling ball rather than a ping-pong ball due to the difference in sizes, or "interaction volumes," of the two balls. (In this analogy, we relate a larger mass directly to a larger attenuation of the moving ball's motion. In reality, a larger mass does not directly relate to attenuating the ball's movement more, but we keep the analogy as an illustrative tool.) Relating the simulation to the bowling-ball-ping-pong-ball analogy, if material 2 has higher scattering and absorption properties (e.g., the bowling balls) than material 1, then the photon packet is more likely to interact with material 2 than material 1. Thus, we rescaled the selection windows of each material to represent the fraction that each material takes of the whole interaction volume instead of the physical volume.
The rescaling goes as follows: let x and y have the same definitions as mentioned above, and let μ t;1 and μ t;2 be the interaction coefficients of materials 1 and 2, respectively. To rescale the volume fraction to the interaction volume of each material in the mixed model, the "cut-off" for the window in selecting between the two materials is cut-off ¼ xμ t;1 xμ t;1 þ yμ t;2 ; (8) so that if the random number is between 0 and cut-off, then the photon packet interacts with material 1, otherwise it interacts with material 2. The propagation distance between interactions, t, given by Eq. (1), is also affected by the difference in interaction crosssections. In Eq. (1), the interaction coefficient, μ t , corresponds to the overall interaction cross-section of a material. A larger μ t corresponds to a larger interaction cross-section, and, in turn, a larger interaction volume. Similar to how Schmitt and Kumar rescaled the total scattering coefficient for the total scattering cross-section, 18 and how Verkruysse et al. rescaled the absorption and scattering coefficients, 17 we redefined μ t as a weighted average of the interaction strengths of the two materials. Thus, the new μ t , μ t;total , has the form μ t;total ¼ xμ t;1 þ yμ t;2 μ t;total ¼ xμ t;1 þ ð1 − xÞμ t;2 ; (9) and is used in place of μ t in Eq. (1).

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 α 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). The linear relation for α 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 α total no longer holds as it did in the layered model.
In terms of the bowling-ball-ping-pong-ball analogy discussed in Sec. 3.3, a material with the larger μ t , or larger interaction cross-section, would be represented by the bowling ball, whereas the one with the lower μ t , or smaller interaction crosssection, 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, α total would not prorate as 50% α bowling and 50% α ping-pong but rather have a higher weighting of α bowling than  α ping-pong . As a result, the α of the material with the higher μ t , in this case, the bowling ball, would dominate the curve for α total , an argument that is also mentioned, but not fully characterized, by Schmitt and Kumar. 18 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," μ t;1 > μ t;2 and α 1 > α 2 .
In another situation, we could have μ t;1 > μ t;2 while α 1 < α 2 . Since α of a given material is dependent on the anisotropic factor, g, and μ 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 shotput balls (material 2). Since the plastic balls are larger, they have a higher interaction cross-section (μ t;1 > μ 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 (α 1 < α 2 ). Since μ t;1 > μ t;2 , α 1 will dominate the curve for α total . The model with a mixture of blood and epidermis is comparable to this analogy (Fig. 5).
The direction in which the curve for α total bends is determined by comparing μ t of the two materials. The material with the higher μ t will have a greater weighting in the total attenuation coefficient, α total . Thus, if the material with the greater μ t also has a greater α, then the curve for α total will bend above the linear relationship (Fig. 4). If the material with the greater μ t has a smaller α, then the curve will bend below the linear relationship (Fig. 5). How much the curve for α total deviates from the linear relationship depends on how much the values of μ t for the two materials differ. In the case of the water-blood mix, μ t;blood ≫ μ t;water (see Table 1), thus, the curve for α total bends far from the linear relationship (Fig. 4). On the other hand, in the blood-epidermis mixed model, μ t;blood > μ t;epidermis , but they are of the same order of magnitude ( Table 1). As a result, although the curve is more weighted to α 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. 17 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 μ a or μ t in the absorption factor (μ a ∕μ 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 μ 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 μ s of epidermis, and in turn, lowers the μ t of epidermis as well [Eqs. (2) and (7)]. Since μ t of epidermis is already lower than that of blood (Table 1), the difference between the μ t s of the two materials becomes even greater, and thus, we should see the nonlinearities become more pronounced based on previous arguments.

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 μ t values of the materials differ by an order of magnitude or more; otherwise, a layered model may be sufficient.