## Abstract

We present an efficient Monte Carlo algorithm for simulation of time-resolved fluorescence in a layered turbid medium. It is based on the propagation of excitation and fluorescence photon bundles and the assumption of equal reduced scattering coefficients at the excitation and emission wavelengths. In addition to distributions of times of arrival of fluorescence photons at the detector, 3-D spatial generation probabilities were calculated. The algorithm was validated by comparison with the analytical solution of the diffusion equation for time-resolved fluorescence from a homogeneous semi-infinite turbid medium. It was applied to a two-layered model mimicking intra- and extracerebral compartments of the adult human head.

©2008 Optical Society of America

## 1. Introduction

Fluorescence methods are valuable tools in medical diagnostics [1]. Often the concentration of the fluorophore embedded in the tissue is of primary interest, whereas its fluorescence spectrum, lifetime or quantum yield may reveal information about its molecular environment. So far, fluorescence-based *in vivo* imaging of human tissue was mainly restricted to superficially accessible tissues like skin [2] and the mucosa of hollow organs [3]. On the other hand, in *vivo* imaging of fluorescence-labeled molecular probes was reported for small animals [4, 5] demonstrating the potential to image tumorigenesis, host-tumor interactions and specific effects of pharmacological interventions [6].

A number of studies focused on the evaluation of fluorophore concentration in a turbid homogeneous medium were based on the application of continuous (cw) light [7, 8]. Hyde et al. [9] and Stasic et. al. [10] derived depth-dependent fluorophore concentrations from spatially-resolved cw measurements, but the application of this method was limited to short interoptode distances (up to 10 mm) and hence to superficial tissue compartments. For tomographic reconstruction of the concentration of the dye distributed non-homogeneously in specific regions of the body frequency domain [11] and time-domain [12] approaches were suggested. In Refs. [13–16] fluorescence lifetime or quantum yield were analyzed as well.

Studies of fluorescence light remitted from deep within large tissue volumes were carried out up to now on canine breast [17] and human breast [18]. Recently, we reported on exciting and recording fluorescence of a dye bolus non-invasively through the skull and scalp of the adult human head, providing evidence that the fluorescence observed originated from the cortex [19].

Monte Carlo (MC) simulations are the method of choice to model fluorescence signals from inhomogeneous turbid media. The classical fluorescence MC algorithm reported by Welch et al. [20, 21] is based on the evaluation of the probability of absorption of a photon in each scattering event. Crilly et al. and Swartling et al. [22, 23] proposed accelerated algorithms based on two independent MC runs. The first run is performed for the optical properties of the medium and fluorophore at the excitation wavelength and yields the spatial distribution of photons in the medium. In the second run a similar simulation is performed at the emission wavelength for the source of photons located at the position of the detector (“reverse emission”). Finally, the algorithm requires the convolution of the results of these two runs. The accuracy of the results provided by this algorithm depends on the size of the spatial and temporal grids used for storing the results of these two simulations.

In the present paper we propose a new efficient algorithm for fluorescence Monte Carlo simulations in layered media. It enables simulation of propagation of excitation as well as fluorescence photons in the tissue. A preliminary version of this code was already successfully used to explain the results reported in [19]. Here the algorithm is discussed and validated in more detail. In addition to the simulation of distributions of times of arrival of fluorescence photons, we derive the spatial distribution of fluorescence generation inside the medium. Our MC algorithm is validated in a series of simulations carried out for fluorescence from a homogeneous semi-infinite medium by comparison with the corresponding solution of the diffusion equation. The algorithm is applied to simulate fluorescence signals from a dye bolus in the human head.

## 2. Monte Carlo simulations

#### 2.1 Theoretical background

In this section we describe our Monte-Carlo approach to simulate the generation and propagation of fluorescence photons in a layered tissue model. The geometry of the model is outlined in Fig. (1). A point source is located at the plane surface of the turbid medium. The photons are collected at a source-detector separation *r*. The absorption and reduced scattering coefficients of the layered medium can be defined independently for each layer, indexed by *j*. For modeling absorption, we applied a variance reduction technique, i.e. we follow the propagation of a photon bundle of gradually decreasing weight at the excitation wavelength *λ _{x}* [24, 25]. In a similar way, the concept of photon bundles was applied to fluorescence. Conversion of a fraction of the absorbed photons into fluorescence leads to the generation of fluorescence photon bundles (at emission wavelength

*λ*) the weights of which decrease due to absorption at

_{m}*λ*.

_{m}For calculating the trajectories of the photon bundles we started from the following assumptions:

