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 λ_m) the weights of which decrease due to absorption at λ_m.

For calculating the trajectories of the photon bundles we started from the following assumptions:

 

Fig. 1. Schematic of propagation of a photon bundle (blue) between source and detector and generation of fluorescence (red) within a layered tissue model (DTOF - distribution of times of flight of diffusely reflected photons, DTA - distribution of times of arrival of fluorescence photons at the detector).

Download Full Size | PDF

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^(q)_ij traveled by the photon bundle 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^(e)_ij are evaluated. Based on l^(q)_ij and l^(e)_ij, multiple sets of absorption parameters can be applied to the same trajectory as will be discussed below.

The intrinsic absorption coefficients of the j-th tissue layer are denoted by µ_ax,j and µ_am,j 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 µ_afx,j and µ_afm,j at λ_x and λ_m, respectively. The coefficient µ_afm,j 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 µ_am,j and µ_afm,j at a fixed λ_m 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.

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 λ_m by the factor

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_i. 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 c=c₀/n (c₀=3·10⁸ m/s, n - refractive index of the tissue) in all layers and at λ_x as well as λ_m, the total travel time from source to detector of both the excitation and the fluorescence photon bundle 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/τ_F). Hence, when considering statistical moments of the DTA, τ_F simply adds to the first moment of the DTA and τ ² _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 λ_m, such change in the location of conversion leads to a slight modification of the weight of the fluorescence photon bundle, i.e.

instead of Eq. (6), where

and l ^(q)* _ij is the distance traveled in jth layer by ith 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_m,i according to Eq. (6).

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_q and z_q are coordinates of the qth 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

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.

 

Fig. 2. Comparison of distributions of times of arrival of fluorescence photons calculated for a homogeneous semi-infinite medium by the analytical solution of the diffusion equation (solid lines) with the results of Monte Carlo calculations (symbols) for several source-detector separations.

Download Full Size | PDF

Table 1. Initial homogeneous optical properties of tissue and fluorescing dye distributed in tissue

View Table

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_F). The comparison of integrals was performed after normalization to the sum of all values in the set of simulations.

 

Fig. 3. Validation of the Monte Carlo code for a homogeneous semi-infinite medium. Moments of DTAs (N_Ftot, <t>_F, V_F) simulated at 6 source-detector separations vs. absorption coefficient of the tissue at the excitation wavelength µ_ax (upper row) and at emission wavelength µ_am (lower row), other optical properties as given in Table 1. The solid lines represent the data obtained from the diffusion model [Eq. (13)], the symbols mark the results of Monte Carlo simulations.

Download Full Size | PDF

 

Fig. 4. Moments of DTAs vs. absorption coefficient of the dye at the excitation wavelength µ_afx (upper row) and at emission wavelength µ_afm (lower row). The lines represent the data obtained from the diffusion model [Eq. (13)], the symbols mark the results of Monte Carlo simulations.

Download Full Size | PDF

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.

 

Fig. 5. Spatial distribution of generation probabilities for homogeneous optical properties (a), for increased µ_afx in the lower layer (b) and increased µ_afx in the upper layer (c). Dimensions in mm. The gray scale is linear. In each panel black corresponds to the individual maximum.

Download Full Size | PDF

 

Fig. 6. Changes in moments (integral, mean time of flight and variance) of the DTAs (red squares, scale on the left axis) and DTOFs (blue circles, scale on the right axis) calculated by Monte Carlo simulations for increasing µ_afx in lower (filled symbols) and upper (open symbols) layers of the tissue model; all other optical properties as given in Table 1.

Download Full Size | PDF

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>_F on the spatial distribution of the fluorophore. In particular, in case of high dye concentration in the lower tissue compartment, <t>_F is expected to be larger than if the dye is distributed homogeneously in the whole medium

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>_F and V_F reveal specific changes caused by increasing µ_afx 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, V_F and <t>_F decrease again due to saturation effects.

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 µ_afx 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(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 µ_afx 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 <t>_F as well as V_F can be observed (see Fig. (7), third row).

 

Fig. 7. Changes of moments (normalized integral, mean time of flight and variance) of the DTAs (red lines) and DTOFs (blue lines) obtained from Monte Carlo simulations for interoptode distances of 2 cm (thick lines) and 3 cm (thin lines). The assumed time courses of µ_afx in lower (cyan symbols) and upper (magenta symbols) layers of the model are presented in the left column: homogeneous bolus (first row), slowly decaying bolus in the upper layer only (second row), fast decaying bolus in the lower layer only (third row).

Download Full Size | PDF

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 µ_afx are difficult to estimate for the in-vivo situation. Therefore, we analyzed how the ratio of amplitudes of changes of µ_afx 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 Δ<t>_F as well as ΔV_F 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.

 

Fig. 8. Changes of moments (normalized integral, mean time of flight and variance) of the DTAs (red lines) and DTOFs (blue lines) obtained from Monte Carlo simulations for interoptode distances of 2 cm (thick lines) and 3 cm (thin lines). The assumed time courses of µ_afx (left column) in lower (cyan symbols) and upper (magenta symbols) layers reflect realistic changes of concentration of the fluorophore in the intra- and extracerebral tissue compartments.

Download Full Size | PDF

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 µ_afm, 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.

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>_R. These result matches well with previous 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.