1. Introduction

From the point of view of biomedical optics, light may interact with biological tissues by absorption, scattering, reflection/refraction, quasi-elastic scattering (e.g., laser-Doppler effect), etc. Among these types of interaction, absorption has a particular important role to play for biological applications. In fact, depending on the intensity of the utilized light sources and on the tissue composition (e.g., a particular chromophores content), absorption interaction may manifest itself in different ways related to different biomedical applications. In this sense, we can cite for example photochemical, photothermal or photomechanical interactions [1].

In practice, photochemical interactions allow one to locally activate specialized molecules, previously injected in the body/tissue, permitting, e.g., to treat tumoral tissues (see, e.g., Ref. [2] and similar annual conferences).

Photothermal interactions appear when the light intensity is sufficiently high to increase tissue temperature (i.e., local blood flow is not high enough to cool the tissue and maintain it to a normal temperature). Therefore, targeted pathological tissues are treated by “burning them out”. Note that, photothermal light-tissue interactions may, e.g., also be exploited for tissue cutting, coagulating or welding [3].

Concerning photomechanical interaction, a clever choice of light intensity and duration allows one to induce, thanks to a moderate local increase in temperature, a tissue “dilation” (pressure increase). This sudden increase in pressure produces sound waves that can be recorded and exploited to reconstruct high quality 3D images of the tissues vessels (e.g., arterioles, capillaries) [4].

Other examples might be given, but the main aim here is to highlight the fact that all the above phenomena depend from one common optical quantity: the value of the fluence rate, $\Phi (\textbf {r})$, at position $\textbf {r}$ in the medium. In fact, in an intuitive manner and without to go into details, we can recall that photochemical interactions can be described, e.g., with the following model [1]:

From Eqs. (1), (2) and (3) it clearly appears that all these phenomena depend on $\Phi (\textbf {r})$ and that the quantities of interest are proportional to $\Phi (\textbf {r})$. The above few examples also highlight the reason why $\Phi (\textbf {r})$ represents a fundamental quantity in biomedical optics, and why great effort has been made in the past to simulate $\Phi (\textbf {r})$, e.g., with Monte Carlo (MC) models [5].

Since its creation until nowadays, MC modelization of $\Phi (\textbf {r})$ has mainly been performed by using the method originally conceived by Wang et al. [5] and largely employed in biomedical optics [6–11]. As it is well known, this method is based on the determination of the fraction of absorbed light power in a small tissue volume. In this work we will refer to this method as “classical method”. However, this method is not the only possible one allowing to assess $\Phi (\textbf {r})$. For this reason, and considering the importance of $\Phi (\textbf {r})$ in biomedical optics, we would like to propose here an alternative approach based on photons pathlengths instead of absorbed power.

To this aim, in the present contribution we will exploit a method currently utilized in nuclear physics [12–16] but, to the best of our knowledge, never adopted in near infrared spectroscopy (NIRS) which is the focus of this work. The definition of fluence provided by the International Commission on Radiological Units and Measurements (ICRU) since 1971 is based on the number of the particles crossing an infinitesimally small sphere (for a recent glossary see Ref. [17]). The works by Chilton, Carlsson, and Papiez and Battista [13–16] were aimed to overcome the weaknesses of the ICRU definition, such as for instance the specific shape of the sampling volume, proposing a second definition of fluence in terms of the average value of the pathlength spent by particles within any infinitesimal volume on the macroscopic scale, i.e., as the pathlength per unit volume. Thus, in the following section we will demonstrate, within the frame of the steady state radiative transfer equation (RTE), the link existing between $\Phi (\textbf {r})$ and the length of the random paths [18] traveled by emitted photons in a medium. Beyond its theoretical interest, this method can also be applied as an alternative approach to assess $\Phi (\textbf {r})$ in MC simulations. For the sake of completeness the classical method will also be derived. For simplicity we will call the pathlength based method, PBM, and the classical absorption based method, ABM.

2. Theory

In the present section we will analytically derive the relationship existing between $\Phi (\textbf {r})$ and the photons’ pathlengths. The classical relation linking $\Phi (\textbf {r})$ and the absorbed light power is also derived.

2.1 Calculation of the average fluence in a sub volume with non-uniform absorption coefficient: PBM