- Scattering is isotropic and can be characterized by a single reduced scattering coefficient. For interoptode distances of a few centimeters as considered here, a completely randomized scattering direction is a good approximation.
- The reduced scattering coefficients of each layer
*j*are equal at the excitation (*λ*) and emission (_{x}*λ*) wavelengths, i.e._{m}*µ*′_{sx,j}≈*µ*′_{sm,j}≈*µ*′_{s,j}. This assumption is justified since excitation and emission wavelengths usually differ by a few 10 nm only. In the near-infrared region the reduced scattering coefficient of tissues typically depends only slightly on wavelength, i.e. for such small differences in wavelength the reduced scattering coefficient differs by less than 10% [26]. - All fluorescence photon bundles generated along the path of a particular bundle of excitation photons follow this very trajectory. This assumption is justified since the ensembles of all excitation and fluorescence trajectories from a given location within the medium to the detector are virtually the same provided scattering is isotropic and equal at
*λ*and_{x}*λ*._{m}

First the trajectory is calculated without considering any absorption. For each scattering event, the path length to the next scattering event is determined by

where *R* is a random number within the interval (0,1). When the path between two subsequent scattering events crosses the boundary between two layers with *j*=*A* and *j*=*B* respectively, the portion of the free path length in medium *B* is recalculated as

where *L _{B}* is the corresponding value before considering different scattering coefficients. Fresnel reflections at the boundaries between the layers as well as at the surface of the tissue are neglected to avoid prolongation of the computation time. After each scattering event

*q*the cumulative path lengths

*l*traveled by the photon bundle

^{(q)}_{ij}*i*within the layers

*j*up to this scattering event

*q*are stored for further processing after completion of the trajectory calculation. At the point of detection

*e*, i.e. upon crossing the surface of the tissue, the total path lengths for this bundle in all layers

*l*are evaluated. Based on

^{(e)}_{ij}*l*and

^{(q)}_{ij}*l*, multiple sets of absorption parameters can be applied to the same trajectory as will be discussed below.

^{(e)}_{ij}The intrinsic absorption coefficients of the *j*-th tissue layer are denoted by *µ _{ax,j}* and

*µ*at the excitation and emission wavelengths, respectively. The fluorescent dye (fluorophore) that is distributed in the tissue contributes to absorption which is described by the concentrationdependent absorption coefficients

_{am,j}*µ*and

_{afx,j}*µ*at

_{afm,j}*λ*and

_{x}*λ*, respectively. The coefficient

_{m}*µ*describes reabsorption of fluorescence photons by the fluorophore. It should be noted, however, that a possible reemission process after reabsorption is omitted. This assumption is justified as long as the fluorescence quantum yield of the fluorophore is small which is true for a dye like indocyanine green (ICG). Moreover, considering the coefficients

_{afm,j}*µ*and

_{am,j}*µ*at a fixed

_{afm,j}*λ*only, we neglect that fluorescence photons in different parts of the fluorescence spectrum might experience different reabsorption by the fluorophore as well as absorption by the tissue.

_{m}Bundles of excitation photons of unity initial weight are launched from a point source at the surface of the tissue model. The conversion probability *P*
^{(q)}
_{c,i} along the path between two consecutive scattering events numbered by *q*-1 and *q* is given by the probability of absorption by the fluorophore and its fluorescence quantum yield *Φ*.

where we assume *Φ _{j}*=

*Φ*for all layers, i.e. the same microenvironment in all compartments. For

*q*=1, i.e. on the path between source and first scattering event,

*l*

^{(q-1)}

_{ij}equals 0.

Up to scattering event *q*, the weight of excitation bundle *i* has decreased to

by the absorption of the tissue as well as of the dye at *λ _{x}*. The initial weight of the fluorescence photon bundle emanating from scattering event

*q*is thus

*W*

^{(q)}

_{x,i}

*P*

^{(q)}

_{c,i}. On its travel to the surface of the tissue it decreases due to absorption at

*λ*by the factor

_{m}After summation over all scattering events, the total contribution of the excitation photon bundle *i* to the fluorescence signal is

Such contributions are sorted according to the distance *r _{i}* between the source position and the point of exit at which the bundle

*i*escapes from the tissue as well as according to its total travel time

*t*. Fluorescence histograms, i.e. distributions of times of arrival of fluorescence photons (DTA), are collected for various detection areas. The total travel time includes the time traveled by the photon bundle at the excitation wavelength up to the point of conversion plus the time traveled by the fluorescence photon bundle from the point of conversion to the point of detection (for simplicity the fluorescence lifetime has been neglected as discussed below). Assuming the same velocity of light

_{i}*c*=

*c*/

_{0}*n*(

*c*=3·10

_{0}^{8}m/s,

*n*- refractive index of the tissue) in all layers and at

*λ*as well as

_{x}*λ*, the total travel time from source to detector of both the excitation and the fluorescence photon bundle

_{m}*i*is

The fluorescence lifetime *τ _{F}* has been neglected so far. It would lead to an exponentially distributed delay in the arrival of the fluorescence photon bundles at the detector. Consequently, a finite fluorescence lifetime that is the same in all layers can be easily taken into account by convolving the DTA obtained as described before with exp(-

*t*/

*τ*). Hence, when considering statistical moments of the DTA,

_{F}*τ*simply adds to the first moment of the DTA and

_{F}*τ*

^{2}

_{F}to its variance.

So far, the weights given in Eqs. (4) and (5) were calculated assuming that the conversion of excitation into fluorescence photons occurs exactly in the scattering point (indexed by *q*). In fact, the conversion may appear at any location in between the scattering points *q-1* and *q*. In case of different total absorption coefficients at *λ _{x}* and

*λ*, such change in the location of conversion leads to a slight modification of the weight of the fluorescence photon bundle, i.e.

_{m}instead of Eq. (6), where

and *l*
^{(q)*}
_{ij} is the distance traveled in *j*th layer by *i*th photon bundle between the actual conversion point and *q*-th scattering event: 0<*l*
^{(q)*}
_{ij}<(*l*
^{(q)}
_{ij}-*l*
^{(q-1)}
_{ij}). The exponential term in Eq. (8) is close to unity provided that
$\sum _{j}\frac{\Delta {\mu}_{a,j}}{{\mu}_{s,j}\prime}$
is small. The latter is true for typical dye concentrations and tissue optical properties, since the absorption coefficients of the tissues of the head as well as of a dye at a concentration of several 0.1 µmol/l are on the order of 0.1 cm^{-1} at most, and *µ′ _{s}* is on the order of 10 cm

^{-1}. Therefore we decided to calculate

*W*according to Eq. (6).

_{m,i}In parallel to the calculation of contributions to DTAs of fluorescence photons the contribution to the distribution of times of flight (DTOF) of diffusely reflected photons is analyzed for each excitation photon bundle. This contribution is given by the survival weight of the ith photon bundle at *λ _{x}* when arriving at the detector

which depends on the absorption coefficients of the tissue *µ*
_{ax,j} and the dye *µ*
_{afx,j} and the total path lengths in all layers. These weights are again accumulated according to *r* and *t*.

For the simulation of complete sets of DTAs and DTOFs taking into account varying absorption coefficients, we applied the “white Monte Carlo” [27] approach. According to Eqs. (1–4) and (10) the weights of the excitation and fluorescence photons bundles as well as the probabilities of conversion to fluorescence can be evaluated when the simulation of the trajectory of bundle *i* is completed and the values of *l*
^{(q)}
* _{ij}* and

*l*

^{(e)}

_{ij}are known. Thus multiple sets of absorption coefficients may be applied to the same trajectory. It should be noted that the multiple DTAs and DTOFs obtained in this way are not statistically independent. However, this is not relevant if the number of photon bundles in the simulation is large enough to ensure low statistical uncertainty of the quantities derived. The “white Monte Carlo” approach allows for a formidable saving of computation time. It is applicable if the geometry, i.e. the thickness of layers, as well as the scattering coefficients

*µ*′

_{s, j}in each layer remain constant.

To analyze the spatial profile of fluorescence generation at a fixed source-detector separation, we consider the contributions to the fluorescence signal that are generated within a given voxel *V* of size *Δx*, *Δy*, *Δz* centered around (*x,y,z*) by all photon bundles which pass this voxel and have a total travel time within the interval (*t*-Δ*t*/2, *t*+Δ*t*/2)

with

$q:{x}_{q}\in (x-\frac{\Delta x}{2},x+\frac{\Delta x}{2}),{y}_{q}\in \left(y-\frac{\Delta y}{2},y+\frac{\Delta y}{2}\right),{z}_{p}\in \left(z-\frac{\Delta z}{2},z+\frac{\Delta z}{2}\right)$

where *x _{q}*,

*y*and

_{q}*z*are coordinates of the

_{q}*q*th scattering event. The first summation extends over all scattering events

*q*that occur within the given voxel. We call the spatial distribution defined by Eq. (11) “generation probability distribution”. The total spatial integral of the expression in Eq. (11) as a function of time equals the DTA for a given source-detector separation. By restricting the integration to certain compartments (layers) of the medium, their individual contributions to the DTA of fluorescence photons can be calculated.

#### 2.2 Monte Carlo code

The Monte Carlo code used before for simulation of DTOFs of diffusely reflected light in a layered turbid medium [25, 28] was improved and extended to include the simulation of time-dependent fluorescence according to the algorithm described above. In addition, the option to perform the simulations for a set of absorption coefficients simultaneously was implemented. To further reduce computation time, two procedures were followed, (i) simulation of a photon trajectory was stopped when its total propagation time exceeded the maximum time considered for calculation of DTAs, (ii) photon bundles were collected simultaneously in a set of concentric, consecutive, ring-shaped detectors. The latter is possible because of the rotational symmetry of the layered model around the source. The code allows for free definition of the following parameters: thickness of layers, their optical properties, duration of time windows, radius and width of the ring detectors as well as spatial grid for calculation of the generation probability.