Let’s define a medium of volume $V$ with absorption coefficient $\mu _{a0}(\textbf {r})$ and scattering coefficient $\mu _{s0}(\textbf {r})$ [19] with $\textbf {r} \in V$ and external surface $\Sigma$. The medium is illuminated by both internal and external light sources emitting a total power $\mathcal {P}_e$. Let also be a subvolume $V_i \subseteq V$ and $g_i(\ell _i;\mu _{s0}(\textbf {r}), \mu _{a0}(\textbf {r}), p(\theta,\varphi ), n_{in}(\textbf {r}),n_{out}(\textbf {r}_\Sigma ),{\cal G}_{medium},{\cal G}_{sources})$ the probability density function (pdf) for an emitted photon to travel the pathlength $\ell _i\ge 0$ inside $V_i$, before being absorbed or exiting the medium. By definition: $\int _0^{+\infty }g_i(\ell _i;.) d\ell _i=1$. The pdf $g_i(\ell _i;.)$ depends on the distribution of the optical properties of the medium $\mu _{s0}(\textbf {r})$, $\mu _{a0}(\textbf {r})$, $p(\theta,\varphi )$, $n_{in}(\textbf {r})$ and $n_{out}(\textbf {r}_{\Sigma })$, i.e., for $\textbf {r} \in V$, the scattering coefficient, the absorption coefficient, the phase function, the refractive index of the medium and the refractive index on the external surface, $\textbf {r}_{\Sigma } \in \Sigma$, respectively. The symbols ${\cal G}_{medium}$ and ${\cal G}_{sources}$ formally represent the dependence from the geometry of the medium and from the distribution of the light sources, respectively. For the sake of clarity, in the remaining text we will use the simplified notation

Note that $g_i(\ell _i; \mu _{a0}(\textbf {r}))$ must also take into account photons that never entered the volume $V$ due to reflections at the external boundary. These photons obviously have $\ell _i=0$ inside $V_i$, but must be taken into account when assessing the mean pathlength of all emitted photons.

The pdf $g_i(\ell _i; \mu _{a0}(\textbf {r}))$ multiplied by $d\ell _i$ represents the probability that a photon travels a length between $\ell _i$ and $\ell _i+ d\ell _i$ inside $V_i$. Thus, $\mathcal {P}_e g_i(\ell _i; \mu _{a0}(\textbf {r}))d\ell _i$ is the fraction of the emitted radiation $\mathcal {P}_e$ that travels a length between $\ell _i$ and $\ell _i+ d\ell _i$ inside $V_i$. It follows that the magnitude of the emitted radiation that travels any length $\ell _i$ in $V_i$ (including $\ell _i=0$), is

Part of the emitted power, $\mathcal {P}_e$, will be absorbed in the volume $V_i$.

According to the RTE the total absorbed power inside the volume $V_i$ can be expressed as [19]

The same formula is valid for the total power absorbed $\mathcal {P}_{A_{}}$ inside the whole volume $V$ by integrating over $V$.

Let’s now express the increase in absorbed power in $V_i$, $\delta \mathcal {P}_{A_{i}}(\delta \mu _{a_i})$, due to an infinitesimal and constant increase, $\delta \mu _{a_i}$, of the absorption inside $V_i$. Due to the absorption increase the total absorbed power inside $V_i$ becomes

In the limit $\delta \mu _{a_i} \rightarrow 0$ we have that the significant terms of $\delta \mathcal {P}_{A_{i}}$ are

It can be noted that there are two contributions to the variation of the total absorbed power, one related to the absorption variation $\delta \mu _{a_i}$ and a second one related to the variation of fluence $\delta \Phi$ generated by the absorption variation. Thus, we have a variation of absorbed power between a baseline (absorption $\mu _{a0}(\textbf {r})$) and a final state (absorption $\mu _{a0}(\textbf {r})+\delta \mu _{a_i}$). While the first contribution pertains only to the region $V_i$, the second term $\delta \Phi$ is related to absorption variations due to changes in the fluence inside the whole medium and not only in $V_i$. It can be finally noted that the same formula [Eq. (9)] can be used to express the total variation of absorbed power $\delta \mathcal {P}_{A_{}}$ in $V$ by replacing $V_i$ with $V$.