The simulations described below were focused on two cases, (i) a homogeneous semi-infinite medium and (ii) a two-layered medium with an upper layer of 1 cm thickness on top of a semi-infinite, lower compartment. The initial set of optical properties of the tissue and the fluorophore for the homogeneous case is given in Table 1. The simulations were carried out for a total of 150 million photon bundles launched and 10 million photon bundles detected in total at source-detector separations between 5 mm and 62 mm. The distributions (DTAs and DTOFs) were calculated for a set of 19 ring detectors, each 3 mm wide, and for 150 time windows of 25 ps each covering a total time interval from 0 to 3.75 ns. With the optical properties given in Table 1, such simulation took about 10 h on a PC with a Pentium III processor (1.6 GHz). For a set of 100 values of *µ _{afx}*, the computation time increased by a factor of 10 only.

For the calculation of the spatial profiles of fluorescence generation probability according to Eq. (11) it was necessary to use a small detection spot instead of a ring detector. In order to harness all photon bundles detected within a whole ring, the trajectory of each photon bundle was rotated around an axis perpendicular to the surface of the medium and pointing through the source position. After rotation, the *y* coordinate of the escape point was zero. Equation (11) was evaluated with the trajectories after rotation, i.e. with the new coordinates of all scattering events. Furthermore, for the presentation of the spatially dependent generation probability the values of the 3D voxel matrix were projected onto the source-detector plane (x-z plane) perpendicular to the surface of the medium.

## 3. Validation of the Monte Carlo code

For validation of our algorithm, Monte Carlo simulations were compared with simulations based on the analytical expression for the time-resolved fluorescence photon fluence rate [29] as derived from the diffusion uationation for a homogeneous semi-infinite medium and zero boundary conditions

$$\mathrm{where}\phantom{\rule{.2em}{0ex}}{z}_{0}=\frac{1}{{\mu}_{s}^{\prime}},\phantom{\rule{.2em}{0ex}}D\cong \frac{1}{\left(3{\mu}_{s}^{\prime}\right)}\phantom{\rule{.2em}{0ex}}\mathrm{and}\phantom{\rule{.2em}{0ex}}{\beta}_{m}=\left({\mu}_{\mathrm{am}}+{\mu}_{\mathrm{afm}}\right)c,{\beta}_{x}=\left({\mu}_{\mathrm{ax}}+{\mu}_{\mathrm{afx}}\right)c.$$

The comparison of DTAs simulated by the Monte Carlo algorithm and by the diffusion model for the set of optical properties given in Table 1 is presented in Fig. (2). The diffusion data were calculated for values of *r* equal to the mean radius of the corresponding ring. The raw Monte-Carlo DTAs were divided by the area of the corresponding detection rings. Subsequently each of the two data sets of DTAs, i.e. MC and diffusion, was normalized to the sum of the integrals of all DTAs of this set within the time interval from 1 ns to 3.5 ns. As can be seen in Fig. 2, the results of MC simulations match well with the DTAs obtained from the diffusion model with respect to the shape of the DTAs as well as to their relative amplitudes for different source-detector separations.

In another set of simulations the absorption coefficients were varied over a large range. The comparison of the resulting DTAs with corresponding results obtained from Eq. (13) was carried out by analyzing their statistical moments. We considered three moments of the DTAs, zeroth moments, i.e. integral (sum of photons in the DTA, *N _{F,tot}*), first centralized moment, i.e. mean time of arrival of fluorescence photons (<

*t*>

_{F}) and second centralized moment, i.e. variance of DTA (

*V*). The comparison of integrals was performed after normalization to the sum of all values in the set of simulations.

_{F}In four series of simulations the absorption coefficients of the dye and the tissue at the excitation and emission wavelengths were varied in the range between 0.00001mm^{-1} and 0.655 mm^{-1}. For each series of simulations only one absorption coefficient was changed whereas all other coefficients were kept fixed to the values shown in Table 1. Results of the comparison are presented in Figs. (3) and (4). It can be seen that the moments of DTAs of fluorescence photons calculated by the MC method match with the moments calculated from Eq. (13). Discrepancies between both models appear at short interoptode distances and high absorption coefficients. This observation can be explained by the fact that the diffusion approximation approaches its limits in these cases.

## 4. Monte Carlo simulations in a two-layered model

In Fig. (5) spatial distributions of generation probabilities are presented for a source-detector separation of 32 mm. The case of dye distributed homogeneously in the medium with the optical properties as given in Table 1 [Fig. (5a)] was compared with two specific inhomogeneous cases. The dye has a higher concentration (10 times higher *µ _{afx}*) in the lower layer [Fig. (5b)] or in the upper layer [Fig. (5c)], respectively. The thickness of the upper layer mimicking extracerebral tissue was taken as 10 mm. When the dye of high concentration is distributed in the lower layer most of the fluorescence photons are generated in this compartment [s. Fig. (5b)]. In contrast, for the higher dye concentration in the upper layer, fluorescence photons are preferentially generated in the upper layer [s. Fig. (5c)]. This analysis suggests a dependence of mean time of arrival of fluorescence photons <

*t*>

*on the spatial distribution of the fluorophore. In particular, in case of high dye concentration in the lower tissue compartment, <*

_{F}*t*>

*is expected to be larger than if the dye is distributed homogeneously in the whole medium*

_{F}An asymmetry of the probability distributions in respect to the midpoint of the optodes can be observed. This effect is caused by the different absorption coefficients for excitation and fluorescence photons. Because the absorption coefficient of the medium and the dye at the emission wavelength is much lower, the photons that were absorbed and converted to fluorescence close to the source position have a high probability to arrive at the detector.

The dependence of moments of DTAs on the concentration and spatial distribution of the dye was analyzed in more detail [s. Fig. (6)]. Changes in moments of simulated DTAs and DTOFs corresponding to various values of *µ _{afx}* in the lower and upper layers of the model were calculated. Analysis of DTAs shows that <

*t*>

*and*

_{F}*V*reveal specific changes caused by increasing

_{F}*µ*in the lower and upper layers. We observe an initial increase of mean time of arrival and variance of DTAs when

_{afx}*µ*increases in the lower compartment and a decrease of these moments when the absorption by the dye increases in the upper layer. These observations are in agreement with the qualitative analysis of generation probability discussed above. However, for large values of

_{afx}*µ*in the lower layer and large source-detector separations,

_{afx}*V*and <

_{F}*t*>

*decrease again due to saturation effects.*

_{F}In contrast, for diffuse reflectance we observe a decrease in integral, mean time of flight as well as variance of simulated DTOFs with increasing *µ _{afx}* in any of the two layers. These findings are in accordance with our previous results related to the analysis of diffuse reflectance [28, 30, 31].

The next series of Monte Carlo simulations was related to the analysis of our recent *in-vivo* measurements of DTAs after intravenous injection of a bolus of indocyanine green [19]. We considered the temporal changes of moments of DTAs due to inflow and washout of a bolus of a fluorophore by assuming different time courses of changes of *µ _{afx}* in the upper and lower layers to mimic different kinetics of the dye in the extra- and intracerebral tissue compartments. As it was previously shown by perfusion MRI as well as optical measurements the inflow and washout of the dye in cerebral tissue is faster than in the extracerebral compartment consisting of skin and skull [28, 32, 33]. The assumed time courses of

*µ*are shown in the left column of panels in Figs. (7) and (8). It should be noted that we always assumed a constant, homogeneous background dye concentration. In the first row of Fig. 7 we present results obtained for the same changes of

_{afx}*µ*(

_{afx}*T*) in both layers. In this case both mean time of flight and variance of the DTAs decrease after arrival of the bolus. In the second (third) row changes of moments are considered when

*µ*(

_{afx}*T*) increases in the lower (upper) layer only. In case of increase of

*µ*in the upper layer the pattern of changes in the moments is very similar to the homogeneous case (cf. Fig. (7), first and second rows). In contrast, when the change of

_{afx}*µ*appears in the lower layer, an increase of <

_{afx}*t*>

*as well as*

_{F}*V*can be observed (see Fig. (7), third row).

_{F}Furthermore, we modeled the realistic case that the bolus appears in both layers, but its inflow and washout is faster in the lower layer (brain) compared to the upper layer (extracerebral tissue) [s. Fig. (8)]. The assumed time courses of changes of *µ _{afx}*(

*T*) in both layers are shown in the left column of panels in Fig. 8. In addition, we modeled a realistic reabsorption of fluorescence by the dye by assuming that

*µ*(

_{afm}*T*)=

*µ*(

_{afx}*T*)/2. The top row of Fig. (8) shows the results for equal maximum dye concentrations in both layers. We found that the patterns of changes in moments calculated by Monte Carlo simulations are very similar to those observed in

*in-vivo*experiments [19]. After an initial increase of mean time of flight and variance these moments decrease and tend to return to their initial levels corresponding to a low, homogenous dye concentration. Unfortunately, the absolute values of amplitudes of changes of

*µ*are difficult to estimate for the

_{afx}*in-vivo*situation. Therefore, we analyzed how the ratio of amplitudes of changes of

*µ*in the extracerebral and intracerebral compartments influences the observed patterns of moments. Assuming a constant amplitude of the bolus in the lower compartment (brain) we changed the amplitude of the bolus in the upper layer (extracerebral tissue). In general, the observed patterns of changes in each of the moments show a similar behaviour (cf. rows in Fig. 8). Moreover, the patterns of changes of moments observed at both source-detector separations are similar. Specifically, the amplitude of the positive peaks in Δ<

_{afx}*t*>

*as well as Δ*

_{F}*V*seem to be rather invariant for the different extracerebral concentrations and source-detector separations. However, it can be noted that the dynamics of changes of mean time of flight and variance in the transition from positive to negative values is influenced significantly by the ratio of amplitudes of the boli in both layers.

_{F}## 5. Discussion and conclusions