The first term of Eq. (9), for the properties of the microscopic Beer–Lambert–Bouguer’s law (mBLBL) (valid within the RTE; see Appendix A), can be expressed as

We note that this term is the change of absorbed power inside $V_i$ due to the change in the absorption coefficient when the fluence is fixed at the baseline level $\Phi (\textbf {r};\mu _{a0}(\textbf {r}))$.

By using Eq. (10) and by expanding the exponential in Taylor series, we obtain

It can be noted that

By equating the right side of Eq. (13) with the right side of Eq. (11), in the limit $\delta \mu _{a_i} \rightarrow 0$ we obtain

We use the obtained relation to express the fluence in terms of the mean pathlength

Note that, beyond the present purposes, Eq. (16) can also be used to analytically calculate $\langle \ell _i \rangle$ when an analytical solution of the fluence rate inside $V_i$ is known. This can be for instance the case of many homogeneous geometries for which analytical solutions of the diffusion equation for the fluence are available.

2.2 Calculation of the average fluence in a sub volume with uniform absorption coefficient: classical theory, ABM

We derive here the classical theory (utilized in MC simulations) by showing that it can be obtained as a particular case of Sec. 2.1, when $\mu _{a}(\textbf {r})$ is constant over $V_i$. In fact, in the case $\mu _{a}(\textbf {r}) = \mu _{a_i}$ for $\textbf {r} \in V_i$, i.e., $\mu _{a_i}$ is a constant, Eq. (6) can be rewritten as

By introducing Eq. (17) in Eq. (16) we also obtain

Again we note that $\langle \ell _i \rangle$ appearing in Secs. 2.1 and 2.2 represents the mean pathlength traveled in $V_i$ by all the photons emitted by the light sources (thus also, e.g., truncated paths due to absorption inside $V_i$).

Equation (18) also offers an alternative way to calculate the absorbed power inside a volume $V_i$ by knowing the absorption coefficient of the volume, the emitted power and the mean path length $\langle \ell _i \rangle$.

In summary, in this section we have demonstrated how the mean fluence rate $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$, in a volume $V_i \subseteq V$, can be derived from the mean pathlength $\langle \ell _i \rangle$ of all emitted photons [Eq. (16)]. The $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$ method represents an alternative to the classical $\langle \Phi _i \rangle _{\mathcal {P}_{A_i}}$ method [Eq. (17)], with the advantage that $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$ is also applicable in the case $\mu _{a}(\textbf {r})$ is not constant for $\textbf {r} \in V_i$.

3. Results

3.1 Analytical example 1: non-absorbing medium subject to uniform Lambertian illumination

In this example we would like to highlight the interesting fact that the PBM also works in the particular case where $\mu _{a_i}=0$. For example, let be a generic finite medium of volume $V$ with $\mu _{a}(\textbf {r})=0$, $\forall \textbf {r} \in V$ (and thus $\mu _{a_i}=0$) and non nil scattering coefficient. The medium $V$ has refractive index $n_{in}$, and the medium outside $V$ refractive index $n_{out}$. For the present example let chose $V_i=V$. The surface of the medium is illuminated by a uniform and isotropic light source (Lambertian illumination) of total power $\mathcal {P}_e$. Thus, thanks to the invariance property (IP) [20–23], we immediately know that [21,22]

3.2 Analytical example 2: infinite medium with uniform absorption coefficient and a point source

Let’s derive now $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$ for an homogeneous infinite medium and a point light source situated at the axis origin, with constant non nil absorption, $\mu _a$, and scattering, $\mu _s(\textbf {r})$, coefficients. $V_i$ will represent the whole medium. Thus, from Eqs. (16) and (14)

From the constant $\mu _a$ and Eq. (6) (for the case $\delta \mu _{a_i}=0$), the exact solution of the above integral (with $V_i$ infinite) is

Equation (24) represents the mean pathlength traveled by all the photons starting from the source to the point where they are (“forcefully”) absorbed. In fact, if we take the limit $\mu _a \rightarrow 0$ then $\langle \ell _i \rangle$ becomes infinite, i.e., the photons are never stopped by absorption. In this case, it is easy to check that the assessment of $\langle \ell _i \rangle$ [Eq. (24)] by means of the fluence rate is exact. It is sufficient to note that in this particular example it is not important to consider the fact that photons undergo scattering events. The essential point is to know the pdf, $g(\ell _i)=\mu _a\exp (-\mu _a \ell _i)$, that a photon is absorbed after a pathlength $\ell _i$. This allows us to obtain the mean pathlength by means of a trivial calculation, i.e., $\langle \ell _i \rangle = \int _0^{+\infty }\ell _i g(\ell _i)d\ell _i=\mu _a^{-1}$; confirming Eq. (24). Note that these results hold for any type of source distribution and phase scattering function. It can be finally noticed that the above expression for the pathlength traveled by all the photons is different from the average pathlength traveled by the radiation received by a point receiver measuring the fluence at a distance $r$ from the source that is the quantity usually employed in diffuse optics [19,26].

3.3 Numerical example (MC): 100-layered laterally infinite slab illuminated by a plane wave

In this example we will assess, by MC simulations, the mean fluence rate with the ABM and the PBM. This will allow us to highlight some characteristics related to the numerical convergence of these two types of calculation.

 

Fig. 1. Summary of the mean fluence values, computed with the PBM [$\langle \Phi _i \rangle _{\langle \ell _i \rangle }$] and ABM [$\langle \Phi _i \rangle _{\mathcal {P}_{A_i}}$] methods, for each of the 100 layers and 18 different combinations of $\mu _a$ and $\mu _s$ (see text) for a total of 1800 points.

Download Full Size | PDF

We will treat the case of a laterally infinite 100-layered slab. For this, let we start with a laterally finite 100-layered slab of 10 mm thickness having top surface $S$ and composed by 100 layers of $L=0.1$ mm thickness and volume $V_i=S L$ each. All the layers have the same optical parameters $\mu _a$, $\mu _s$, refractive index $n=1.4$ and isotropic phase function. The external medium has refractive index $n_{out}=1$. An external plane wave of unitary photon flux normally impinging to the slab top surface was utilized as a light source (i.e., an impinging flux vector of intensity $I_e=\frac {\mathcal {P}_e}{S}=1$ W mm$^{-2}$, for more details see Refs. [21,22]). This is equivalent to have $\mathcal {P}_e/S = 1$ W mm$^{-2}$ (constant ratio even when $S$ goes to infinity). Thus, for the chosen light source and geometry we can rewrite Eq. (16) for any given layer “$i$” $(i=1,2,\ldots 100)$ as:

In Fig. 1 we have reported $\langle \Phi _i \rangle _{\mathcal {P}_{A_i}}$ as a function of $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$, obtained in this example. As expected, we observe a linear relationship, where the observed random distance from the identity line is due to the statistical noise intrinsic in the MC method. Obviously, if we choose $N\rightarrow +\infty$, all the points will eventually fall on the identity line.

To better evaluate the reliability of the PBM, let’s compare it to the classical ABM by defining

 

Fig. 2. Consistency of the $\langle \Phi _i \rangle _{\mathcal {P}_{A_i}}$ and $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$ MC data as a function of the layer number # (see text). The range for the abscissa (ordinate) axis is the same for all the panels.

Download Full Size | PDF

From Fig. 1 it is not possible to determine if the distance of the data from the identity line is mainly due to a statistical variability of the ABM or PBM or both. Moreover, no information appears on the role of different $\mu _s$ values. To better clarify this point in Fig. 3 we have reported the relative error on the fluence, $RE_{\langle \ell _i \rangle }$, for the PBM method, i.e.,

To better appreciate the convergence property of the two methods versus $N$ in Fig. 4 the relative error on the fluence is plotted for $N=10^6$. Comparing Fig. 4 ($N=10^6$) with Fig. 3 ($N=10^8$) we notice a reduction of the relative errors on both $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$ and $\langle \Phi _i \rangle _{\mathcal {P}_{A_i}}$ by a factor $\approx 10$ when $N$ is changed from $10^6$ to $10^8$. This behavior is an indicator of the convergence speed of the two methods versus $N$ and it shows that increasing the number of injected photons ABM and PBM have a similar improvement in the convergence of the calculation.

 