In the present paper we described a novel algorithm for Monte Carlo simulations of fluorescence from layered turbid media, its implementation, its validation for a homogeneous semi-infinite medium and its application for a two-layered model of the adult head. It should be noted that the validation was restricted to the homogeneous case since to our knowledge there is no closed-form solution of the diffusion equation available to describe time-resolved fluorescence in a layered turbid medium.

The major limitation of the algorithm presented here is the assumption of equal scattering coefficients at the excitation and at the emission wavelengths. This assumption is reasonable for living tissues since the change of *µ*′* _{s}* with wavelength is fairly small [26]. Additionally, we assumed that scattering in the medium under investigation is isotropic. This assumption remains reasonable for interoptode distances much larger than the mean free path. On the other hand, these assumptions account for the particular efficiency of the algorithm with respect to computation time since every photon launched into the tissue contributes to the distribution of times of arrival of fluorescence photons. As a result, fast simulation of DTAs is possible even for large interoptode distances. In previously published accelerated algorithms [23], a grid in time and space for simulation of DTAs has to be defined and convolution of histograms is needed. This convolution approach requires the rotational symmetry of the problem. Moreover, the sizes of the grids significantly influence the accuracy of the results of simulations. Our algorithm avoids both limitations.

Similar to other algorithms [27, 34], our code allows to apply the “white Monte Carlo” approach. In this way, fast acquisition of results for various combinations of absorption coefficients in the layers of the model is possible. The time of computations for multiple sets of values of the absorption coefficients of the dye and the tissue at excitation and emission wavelengths is considerably reduced as the most time-consuming part of the code, i.e. the generation of the trajectory of a photon bundle, is executed only once. This option proved particularly advantageous when modeling boli of a fluorophore, i.e. a time course of dye concentration in various layers.

In addition to modeling measurable quantities as, e.g. the DTA of fluorescence photons, we introduced and studied the spatial distribution of fluorescence generation probability. This quantity allowed us to analyze in more detail which parts of the tissue volume contribute to the detected time-resolved fluorescence signal. In this way additional information was gained that is not available from the signals measured on the surface of the tissue. The location of maximum generation probability naturally depends on the fluorophore concentrations in the various layers of the tissue. Moreover, the maximum can be shifted towards source or detector depending on the relative magnitude of *µ _{afx}* and

*µ*, i.e. depending on whether excitation or fluorescence photons are more likely to be absorbed on the trajectory from source to detector. Finally, the spatial distribution of generation probability is time-dependent, i.e. the location of maximum fluorescence generation depends on the total time of flight of the photon bundles.

_{afm}Our efforts to develop an efficient fluorescence Monte-Carlo algorithm were driven by applications of non-invasive *in-vivo* fluorescence detection of ICG boli in the adult human brain, e.g. for perfusion monitoring in stroke patients, and the need to understand the observed signals. We started from a two-layered model simulating intra- and extracerebral tissue and different time courses of boli of the fluorophore in both layers. In this way we could explain the most significant features of the observed in-vivo signals of moments of time-of-flight distributions. The initial increase of mean time of flight and variance is connected with the faster inflow of the dye to the intracerebral compared to the extracerebral compartment. By performing various Monte-Carlo simulations for realistic time courses of dye boli in the brain and overlying tissue, we studied the influence of the ratio of fluorophore concentrations in both compartments as well as of the source-detector separation on the signals detected. By comparing fluorescence signals with signals of diffuse reflectance we found that the magnitudes of changes in moments of DTAs of fluorescence photons are in general larger than the magnitudes of changes of moments of DTOFs of diffusely reflected photons. In addition, mean time of arrival of fluorescence photons <*t*>* _{F}* exhibits more specific dynamic features than mean time of flight of diffusely reflected photons <

*t*>

*. These result matches well with previous*

_{R}*in-vivo*results [19]. Such observations lead us to the conclusion, that monitoring of diffusely remitted fluorescence instead of diffusely reflected photons may be advantageous for the assessment of cerebral perfusion.

So far we followed the simplest possible approach to model the human head, i.e. considering two layers to distinguish between an intra- and extracerebral compartment only. More refined modeling based on a larger number of layers with realistic thickness and optical properties is possible with the existing Monte-Carlo code. The algorithm could be also generalized to include more complex (not only layered) inhomogeneous tissue geometries. Thus it could be useful in simulations of fluorescence-based molecular imaging methods for thick tissues in general. Such methods include the detection of leakage of a blood-pool agent from vessels in neoplastic or inflammatory tissue as well as imaging of the distribution of specific fluorescence-labeled molecular probes in tissue and hence are potentially useful for detection and localization of diseased tissue [35] like in stroke or tumor diagnostics.

## Acknowledgment

The work was supported by the Polish Ministry of Science and Higher Education within project P-N/012/2006 and German Federal Ministry of Education and Research, grants BMBF 01 GZ 0711 (Bilateral German-Polish research projects in neurosciences) and BMBF 01 GO 0208 (Berlin NeuroImaging Center).