Fig. 3. Relative errors $RE_{\mathcal {P}_{A_i}}$ and $RE_{\langle \ell _i \rangle }$, on the calculated fluence $\langle \Phi _i \rangle _{\mathcal {P}_{A_i}}$ and $\langle \Phi _i \rangle _{\langle \ell _i \rangle }$, computed with the ABM method and the PBM method, respectively, as a function of the layer number #. The range for the abscissa (ordinate) axis is the same for all the panels and the results pertain to $N=10^8$.

Download Full Size | PDF

The choice to report the relative error in the comparison between ABM method and PBM method is motivated by the fact that this quantity characterizes the reliability of the simulated data [27]. In general, a MC result should show a relative error lower than 0.1 to produce reliable confidence intervals of the simulated data [27]. About the reliability of the presented results we also notice that the inner core of the code has been accurately verified by previous investigations [25,28].

Finally, we note that, although the presented MC results were obtained by using a Henyey-Greenstain scattering function with asymmetry factor $g=0$, similar results (i.e., the faster convergence of the PBM method) are expected for any other kind of scattering function and distribution of the optical properties ($\mu _a$, $\mu _s$).

 

Fig. 4. The same as Fig. 3 but for $N=10^6$.

Download Full Size | PDF

4. Conclusions

In conclusion, we have analytically and numerically derived a useful formula, i.e., Eq. (16), valid in complete generality within the steady state RTE, to obtain the mean fluence rate, in a given sub-volume of a propagating medium, with a method (PBM) alternative to the classical one (ABM) usually utilized in NIRS. This method is based on the calculation of the mean pathlength traveled by all emitted photons inside a sub-volume of any shape and size. As it can be deduced from the theory and the presented examples, the PBM can be easily implemented in any MC simulation for photon transport.

A potential advantage of the PBM (as shown, e.g., in Sec. 3.3) is that it might improve the statistics and thus the convergence of the MC simulations. In fact, when measuring the mean pathlength (PBM), all the photons (paths) that cross the region of interest are taken into account. While, with the classical approach (ABM), only the “absorbed” photons are considered. In the latter case, this may produce (especially for lower values of the absorption coefficient) a poorer MC statistics for the same number of launched photons. However, it remains to verify if the use of other MC methods for photon propagation (Microscopic Beer-Lambert law, Albedo-weight, Absorption and Scattering pathlength rejection method [29]) and especially their GPU implementations, will still confirm these results. This interesting observation will certainly be a matter of future investigations. Further, the PBM method can be also applied when the absorption coefficient is not constant or nil inside the considered sub-volume, whilst the applicability of the ABM method is not possible under these conditions. The reliability of the presented results has been also accurately investigated. For example, one crucial point of a MC simulation is to deal properly with the reflections at the internal and external boundaries. This was done in our previous works where the MC code was tested against the analytical solution of steady-state RTE for layered media illuminated by a uniform Lambertian source [22,25,28].

In comparison with the previous literature, we have fully re-framed the $\langle \ell _i \rangle$-$\Phi (.;.)$ relationship within the RTE. In practice, the derivation proposed in this paper offers a more general view of the possible applications of the PBM method in biomedical optics for those cases where the information on the fluence is relevant.

The definition of fluence given by ICRU [17] is different from that currently used within the RTE [19]; however, it is worth to note that even in the ICRU report is mentioned the relationship $\langle \ell _i \rangle$-$\Phi (.;.)$ that in the present work is obtained within the RTE. This to emphasize that the awareness of this relation is present in many fields such as dosimetric calculations where “the fluence is frequently expressed in terms of the lengths of the particle trajectories” even when only an intuitive proof of its validity is given.

Thus, the PBM method appears to be based on a fundamental law that maybe should be better considered in biomedical optics. In this sense, we hope that the present contribution will serve to fill a missing part on this interesting topic, and that it will also serve as a stimulus to develop original MC investigations related to biological tissues.

We want to finally stress that the obtained relation [Eq. (16)] can be used to analytically calculate $\langle \ell _i \rangle$ when an analytical solution of the fluence rate inside $V_i$ is known (e.g., diffusive media). Although $\langle \ell _i \rangle$ is the mean pathlength of all emitted photons and not the mean pathlength of the received photons usually employed in biomedical optics, it can return important information for the understanding of photon migration through diffusive media such as biological tissues. It can be at last expected that the validity of the proposed relation between fluence and mean pathlength should be also verified in a regime of anomalous transport [30].