## References and Links

**1. **J. Lakowicz, *Principles of Fluorescence Spectroscopy* (Plenum US, 1999).

**2. **K. Svanberg, I. Wang, S. Colleen, I. Idvall, C. Ingvar, R. Rydell, D. Jocham, H. Diddens, S. Bown, G. Gregory, S. Montan, S. Andersson-Engels, and S. Svanberg, “Clinical multi-colour fluorescence imaging of malignant tumours--initial experience,” Acta Radiol. **39**, 2–9 (1998). [PubMed]

**3. **M. Ortner, B. Ebert, E. Hein, K. Zumbusch, D. Nolte, U. Sukowski, J. Weber-Eibel, B. Fleige, M. Dietel, M. Stolte, G. Oberhuber, R. Porschen, B. Klump, H. Hortnagl, H. Lochs, and H. Rinneberg, “Time gated fluorescence spectroscopy in Barrett’s oesophagus,” Gut **52**, 28–33 (2003). [CrossRef]

**4. **V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. **8**, 757–760 (2002). [CrossRef] [PubMed]

**5. **A. Becker, C. Hessenius, K. Licha, B. Ebert, U. Sukowski, W. Semmler, B. Wiedenmann, and C. Grotzinger, “Receptor-targeted optical imaging of tumors with near-infrared fluorescent ligands,” Nat. Biotechnol. **19**, 327–331 (2001). [CrossRef] [PubMed]

**6. **R. Weissleder, “Scaling down imaging: Molecular mapping of cancer in mice,” Nat. Rev. Cancer **2**, 11–18 (2002). [CrossRef] [PubMed]

**7. **R. Weersink, M. S. Patterson, K. Diamond, S. Silver, and N. Padgett, “Noninvasive measurement of fluorophore concentration in turbid media with a simple fluorescence/reflectance ratio technique,” Appl. Opt. **40**, 6389–6395 (2001). [CrossRef]

**8. **K. R. Diamond, T. J. Farrell, and M. S. Patterson, “Measurement of fluorophore concentrations and fluorescence quantum yield in tissue-simulating phantoms using three diffusion models of steady-state spatially resolved fluorescence,” Phys. Med. Biol. **48**, 4135–4149 (2003). [CrossRef]

**9. **D. E. Hyde, T. J. Farrell, M. S. Patterson, and B. C. Wilson, “A diffusion theory model of spatially resolved fluorescence from depth-dependent fluorophore concentrations,” Phys. Med. Biol. **46**, 369–383 (2001). [CrossRef] [PubMed]

**10. **D. Stasic, T. J. Farrell, and M. S. Patterson, “The use of spatially resolved fluorescence and reflectance to determine interface depth in layered fluorophore distributions,” Phys. Med. Biol. **48**, 3459–3474 (2003). [CrossRef] [PubMed]

**11. **E. M. Sevick-Muraca, J. S. Reynolds, J. Lee, D. Hawrysz, A. B. Thompson, R. H. Mayer, R. Roy, and T. L. Troy, “Fluorescence lifetime imaging of tissue volumes using near- infrared frequency domain photon migration,” Photochem. Photobiol. **69**, 66S–66S (1999).

**12. **J. Wu, J. Wang, L. Perelman, I. Itzkan, R. Dasari, and F. Ms, “Time-resolved multichannel imaging of fluorescent objects embedded in turbid media,” Opt. Lett. **20**, 489–491 (1995). [CrossRef] [PubMed]

**13. **R. H. Mayer, J. S. Reynolds, and E. N. Sevick-Muraca, “Measurement of the fluorescence lifetime in scattering media lay frequency-domain photon migration,” Appl. Opt. **38**, 4930–4938 (1999). [CrossRef]

**14. **D. Hattery, V. Chernomordik, M. Loew, I. Gannot, and A. Gandjbakhche, “Analytical solutions for time-resolved fluorescence lifetime imaging in a turbid medium such as tissue,” J. Opt. Soc. Am. A **18**, 1523–1530 (2001). [CrossRef]

**15. **D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. SevickMuraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. **36**, 2260–2272 (1997). [CrossRef] [PubMed]

**16. **C. L. Hutchinson, T. L. Troy, and E. M. Sevick Muraca, “Fluorescence-lifetime determination in tissues or other scattering media from measurement of excitation and emission kinetics,” Appl. Opt. **35**, 2325–2332 (1996). [CrossRef] [PubMed]

**17. **J. S. Reynolds, T. L. Troy, R. H. Mayer, A. B. Thompson, D. J. Waters, K. K. Cornell, P. W. Snyder, and E. M. Sevick-Muraca, “Imaging of spontaneous canine mammary tumors using fluorescent contrast agents,” Photochem. Photobiol. **70**, 87–94 (1999). [CrossRef] [PubMed]

**18. **A. Corlu, R. Choe, T. Durduran, M. A. Rosen, M. Schweiger, S. R. Arridge, M. D. Schnall, and A. G. Yodh, “Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,” Opt. Express **15**, 6696–6716 (2007). [CrossRef] [PubMed]

**19. **A. Liebert, H. Wabnitz, H. Obrig, R. Erdmann, M. Moller, R. Macdonald, H. Rinneberg, A. Villringer, and J. Steinbrink, “Non-invasive detection of fluorescence from exogenous chromophores in the adult human brain,” Neuroimage **31**, 600–608 (2006). [CrossRef] [PubMed]

**20. **A. J. Welch and R. Richards-Kortum, “Monte Carlo simulation of the propagation of fluorescent light” in *Laser induced Interstitial Thermotherapy*, G. Muller and A. Roggan, eds., (1995), p. 174–189.

**21. **A. J. Welch, C. Gardner, R. Richard-Kortum, E. Chan, G. Criswell, J. Pfefer, and S. Warren, “Propagation of fluorescent light,” Lasers. Surg. Med. **21**, 166–178 (1997). [CrossRef] [PubMed]

**22. **R. J. Crilly, W. F. Cheong, B. Wilson, and J. R. Spears, “Forward-adjoint fluorescence model: Monte Carlo integration and experimental validation,” Appl. Opt. **36**, 6513–6519 (1997). [CrossRef]

**23. **J. Swartling, A. Pifferi, A. M. K. Enejder, and S. Andersson-Engels, “Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues,” J. Opt. Soc. Am. A **20**, 714–727 (2003). [CrossRef]

**24. **M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. van der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol. **38**, 1859–1876 (1993). [CrossRef] [PubMed]

**25. **J. Steinbrink, H. Wabnitz, H. Obrig, A. Villringer, and H. Rinneberg, “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. **46**, 879–896 (2001). [CrossRef] [PubMed]

**26. **J. R. Mourant, J. P. Freyer, A. H. Hielscher, A. A. Eick, D. Shen, and T. M. Johnson, “Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics,” Appl. Opt. **37**, 3586–3593 (1998). [CrossRef]

**27. **A. Pifferi, R. Berg, P. Taroni, and S. Andersson-Engels, “Fitting of time-resolved reflectance curves with a Monte Carlo model,” in *Advances in Optical Imaging and Photon Migration*, R. R. Alfano and J. G. Fujimoto, eds. (1996), pp. 311–314.

**28. **A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Moller, R. Macdonald, A. Villringer, and H. Rinneberg, “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. **43**, 3037–3047 (2004). [CrossRef] [PubMed]

**29. **M. Patterson and B. Pogue, “Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues,” Appl. Opt. **33**, 1963–1974 (1994). [CrossRef] [PubMed]

**30. **A. Liebert, H. Wabnitz, J. Steinbrink, M. Moller, R. Macdonald, H. Rinneberg, A. Villringer, and H. Obrig, “Bed-side assessment of cerebral perfusion in stroke patients based on optical monitoring of a dye bolus by time-resolved diffuse reflectance,” Neuroimage **24**, 426–435 (2005). [CrossRef] [PubMed]

**31. **H. Wabnitz, M. Moeller, A. Liebert, A. Walter, R. Erdmann, O. Raitza, C. Drenckhahn, J. P. Dreier, H. Obrig, J. Steinbrink, and R. Macdonald, “A time-domain NIR brain imager applied in functional stimulation experiments,” in *Photon Migration and Diffuse-Light Imaging II*, K. Licha and R. Cubeddu, eds., (2005), p. 58590H.

**32. **M. Kohl-Bareis, H. Obrig, K. Steinbrink, K. Malak, K. Uludag, and A. Villringer, “Noninvasive monitoring of cerebral blood flow by a dye bolus method: Separation of brain from skin and skull signals,” J. Biomed. Opt. **7**, 464–470 (2002). [CrossRef] [PubMed]

**33. **T. S. Leung, I. Tachtsidis, M. Tisdall, M. Smith, D. T. Delpy, and C. E. Elwell, “Theoretical investigation of measuring cerebral blood flow in the adult human head using bolus Indocyanine Green injection and near-infrared spectroscopy,” Appl. Opt. **46**, 1604–1614 (2007). [CrossRef] [PubMed]

**34. **Q. Liu and N. Ramanujam, “Scaling method for fast Monte Carlo simulation of diffuse reflectance spectra from multilayered turbid media,” J. Opt. Soc. Am. A **24**, 1011–1025 (2007). [CrossRef]

**35. **J. Steinbrink, A. Liebert, H. Wabnitz, R. Macdonald, H. Obrig, A. Wunder, R. Bourayou, T. Betz, J. Klohs, U. Lindauer, U. Dirnagl, and A. Villringer, “Towards non-invasive molecular fluorescence imaging of the human brain,” Neurodegenerative **5**, 296–303 (2008). [CrossRef]