## 1.

## Introduction and Context

Modeling of absorption, diffusion, and fluorescence pheno-mena in a biological tissue allows the study of the distribution of light that propagates in it, while taking into account its intrinsic optical properties.^{1, 2} These physico-optical characteristics depend chiefly on the concentrations, sizes, morphologies, spatial orientations, biochemical composition, and metabolic activity of the tissue compounds as well as their structural organization (cell and conjonctive tissue layers varying according to the type of organ, for instance).^{3} These elements are fundamental characteristics of the gold standard histological diagnosis of healthy and pathological states of biological tissues. Among the main light-absorbing tissues (chromophores), melanin, hemoglobin, and water are to be found.^{4} The elastic scattering of light in tissues is due to variations of the indices of refraction between intra- and extracellular elements (nuclear cells, cell organelles, etc).^{5} Finally, multiple fluorescing macromolecules also exist naturally in living tissues such as proteins (collagen, elastin, keratin), enzymes and coenzymes (various forms of flavins, NADH), or porphyrins.^{1, 6} These intrinsic fluorophores have relatively different but overlapping excitation and emission spectra cove-ring spectral bandwidths of
$280\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}450\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
and
$350\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}700\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, respectively. These fundamental elements are targeted by optical methods for *in-vivo* diagnosis and biological tissue characterization.^{7, 8, 9}

For a given tissue (especially epithelial ones), the presence and concentration of many fluorophores vary with the depth, and according to pathological changes, including cancer. Several studies show that a classification of the states of healthy, hyperplastic, and dysplastic tissues is made possible by using different wavelength bands, the shape of spectra, and the relative amplitudes of emission intensities (including the contrast between porphyrins and NADH/FAD, i.e.,
${\mathrm{I}}_{610\u2013690}$
versus
${\mathrm{I}}_{440\u2013535}$
),^{10, 11} or by in-depth access to the intrinsic fluorescence spectra (the fluorescence spectra acquired at the surface are modulated by the absorption from chromophores or fluorophores and the diffusion of the medium).^{12, 13, 14} For example, Biswal
^{15} show that the profile of autofluorescence spectra obtained with
$440\text{-}\mathrm{nm}$
excitation is clearly altered by the characteristics of porphyrin absorption dips at 540 and
$580\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. These studies also highlight that the measured bulk tissue autofluorescence results from the superposition of emissions of these molecules, whose absorption (excitation) and emission spectra overlap. For example, the absorption spectrum of elastin
$\left(320\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}420\phantom{\rule{0.3em}{0ex}}\mathrm{nm}\right)$
overlaps with the second peak of the NADH absorption spectrum and a portion of the one of collagen, as well as those of porphyrins and the first peak FAD
$\left(360\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}420\phantom{\rule{0.3em}{0ex}}\mathrm{nm}\right)$
, so that the emission spectra of fluorescent proteins (elastin, collagen, keratin,
$360\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}460\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
) and NADH
$\left(400\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}530\phantom{\rule{0.3em}{0ex}}\mathrm{nm}\right)$
respectively overlap with the absorption spectrum of porphyrins and the second absorption peak of FAD.^{16} The consideration and analysis of these interactions come together to properly model the interactions of absorption/multiple excitation and spectral emission of several fluorophores.

Consequently, knowledge and/or identification of tissue optical properties is required in many clinical applications such as *in-vivo* diagnosis and noninvasive follow up, photodynamic therapy (PDT), or laser surgery.^{17} A large part of photodiagnostic methods target cancerology applications and especially studies of epithelial tissues (hollow organs, skin) where 85% of the cancers occur. The implementation of spectroscopic methods *in situ* can easily be performed by means of multiple-fiber probes (i.e., called “point spectroscopy” or “optical biopsy”).^{6, 18} The diagnostic efficiency of such approaches has already been validated separately for autofluorescence, diffuse reflectance, and elastic scattering spectroscopies on numerous tissues (skin,^{18, 19} cervix,^{20} oral mucosa,^{21, 22} bladder,^{23} and breast^{24}). However, it can still be improved by combining those methods, as shown for the cervix,^{8, 25} breasts,^{26} head and neck,^{14} bronchial tree,^{27} and oral cavity.^{28}

Models of light propagation in tissues are of importance and needed in a number of works where the extraction of local intratissular information based on global surface measurements is a challenge (for example, fluence spatial distribution, oxygen or chromophore concentration, intrinsic fluorescence quantification, and optical coefficient identification).^{2, 3, 9, 29, 30, 31, 32}

The Monte Carlo method has been widely used in the past
$20\phantom{\rule{0.3em}{0ex}}\text{years}$
in many studies for simulating the propagation of light in turbid and complex media, such as biological tissues.^{33} It consists in propagating a packet of elementary energy (photon or group of photons) step by step, with random sampling of the length and direction of the displacement steps, as well as the proportion of energy absorbed in each elementary voxel of the medium, with reference to probability distribution functions depending on the optical coefficients of the medium.

The advantage of this statistical simulation method lies in first providing accurate results, whatever the interfiber or source-detector distances, and second in providing the detailed spatial distribution of the fluence while taking into account a large number of parameters: the intensity spectrum of the excitation source, the locations and geometries of both illumination and light collection, the local coefficients of absorption and diffusion, the refractive indices and dioptres, etc.^{34}

Wilson and Adam^{35} and Flock
^{36} were among the first to introduce the use of the Monte Carlo method for modeling light-tissue interactions and propagation in multilayer structures. They validated their algorithm on highly diffusing phantoms at a unique wavelength of
$633\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. In 1993, Wu, Feld, and Rava^{37} proposed a monolayer theoretical model coupling absorption, diffusion, and fluorescence, and taking the quantum yield as well as an isotropic emission of fluorescence into account. An experimental validation performed on arterial tissues highlighted that intrinsic fluorescence information may be extracted from global colocalized measurements of fluorescence and diffuse reflectance.

Between 1992 and 1995, Wang, Jacques, and Zheng^{34, 38} developed Monte-Carlo for multilayered media (MCML), a reference program in C-code for simulation in multilayer models, freely available on the web, and widely used since then for simulating light absorption and light diffusion between 350 and
$1000\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, in various experimental configurations^{39, 40} and on different types of tissues.^{3, 5, 33, 41}

Monte Carlo methods are widely used to model the fluorescence emitted by the tissues to study the effects of excitation/emission geometries on fluorescence, as well as those of the absorption and diffusion properties of the medium, and to study the links between bulk measured fluo-rescence and intrinsic fluorescence at different depths of tissue.^{42} They have been applied by Keijzer
^{43} and Welch
^{44} to study the human aorta autofluorescence, depending in particular on excitation and emission geometries and on the fluorescence quantum yield; by Zonios
^{45} to connect the spectra measured in clinics and the human colon tissue microstructure; by Qu
^{12} to measure the optical properties to study the distortion spectra of autofluorescence depending on the optical properties of lung tissues; and by Drezek
^{13} to determine the relative contribution of two fluorophores present in the epithelial (NADH) and stromal (collagen) layers of a bilayer tissue model. In 1997, Zeng
^{2} simulated fluorescence excitation and emission in two steps: Monte Carlo calculation of the fluence distribution at one single excitation wavelength of
$442\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, then fluorophore excitation with reference to that distribution and final calculation of the propagation for fluorescence emission (seven-layer model,
$470\text{-}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}750\text{-}\mathrm{nm}$
bandwidth). Following the same approach, Welch in 1997^{44} applied Monte Carlo to a multilayer model for simulating the propagation of an excitation light at
$476\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
and the emission of fluorescence light at
$600\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
to study the factors that influence the shape of the autofluo-rescence spectra measured. They showed that 90% of the emitted fluorescence measured comes from thicknesses of tissue up to
$0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
(at
$476\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
). This algorithm was then modified by Drezek in 2001, and applied to cervical tissues (monoexcitation at
$380\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
).^{13} In 2001, Müller
^{29} validated a correction method of the absorption and diffusion contributions to recover an intrinsic fluorescence spectrum, taking into account the spectral dimension of the optical coefficients. However, the algorithm proposed was not validated for fluorophores with high quantum efficiency and high absorption, and secondary emission of fluorescence was not implemented. In 2003, Liu, Zhu, and Ramanujam^{33} were the first to describe precisely an algorithm of statistical simulation of fluorescence with one excitation wavelength and one emission wavelength. They also gave the expressions of the optical coefficients of the model with reference to absorption and diffusion by fluorescence. A rigorous experimental validation was carried out with fluorescence and diffuse reflectance on phantoms for an excitation wavelength of
$460\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
and an emission one at
$520\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. Ma
^{32} compared the results obtained with a simplified analytical model of fluorescence and a simplified method of Monte Carlo for simulating time-resolved fluorescence. They observed that the Monte Carlo method gives better results when fluorophores lie in areas close to the source and to the detection fibers, as well as for short interfiber distances. A Monte Carlo code was developed by Vishwanath and Mycek^{46} in 2005 to simulate excitation and time-resolved fluorescent light propagation at multiple wavelengths in a bilayered turbid model. Experimental validation was conducted on two-phase (liquid and solid) layered tissue-simulating phantoms. In 2006, Chandra
^{47} obtained good agreement in comparing reflectance and fluorescence spectroscopic measurements made on cartilage tissue samples with simulation results using their time-resolved photon propagation code for multifluorophore emission.

The Monte Carlo technique developed by Churmakov
^{48, 49} simulates the spatial distribution of fluorescence excitation by combining a scattering-center matrix (calculating scattering paths) and an absorption calculation based on the path length applied to a multilayer model of human skin. In 2008, Palmer and Ramanujam^{50} experimentally validated a new extraction method of intrinsic fluorescence with a Monte Carlo model based on the work of Swartling
^{51} and Graff
^{52} The model is validated using liquid synthetic phantoms with a single fluorophore in light-induced fluorescence spectroscopy (
${\lambda}_{\mathit{exc}}=330\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
,
${\lambda}_{\mathit{em}}=350\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}500\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
) and diffuse reflectance spectroscopy
$\left(330\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}600\phantom{\rule{0.3em}{0ex}}\mathrm{nm}\right)$
. The proposed model optimizes the fast-track approach based on the reciprocity principle to determine the absorption energy density grid, applying scaling relations to the two grids according to the set of optical parameters and using quasidiscrete Hankel transform for the convolution of the two grids (faster). Liebert
^{53} also developed a new Monte Carlo algorithm that simulates photon fluorescence generation and diffusion in a multilayer time-resolved model. The principle of absorption and scatte-ring of a packet of photons computes the probability of absorption fluorescence as a conversion (at an excitation wavelength) of a specific percentage of this packet; the emitted photons are then propagated at the emission wavelength. This methodology implies an isotrope distribution, values of reduced scattering coefficients identical at both emission and excitation wavelengths, and does not take the reabsorption into account.

Thus, one of the current difficulties resides in the coupled simulation of absorption, diffusion, and multiple fluorescence in a multilayer medium, over a broad spectral band. Furthermore, these works also confirm the necessity for a rigorous experimental validation of the simulation algorithms prior to their application to biological tissues.

Several of the algorithms that have been developed mostly consider multilayer models with a maximum of one fluorophore per layer and without taking into account any effect of reabsorption or re-emission. The work presented in this work is an incremental contribution to the continuation of the works of Liu, Zhu, and Ramanujam^{33} and Churmakov
^{48, 49} in particular by providing a solution taking into account both the spectral dimension (full spectra of absorption and emission of the fluorophores) and reabsorption/re-emission by several fluorophores whose absorption and emission spectra overlap.

Consequently, this work presents the principle of an algorithm developed for simulating the spectral absorption and emission of several fluorophores in a multilayer model. The experimental validation of the program implemented in C is performed with reflectance data measured at different distances between excitation and collection fibers on various mono- and multilayer phantoms containing absorbing, diffusing, and fluorescing compounds.

## 2.

## Materials and Methods

## 2.1.

### Monte Carlo Simulations

## 2.1.1.

#### Basic principle and fundamental algorithm

The simulation of spatially resolved light transport in multilayered models is already well documented in the literature.^{34, 54, 55} We used these bases to develop our algorithm for multiple fluorescence.

Our modeling configuration of a semi-infinite medium with $N$ layers parallel to each other is considered, with one excitation optical fiber and several emission or collecting fibers (Fig. 1 ). Each layer is described by the following para-meters: the thickness $e$ (mm), the absorption coefficient ${\mu}_{a}\left(\lambda \right)$ $\left({\mathrm{cm}}^{-1}\right)$ , the scattering coefficient ${\mu}_{s}\left(\lambda \right)$ $\left({\mathrm{cm}}^{-1}\right)$ , the refractive index $n$ , and the anisotropy factor $g\left(\lambda \right)$ .

In short, the principle of random migration of photons consists in displacing a large number of photons or photon packets (usually a few hundreds of thousands at least) in directions and over distances drawn at random with reference to density probability functions characterizing the medium $({\mu}_{a},{\mu}_{s},g)$ .

First, the source photons with an initial “weight” equal to unity are injected at origin coordinates from a random position inside the diameter of the excitation fiber. Each new 3-D direction of propagation is determined by statistical sampling of:

• the deflection angle $\theta $ in $[0,\pi ]$ , i.e., the angle between the previous and new direction lines in the same plane, is defined by its cosine based on the Henyey-Greenstein scattering phase function

^{56}as:## Eq. 1

$\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta =\frac{1}{2g}[1+{g}^{2}-{\left(\frac{1-{g}^{2}}{1-g+2g\xi}\right)}^{2}],$with $\xi $ being a random variable uniformly distributed between 0 and 1.

• the azimuthal angle (angle with the previous plane) $\xi $ in $[0,2\pi ]$ is defined as:

## 2.1.2.

#### Algorithms previously developed for fluorescence simulation

Liu, Zhu, and Ramanujam^{33} experimentally verified Monte Carlo modeling of fluorescence and diffuse reflectance measurements in turbid, tissue-like phantom models. By consi-dering the optical parameters of a fluorophore
${\mu}_{af}$
(absorption coefficient of the fluorophore) and
$\varphi $
(fluorescence quantum yield), the authors defined the absorption probability for fluorescence as follows:

Churmakov
^{48} simulated the spatial fluorescence distribution within the human skin. In their simulation algorithm, each photon packet produces only one fluorescent event. Their simulation results suggest that the distribution of autofluo-rescence is significantly suppressed in the near-infrared spectral region. In 2004, Chang
^{57} compared the results of Monte Carlo simulations with an analytical model of fluorescence in a two-layer model applied to normal and preneoplastic epithelial tissue. Their results show that the analytical model provides a good description of fluorescence for the optical properties of normal and precancerous cervical epithelium and stroma. Finally, Matuszak, Sawow, and Wasilewska-Radwanska^{58} implemented the simulation of fluorescence spectra in two steps: the simulation of the excitation spectrum to obtain the intensity of light inside cells containing a fluorophore, and then the simulation of emitted photons (emission spectra).

## 2.1.3.

#### Simulation algorithm developed for multiple fluorescence (proposed solution)

Each layer of the medium (Fig. 1) is considered homogeneous with a thickness $e$ , a refractive index $n$ , an anisotropy factor $g$ , and absorption and scattering coefficients, respectively ${\mu}_{a}$ and ${\mu}_{s}$ . In our model, each layer can also contain one or more fluorophores expressed by the following complementary parameters: the absorption coefficient ${\mu}_{af}\left(\lambda \right)$ , the fluorescence quantum yield $\varphi \left(\lambda \right)$ , the emission spectra ${S}_{\mathit{em}}^{{f}_{i}}\left(\lambda \right)$ , and the molecular concentration of the fluorophores.

The implementation of multiple fluorescence involves three principal steps: first the absorption of a photon by a fluorophore as a function of its absorption spectra and concentration in the medium, second the choice of an emission wavelength with reference to its emission spectra, third the generation of the fluorescence photon in a new direction. These various stages are detailed in the flowchart of Fig. 2, which is part of the main classical Monte Carlo algorithm.

Considering an absorbing and diffusing medium containing ${n}_{f}$ fluorophores with absorption coefficients ${\mu}_{a{f}_{i}}\left(\lambda \right)$ and associated quantum yields ${\varphi}_{i}\left(\lambda \right)$ $(i\u220a[1,{n}_{f}])$ , the stepsize $\Delta s$ is calculated in the following way:

## 5.

At each absorption event, a first test ${T}_{1}$ is carried out to decide whether it is statistically due to a fluorescent compound (Fig. 2). The probabilities ${p}_{\mathit{abs}}$ that a photon is absorbed by a “pure” absorber and ${p}_{{f}_{i}}$ that a photon is absorbed by a fluorophore ${f}_{i}$ can respectively be defined as:

## Eq. 6

${p}_{\mathit{abs}}=\frac{{\mu}_{a}}{{\mu}_{t}}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{p}_{{f}_{i}}=\frac{{\mu}_{a{f}_{i}}}{{\mu}_{t}}.$## Eq. 7

${\xi}_{\mathrm{norm}}=\xi \times \left(\frac{{\mu}_{a}+{\sum}_{i}{\mu}_{a{f}_{i}}}{{\mu}_{t}}\right),$## Eq. 8

${p}_{\mathit{abs}}+\sum _{i=0}^{{n}_{f}-1}{p}_{{f}_{i}}\u2a7d{\xi}_{\mathrm{norm}}\u2a7d{p}_{\mathit{abs}}+\sum _{j=i+1}^{{n}_{f}}{p}_{{f}_{i}}\phantom{\rule{1em}{0ex}}({p}_{{f}_{0}}=0,\phantom{\rule{0.3em}{0ex}}{n}_{f}\u2a7e2).$In the first case (photon absorbed by a “pure” absorber), the photon’s weight becomes:

## Eq. 9

${\omega}_{a}={\omega}_{b}(1-\frac{{\mu}_{a}}{{\mu}_{a}+{\sum}_{i}{\mu}_{a{f}_{i}}+{\mu}_{s}}),$In the second case (photon absorbed by a fluorophore ${f}_{i}$ ), a second test ${T}_{2}$ is performed to decide whether conversion into fluorescence happens, depending on the fluorophore quantum yield (Fig. 2). Again, the principle consists in generating a random number $\xi $ uniformly distributed between 0 and 1.

• If this random number is higher than the quantum yield value $\varphi $ $(\varphi <\xi \u2a7d1)$ , then the incident photon is not converted into fluorescence. It continues its way until the next event, after having updated its energy level with reference to ${\mu}_{a{f}_{i}}$ as follows:

• If $0\u2a7d\xi \u2a7d\varphi $ , then the incident photon is fully absorbed and a new photon of fluorescence is generated at the same coordinates in the medium. The first propagation direction of the new emission photon is randomly sampled between 0 and $360\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ .

The emission wavelength of the new photon is also randomly sampled according to a probability distribution corres-ponding to the emission spectrum of the fluorophore, whatever the absorption wavelength. For example, if a fluorophore has an emission peak at a given wavelength, then the highest the probability emission of a fluorescence photon would be at this wavelength.

The values of the emission spectrum ${S}_{\mathit{em}}^{{f}_{i}}\left(\lambda \right)$ for each fluorophore ${f}_{i}$ being known, the following probability of emission of a fluorescence photon at a wavelength ${\lambda}_{j}$ is associated:

## Eq. 11

$p\left[\text{photon}\left({\lambda}_{j}\right)\right]={v}^{\prime}\left({\lambda}_{j}\right)=\frac{{S}_{\mathit{em}}\left({\lambda}_{j}\right)}{{\sum}_{j=1}^{{n}_{p}}{S}_{\mathit{em}}\left({\lambda}_{j}\right)},$A photon of fluorescence can cause another event of absorption fluorescence (the “cascade” effect). In addition, an incident photon absorbed at a wavelength ${\lambda}_{\mathit{abs}}$ cannot be converted into an emission photon with a wavelength lower than ${\lambda}_{\mathit{abs}}$ .

In the program based on the flowchart in Fig. 2, a specific file was created for each fluorophore in each layer, in which the values of the absorption and emission spectra, quantum yield, and absorption coefficient of fluorophores are stored. The simulation is made in the ascending way, i.e., from the lowest wavelength to the highest one, with $10\text{-}\mathrm{nm}$ sampling steps.

Finally, the function of cumulative distribution ${C}_{{v}^{\prime}}$ such as

## Eq. 12

$\{\begin{array}{ll}\text{if}\phantom{\rule{0.3em}{0ex}}\xi \u2a7d{C}_{{v}^{\prime}}\left({\lambda}_{1}\right)& \Rightarrow \text{emission}\phantom{\rule{0.3em}{0ex}}\text{at}\phantom{\rule{0.3em}{0ex}}\text{wavelength}\phantom{\rule{0.3em}{0ex}}{\lambda}_{1}\\ \text{if}\phantom{\rule{0.3em}{0ex}}{C}_{{v}^{\prime}}\left({\lambda}_{j-1}\right)<\xi \u2a7d{C}_{{v}^{\prime}}\left({\lambda}_{j}\right)& \text{for}\phantom{\rule{0.3em}{0ex}}j\u2a7e2,\\ & \Rightarrow \text{emission}\phantom{\rule{0.3em}{0ex}}\text{at}\phantom{\rule{0.3em}{0ex}}\text{wavelength}\phantom{\rule{0.3em}{0ex}}{\lambda}_{j}\end{array}\phantom{\}},$## 2.2.

### Instrumentation Setup for Experimental Validation

We used a spatially resolved reflectance spectroscopy configuration. It consists in a laser diode (Laser 2000, France) at ${\lambda}_{\mathit{exc}}=(410\pm 5)\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ for the fluorescence excitation of all phantoms except phantom ${L}_{3}$ (see Sec. 2.3) for which a deuterium-tungsten halogen light source (DH2000, Ocean Optics, Dunedin, Florida) bandpass filtered at $(650\pm 20)\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ was used. Backscattered fluorescence emission spectra were measured using an imaging spectrograph (iHR320, HoribaJobin Yvon, France) with a back-illuminated charge-coupled device (CCD) detector, and a diffraction grating ( $150\phantom{\rule{0.3em}{0ex}}\mathrm{gr}\u2215\mathrm{mm}$ , blaze wavelength $500\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ). A multiple fiber optic probe (37 $\mathrm{Si}\u2215\mathrm{Si}$ fibers with $200\text{-}\mu \mathrm{m}$ core diameter and 0.22 numerical aperture) coupled to a specific in-line fiber bundle at the entrance of the spectrograph allows multitrack operations, i.e., the simultaneous acquisition of up to 13 spectra correspon-ding to 13 different distances between illumination and collection fibers $\left(271\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1542\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\right)$ . Results shown in the present work are given for only one distance between illumination and acquisition fibers $\left(271\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\right)$ , but the results obtained were verified and validated for all other distances.

To improve signal-to-noise ratio, three spectra were automatically acquired for averaging. All experimental spectra were preprocessed to be free of any spectral distortions induced by instrumentation (spectral correction in wavelength and intensity). To improve reproducibility, a tripod was used to maintain the tip of the fiber optic probe in gentle contact with the surface of the phantom and perpendicular to it.

## 2.3.

### Test Phantoms

Phantoms were prepared with specific fluorophores, absor-ption, and scattering properties and dimensions, all of which are detailed hereafter. Liquid phantoms were made of absorber, scatterer, and one or two fluorophores, ethanol being used as thinner for the liquid phantoms since emission and absorption spectra of the fluorophores used here are given with ethanol as a thinner.^{38} The relevance of the monofluorophore phantoms was to check the good correlation between measurements and simulation results. The relevance of multifluorophore phantoms was to check if the cascade effect was correctly taken into account. Then, multilayer solid phantoms were prepared in two and three layers with one or two fluorophores.

India ink (Super Black India ink, Lefranc/Bourgeois, France) was selected as the absorber. The absorption coefficient of pure India ink is available in Ref. 33, and the absor-ption coefficient of India ink diluted with water (Kabi-Fresenius) at 0.011% was measured using a spectrophotometer (Lambda EZ210, in DCPR Department of Physics Chemistry of the Reactions, Nancy, France).

Intralipids (Sigma Aldrich, Saint Louis, Missouri) were chosen as scatterer. They are often used as scatterers because they do not have a strong absorption band in the visible. They are made of intravenously administered nutrients consisting of an emulsion of phospholipid micelles and water, easily avai-lable and relatively inexpensive.^{36, 59, 60} The optical coefficients of intralipids are available in Ref. 36. Figure 3 represents the absorbance and the scattering coefficient as a function of wavelength for pure India ink and intralipids, respectively.

Three fluorophores were selected for our experiments with reference to their distinct or partially overlapping emission spectra: eosin Y (
${\mathrm{C}}_{20}{\mathrm{H}}_{6}{\mathrm{Br}}_{4}{\mathrm{Na}}_{2}{\mathrm{O}}_{5}$
, Sigma Aldrich), fluore-scein (
${\mathrm{C}}_{20}{\mathrm{H}}_{12}{\mathrm{O}}_{5}$
, Fluka) and cryptocyanine (
${\mathrm{C}}_{25}{\mathrm{H}}_{25}\mathrm{I}{\mathrm{N}}_{2}$
, Sigma Aldrich). Indeed, eosin Y dissolved in ethanol has a quantum yield of 0.67 at
$490\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, fluorescein a quantum yield of 0.79 at
$425\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, and cryptocyanine a quantum yield of 0.007 at
$650\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
.^{38} Moreover, fluorescein is referred to several articles for the creation and validation of equipment.^{61, 62}

Figure 4 represents the absorption spectrum of fluorescein, eosin, and cryptocyanine in ethanol for liquid phantoms,^{38} and in water for solid phantoms, also measured in DCPR (Nancy).

For all the phantoms, the fluorophores were systematically weighed to obtain a quantity of $10\phantom{\rule{0.3em}{0ex}}\mu \mathrm{mol}$ of product. All simulations were performed with 1,000,000 photons launched at 40 different wavelengths between 400 and $800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ with $10\text{-}\mathrm{nm}$ steps (i.e., 40 million photons per spectrum), to ensure convergence of the results.

## 2.2.1.

#### Liquid phantoms

Four liquid phantoms were prepared: ${\mathrm{L}}_{1}$ ( $3.32\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of fluorescein, $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of ethanol), ${\mathrm{L}}_{2}$ ( $6.47\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of eosin, $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of ethanol), ${\mathrm{L}}_{3}$ ( $4.8\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of cryptocyanine, $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of ethanol), and ${\mathrm{L}}_{4}$ ( $3.32\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of fluorescein, $4.8\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of cryptocyanine in $16\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of ethanol as thinner). All solution phantoms were poured into a black container comparable to one expanded polystyrene container, until reaching a height of $5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . The ${\mathrm{L}}_{4}$ phantom mixing two fluorophores was prepared to test the cascading effect itself (or secondary fluorescence emission due to another fluorescence photon).

## 2.2.2.

#### Solid mono- and multilayer phantoms

We used a solution of agarose (Agarose L, VWR Prolabo, France) at 5% ( $0.35\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of agarose powder and $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of distilled water) that has the property to become liquid at $82\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$ and solid around $35\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$ . Contribution of light scattering and fluorescence from agarose was checked to be negligible compared to the contributions of fluorophores at the wavelengths of excitation chosen for our experiments. Mono- and multilayer thicknesses were measured on sample sections with a Vernier micrometer (error: $\pm 0.02\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ).

A first set of bilayered phantoms with one fluorophore ( ${\mathrm{MC}}_{1}$ , ${\mathrm{MC}}_{2}$ , and ${\mathrm{MC}}_{3}$ ) was prepared, having at the bottom (second layer, thickness ${e}_{1}=3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) $3.32\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of fluorescein and $0.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ of agarose powder in $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ distilled water, and at the top (first layer, thickness ${e}_{1}=2.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) either $5\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of India ink diluted at 0.011% $\left({\mathrm{MC}}_{1}\right)$ , or $5\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of intralipids at 10% $\left({\mathrm{MC}}_{2}\right)$ , or a mixture of both $\left({\mathrm{MC}}_{3}\right)$ . For this set of phantoms, the first layer was kept in solution to prevent the optical properties from being modified by heating. A complementary single-layer phantom $\left({\mathrm{MC}}_{4}\right)$ was also made with $3.32\text{-}\mathrm{mg}$ fluorescein, $0.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ of agarose powder, and $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of distilled water (thickness $3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) to compare the spectral res-ponses obtained for the previous bilayered phantoms with the ones of the second layer alone.

A second set of two- or three-layered phantoms was made with two different fluorophores. One bilayer phantom $\left({\mathrm{MC}}_{5}\right)$ was composed of $6.47\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of eosin, $0.35\text{-}\mathrm{g}$ agarose, and $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of water on the top layer (thickness $3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ), and of $3.32\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of fluorescein, $0.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ of agarose powder, and $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of distilled water at the bottom layer (thickness $2.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ). The trilayered phantom $\left({\mathrm{MC}}_{6}\right)$ was composed of a top layer of (thickness: $2.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) India ink diluted at 0.011%; a middle layer of (thickness $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) $6.47\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of eosin, $0.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ of aga-rose powder, and $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of distilled water; and a bottom layer of (thickness $5.4\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) $3.32\phantom{\rule{0.3em}{0ex}}\mathrm{mg}$ of fluorescein, $0.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}$ of agarose powder, and $8\phantom{\rule{0.3em}{0ex}}\mathrm{mL}$ of distilled water.

Table 1 recapitulates the characteristics and composition of the six types of multilayer phantoms tested.

## Table 1

Characteristics of multilayer phantoms.

Phantom | MC1 | MC2 | MC3 | MC4 | MC5 | MC6 |
---|---|---|---|---|---|---|

First layer | $5\text{-}\mathrm{mL}$ India ink diluted at 0.011% | $5\text{-}\mathrm{mL}$ intralipids at 10% | $2.5\text{-}\mathrm{mL}$ solution 1 | $6.47\text{-}\mathrm{mg}$ eosin | ||

$2.5\text{-}\mathrm{mL}$ solution 2 | $0.35\text{-}\mathrm{g}$ agarose powder | solution 1 | ||||

— | $8\text{-}\mathrm{mL}$ distilled water | |||||

(Solution 1) | (Solution 2) | (Solution 4) | ||||

${e}_{1}=2.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ${e}_{1}=3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ${e}_{1}=2.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ||||

Second layer | $3.32\text{-}\mathrm{mg}$ flurescein | |||||

$0.35\text{-}\mathrm{g}$ agarose powder | ||||||

$8\text{-}\mathrm{mL}$ distilled water | — | Solution 3 | Solution 4 | |||

(Solution 3) | ||||||

${e}_{2}=3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ${e}_{2}=2.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ${e}_{2}=2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ||||

Third layer | — | — | — | Solution 3 | — | Solution 3 |

${e}_{3}=3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ | ${e}_{3}=5.4\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

To quantify the difference between experimental and corresponding simulated spectra, we calculated, for all phantoms, the mean of normalized differences as follows:

## Eq. 13

${\text{error}}_{\text{mean}}=\frac{1}{N}\sum _{i=1}^{N}\mid \frac{{S}_{\mathrm{exp}}\left({\lambda}_{i}\right)-{S}_{\mathrm{sim}}\left({\lambda}_{i}\right)}{{S}_{\mathrm{exp}}\left({\lambda}_{i}\right)}\mid ,$## 3.

## Results and Discussion

## 3.1.

### Primary Validation by Comparison to Simulation Results from Literature

The proper operation of our program, including our fluo-rescence algorithm, was first verified for absorption and diffuse reflectance without fluorophores. Therefore, our results were compared to other Monte Carlo simulation results referenced in the literature with identical model parameter-ization.^{34, 54, 63}

The total diffuse reflectance of a monolayer slab of a turbid medium was computed with the following optical properties: absorption coefficient
${\mu}_{a}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, scattering coefficient
${\mu}_{s}=90\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, anisotropy factor
$g=0.75$
, relative refractive index
$n=1$
(i.e., refractive-index-matched boundary), and thickness
$e=0.02\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
. Ten Monte Carlo simulations with 50,000 launched photons per wavelength were completed. The average and standard errors of the total diffuse reflectance are summarized in Table 2 together with the results from Van de Hulst,^{63} Prahl,^{54} and Wang, Jacques, and Zheng.^{34} All results agree with each other.

## Table 2

Comparison of published simulation results of average and standard errors of total diffuse reflectance Rd obtained for a slab of thickness e=0.02cm with a matched boundary ( μa=10cm−1 , μs=90cm−1 , g=0.75 , and n=1 ), and for a semi-infinite medium with a mismatched boundary ( μa=10cm−1 , μs=90cm−1 , g=0 , and n=1.5 ).

Model | Source | Rd Average | Rd Error |
---|---|---|---|

Slab with a matched boundary $(e=0.02\phantom{\rule{0.3em}{0ex}}\mathrm{cm})$ | Ref. 63 | 0.09739 | — |

Ref. 54 | 0.09711 | 0.00033 | |

Ref. 34 | 0.09734 | 0.00035 | |

Our program | 0.09737 | 0.00019 | |

Semi-infinite medium with a mismatched boundary | Ref. 64 | 0.26000 | — |

Ref. 54 | 0.26079 | 0.00079 | |

Ref. 34 | 0.25907 | 0.00170 | |

Our program | 0.26204 | 0.00138 |

In the same way, we computed the average and standard errors of the total diffuse reflectance for a semi-infinite turbid medium that has a mismatched refractive index with the ambient medium and optical properties:
${\mu}_{a}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
,
${\mu}_{s}=90\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
,
$g=0$
(isotropic scattering), and
$n=1.5$
. Ten Monte Carlo simulations with 5000 photons per wavelength were completed. Here as well, our results are in good correlation with the ones from Giovanelly,^{64} Prahl,^{54} and Wang, Jacques, and Zheng,^{34} as shown in Table 2.

As a primary validation, those results confirm that our program with fluorescence algorithm added operates properly in the absence of fluorophore in the simulation model.

## 3.2.

### Validation with Phantoms in Solution

All experiments were performed under monochromatic excitation with the laser diode source at $410\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , except the phantom ${\mathrm{L}}_{3}$ performed with a a deuterium-tungsten halogen light source bandpass filtered at $650\pm 20\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ .

The correlation between simulated and measured emission spectra for each fluorophore was first verified using mono-fluorescence phantoms in solution (phantoms ${\mathrm{L}}_{1}$ , ${\mathrm{L}}_{2}$ , and ${\mathrm{L}}_{3}$ ). Experimental and simulation curves are shown in Fig. 5 . The spectra were normalized with reference to the maximum peak value of the spectra. Experimental and simulated spectra were normalized in the same way. A relatively good correlation was obtained, with error values of 2.2, 1.9, and 2.5% respectively for phantoms ${\mathrm{L}}_{1}$ , ${\mathrm{L}}_{2}$ , and ${\mathrm{L}}_{3}$ .

To test the proper operation of secondary fluorescence absorption (i.e., cascading fluorescence emission due to the absorption of another “already emitted” fluorescence photon), results in single excitation (phantoms ${\mathrm{L}}_{1}$ , ${\mathrm{L}}_{3}$ ) and multiple (or cascading) excitation (phantom ${\mathrm{L}}_{4}$ : fluorescein, cryptocyanine, ethanol) were compared (Fig. 6 ).

We observe that there is good correlation between experimental and simulated fluorescence spectra given in Fig. 6 with a mean error of 5.3%. For the phantom ${\mathrm{L}}_{4}$ and for the shortest distance between illumination and acquisition Fibers $\left(271\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\right)$ , standard deviations for the experimental and simulated spectra (respectively ${\sigma}_{\mathrm{exp}}$ , ${\sigma}_{\mathrm{sim}}$ ) are:

The emission intensity spectra of fluorescein between 480 and $650\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ (peak at $540\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ) and of cryptocyanine at wavelengths above $700\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ (emission peak at $730\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ) can be clearly noticed, and their proportional contributions are well reproduced by the simulation algorithm including cascading fluorescence. With reference to the absorption spectrum of cryptocyanine, we can even check in Fig. 6 that the amplitude of the fluorescein emission spectrum is clearly decreased in the $580\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}700\text{-}\mathrm{nm}$ wavelength band compared to the one of phantom ${\mathrm{L}}_{1}$ in Fig. 5, indicating that absorption by cryptocyanine occurs and therefore that the cascade effect does exist.

## 3.3.

### Validation with Multilayered Phantoms

Experimental and simulated fluorescence spectra for multilayered phantoms ${\mathrm{MC}}_{1}$ (India ink/fluorescein), ${\mathrm{MC}}_{2}$ (intralipids/fluorescein), ${\mathrm{MC}}_{3}$ (India ink/intralipids/fluorescein) and ${\mathrm{MC}}_{4}$ (fluorescein) are given in Fig. 7 . The latter also highlights the way in which the fluorescence spectra are modified by the presence of absorbing and diffusing layers. Experimental spectra were normalized with reference to the maximum peak value of the experimental value of phantom ${\mathrm{MC}}_{4}$ , and simulated spectra with reference to the maximum peak value of the simulated value of phantom ${\mathrm{MC}}_{4}$ .

Looking comparatively at the spectra of ${\mathrm{MC}}_{4}$ and ${\mathrm{MC}}_{1}$ phantoms, it can be noticed that the addition of a purely absorbing top layer (India ink) leads to a global amplitude attenuation, while the overall shape of the spectrum remains similar. By comparing ${\mathrm{MC}}_{4}$ and ${\mathrm{MC}}_{2}$ spectra, we may observe that the presence of a diffusing top layer (intralipids) causes strong attenuation of the collected signal and also a change in the shape of the curve. In the case of phantom ${\mathrm{MC}}_{3}$ , the absorbing and diffusing top layer leads to attenuation and shape modification of the reflected spectra.

For these four different configurations, the results obtained experimentally and by simulation agree well. The lowest mean percent error obtained between experimental and simulated spectra is 2.5% for phantom ${\mathrm{MC}}_{4}$ (fluorescein), 6.0% for ${\mathrm{MC}}_{2}$ (intralipids, fluorescein), 7.1% for ${\mathrm{MC}}_{1}$ (India ink, fluorescein), and the maximum mean percent error is 9.4% for ${\mathrm{MC}}_{3}$ (India ink, intralipids, fluorescein). It can be noticed that the mean percent error is slightly higher for the solid phantom ${\mathrm{MC}}_{4}$ (2.5%) than for the phantom ${\mathrm{L}}_{1}$ in solution (2.2%).

Figure 8 represents the measured and simulated fluore-scence spectra for the bilayers phantom ${\mathrm{MC}}_{5}$ with two fluorophores (eosin on the top layer and fluorescein at the bottom layer) and for the trilayer phantom ${\mathrm{MC}}_{6}$ with the same two fluorophore layers and an upper absorbing layer (India ink/eosin/fluorescein). A few explanations can be given as to the shape changes of spectra in Fig. 8. As presented in Fig. 3a, India ink absorbs twice as much at $400\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ as at $800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , therefore altering the spectrum unequally over the wavelength band. Thus, the eosin fluorescence appears proportionally with more intensity than the fluorescein fluorescence. Consequently, the modification of the shape in Fig. 8 is due to the absorption spectra of India ink and to the fact that this was added in a supplementary top layer. The mean percent error calculated between experimental and simulated spectra is 9.3% for the bilayer phantom ${\mathrm{MC}}_{5}$ with two fluorophores (eosin, fluorescein) and 10.6% for the trilayer phantom ${\mathrm{MC}}_{6}$ (India ink, eosin, fluorescein). It should be noted that the mean percent error increases with the addition of an absorbing layer, and with the number of layers (increased thickness). This trend is confirmed for the 12 other distances between illumination and acquisition fibers $\left(340\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1542\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\right)$ . For the $528\text{-}\mu \mathrm{m}$ interfiber distance, the average error varies from 3.1% for ${\mathrm{MC}}_{4}$ to 10.9% for ${\mathrm{MC}}_{6}$ . We can observe that the mean percent error values increase with the interfiber distances and with the medium complexity, i.e., the number of collected photons. Given the number of absorbing, fluore-scing, or diffusing layers, the increases in average error are probably due to the relative increase of statistical noise inherent to Monte Carlo simulation, because we used the same number of photons for all simulations (40 million photons per spectrum) whatever the phantom complexity.

Execution times vary with the composition of the phantoms, and more particularly with the optical parameters used in the simulation. Using a pentium 4, $2.4\text{-}\mathrm{GHz}$ and 512-Mo RAM computer, execution times (20 million photons per spectra) range from $6\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ for ${\mathrm{MC}}_{4}$ (fluorescein, agarose) to $34\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ for three-layered phantom ${\mathrm{MC}}_{6}$ with two fluorophores (India ink, eosin, fluorescein). Simulation time is doubled between the monolayer monofluorophore phantom ${\mathrm{MC}}_{4}$ and bilayer bifluorophore ${\mathrm{MC}}_{5}$ $\left(15\phantom{\rule{0.3em}{0ex}}\mathrm{min}\right)$ .

## 4.

## Conclusion

Anatomical, physiological, and biochemical changes result in modifications of the absorption and light scattering in pathological tissues compared to healthy ones. These local changes of the optical properties affect the characteristics of reflectance diffuse spectra. The measurement of endogenous and exogenous fluorescence from tissue has become a very useful research and clinical diagnostic tool. Many authors have shown that fluorescence can indicate the presence of mali-gnant tissue, but one of the major limitations is that the fluorescence signal is mixed with absorption and scattering, which is intrinsic to the tissue components. As a consequence, Monte Carlo methods are widely used to model the fluore-scence emitted by the tissues. Part of the difficulty lies in the coupled simulation of absorption, diffusion, and multiple fluorescence in a multilayer medium, over a broad spectral band.

In this work, a Monte Carlo algorithm is developed for simulating the spectral absorption, diffusion, and multiple fluorescence in a multilayer medium, over a broad spectral band. As a primary validation, the correct operation of our program, including our fluorescence algorithm, is first verified in relation to other Monte Carlo simulation results referenced in literature for absorption and diffuse reflectance without fluorescence (same parametrization). The main experimental validation consists in comparing the fluorescence spectra measured on phantoms of various complexities to Monte Carlo simulated spectra of their corresponding models (mono- and multilayered). Four phantoms in solution with three different fluorophores and six mono- and multilayered solid phantoms with one or two fluorophores are used. The good correlation between simulated and experimentally measured results is quantified on the basis of mean error calculations. For phantoms in solution, the mean percent errors vary from 1.9% (one fluorophore) to 5.3% (two fluorophores). For multilayered phantoms, these errors range from 6.0% (two layers), to 9.4% for the trilayer phantom. The error values vary according to the complexity of the phantom, i.e., the model (for a fixed number of photons launched). Finally, the proper operation of “secondary fluorescence” (fluorescence emission due to the absorption of another already emitted fluorescent photon) is demonstrated.

This model proves appropriate for applications to biological tissues, in which the contribution of absorbers, scatterers, and fluorophores vary. More specifically, this work will be used for *in-vivo* tissue diagnosis such as precancer and cancer detection, where the identification of proper values of the optical parameters for these tissues takes into account the associated phenomenon of multiple fluorescence.

## Acknowledgments

The authors would like to thank the Region Lorraine and the Ligue Contre le Cancer (CD 52, 54) for their financial support, and Céline Frochot (DCPR, Nancy) for complementary spectrophotometric measurements.

## References

**,” Lasers Med. Sci., 23 267 –276 (2008). https://doi.org/10.1007/s10103-007-0478-2 0268-8921 Google Scholar**

*In vitro*fluorescence measurements and Monte Carlo simulation of laser irradiation propagation in porcine skin tissue**,” J. Photochem. Photobiol., B, 38 (2–3), 238 –240 (1997). https://doi.org/10.1016/S1011-1344(96)00008-5 1011-1344 Google Scholar**

*Reconstruction of**in vivo*skin auto-fluorescence spectrum from microscopic properties by Monte Carlo simulation**,” Physiol. Meas, 23 741 –753 (2002). https://doi.org/10.1088/0967-3334/23/4/312 0967-3334 Google Scholar**

*Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions***,” Lasers Med. Sci., 30 (3), 191 –200 (2002). https://doi.org/10.1002/lsm.10026 0268-8921 Google Scholar**

*Optimal methods for fluorescence and diffuse reflectance measurements of tissue biopsy samples***,” J. Biomed. Opt., 8 (3), 495 –503 (2003). https://doi.org/10.1117/1.1578494 1083-3668 Google Scholar**

*In vivo*endoscopic tissue diagnostics based on spectroscopic absorption, scattering, and phase function properties**,” Phys. Med. Biol., 42 803 –814 (1997). https://doi.org/10.1088/0031-9155/42/5/005 0031-9155 Google Scholar**

*Ultraviolet and visible spectroscopies for tissue diagnostics: fluorescence spectroscopy and elastic-scattering spectroscopy***,” Technol. Cancer Res. Treat., 2 (6), 505 –514 (2003). 1533-0346 Google Scholar**

*Instrumentation for multi-modal spectroscopic diagnosis of epithelial dysplasia***,” J. Biomed. Opt., 10 (2), 1 –11 (2005). https://doi.org/10.1117/1.1899686 1083-3668 Google Scholar**

*Combined reflectance and fluorescence spectroscopy for**in vivo*detection of cervical pre-cancer**,” J. Biomed. Opt., 11 (5), 1 –11 (2006). 1083-3668 Google Scholar**

*Autofluorescence spectroscopy of epithelial tissues***,” J. Cell Pharmacol., 3 48 –61 (1992). Google Scholar**

*In vivo*fluorescence in clinical oncology: fundamental and practical issues**,” Encyclopedia of Analytical Chemistry, 20 –56 2000). Google Scholar**

*Fluorescence spectroscopy in vivo***,” Opt. Eng., 34 (11), 3334 –3343 (1995). https://doi.org/10.1117/12.212917 0091-3286 Google Scholar**

*Laser-induced fluorescence spectroscopy at endoscopy: tissue optics, monte carlo modelling, and in vivo measurements***,” J. Biomed. Opt., 6 (4), 385 –396 (2001). https://doi.org/10.1117/1.1413209 1083-3668 Google Scholar**

*Understanding the contributions of NADH and collagen to cervical tissue fluorescence spectra: Modeling, measurements, and implications***,” Cancer, 97 (7), 1681 –1692 (2003). https://doi.org/10.1002/cncr.11255 0008-543X Google Scholar**

*Spectroscopic detection and evaluation of morphologic and biochemical changes in early human oral carcinoma***,” Opt. Express, 11 (24), 3320 –3331 (2003). 1094-4087 Google Scholar**

*Recovery of turbidity free fluorescence from measured fluorescence: an experimental approach***,” Photochem. Photobiol., 68 (5), 603 –632 (1998). 0031-8655 Google Scholar**

*In vivo*fluorescence spectroscopy and imaging for oncological applications**,” Handbook of Optical Biomedical Diagnostics, SPIE Press, Bellingham, WA (2002). Google Scholar**

*Pulse and frequency-domain techniques for tissue spectroscopy and imaging***,” Reviews in Fluorescence 2006, 3 359 –386 Springer, New York (2006). Google Scholar**

*Spectroscopy for the assessment of melanomas***,” Photochem. Photobiol., 73 (2), 178 –183 (2001). https://doi.org/10.1562/0031-8655(2001)073<0178:IVFSON>2.0.CO;2 0031-8655 Google Scholar**

*In vivo*fluorescence spectroscopy of nonmelanoma skin cancer**,” Gynecol. Oncol., 99 (3), S116 –S120 (2005). https://doi.org/10.1016/j.ygyno.2005.07.054 0090-8258 Google Scholar**

*Diffuse reflectance patterns in cervical spectroscopy***,” J. Biomed. Opt., 9 (5), 940 –950 (2004). https://doi.org/10.1117/1.1782611 1083-3668 Google Scholar**

*Clinical study for classification of benign, dysplastic, and malignant oral lesions using autofluorescence spectroscopy***,” Oral Oncol., 41 (2), 117 –31 (2005). https://doi.org/10.1016/j.oraloncology.2004.07.007 0964-1955 Google Scholar**

*The status of**in vivo*autofluorescence spectroscopy and imaging for oral oncology: review**,” Int. J. Cancer, 104 477 –481 (2003). https://doi.org/10.1002/ijc.10959 0020-7136 Google Scholar**

*Optimal excitation-emission wavelengths for autofluorescence diagnosis of bladder tumors***,” Lasers Surg. Med., 38 714 –724 (2006). https://doi.org/10.1002/lsm.20356 0196-8092 Google Scholar**

*Diagnosis of breast cancer using diffuse reflectance spectroscopy: comparison of a Monte Carlo versus partial least squares analysis based feature extraction technique***,” Am. J. Obstet. Gynecol., 186 (3), 374 –382 (2002). https://doi.org/10.1067/mob.2002.121075 0002-9378 Google Scholar**

*Tri-modal spectroscopy for the detection and characterization of cervical precancers**in vivo***,” Ann. Surg. Oncol., 11 (1), 65 –70 (2004). https://doi.org/10.1007/BF02524348 1068-9265 Google Scholar**

*Autofluorescence and diffuse reflectance properties of malignant and benign breast tissues***,” Am. College Chest Phys., 129 (4), 995 –1001 (2006). Google Scholar**

*Optical spectroscopy for the classification of malignant lesions of the bronchial tree***,” Lasers Surg. Med., 36 (5), 356 –364 (2005). https://doi.org/10.1002/lsm.20122 0196-8092 Google Scholar**

*Autofluorescence and diffuse reflectance spectroscopy for oral oncology***,” Appl. Opt., 40 (25), 4633 –4646 (2001). https://doi.org/10.1364/AO.40.004633 0003-6935 Google Scholar**

*Intrinsic fluorescence spectroscopy in turbid media: disentangling effects of scattering and absorption***,” Comput. Methods Programs Biomed., 70 179 –186 (2003). https://doi.org/10.1016/S0169-2607(02)00099-8 0169-2607 Google Scholar**

*Computer simulation of the skin reflectance spectra***,” Lasers Surg. Med., 37 (5), 382 –395 (2005). https://doi.org/10.1002/lsm.20248 0196-8092 Google Scholar**

*Autofluorescence characterization for the early diagnosis of neoplastic changes in DMBA/TPA-induced mouse skin carcinogenesis***,” Appl. Opt., 46 (10), 1686 –1692 (2007). https://doi.org/10.1364/AO.46.001686 0003-6935 Google Scholar**

*Comparison of simplified Monte Carlo simulation and diffusion approximation for the fluorescence signal from phantoms with typical mouse tissue optical properties***,” J. Biomed. Opt., 8 223 –236 (2003). https://doi.org/10.1117/1.1559057 1083-3668 Google Scholar**

*Experimental validation of Monte Carlo modeling of fluorescence in tissues in the UV-visible spectrum***,” Comput. Methods Programs Biomed., 47 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar**

*MCML—Monte Carlo modeling of light transport in multi-layered tissues***,” Med. Phys., 10 824 –830 (1983). https://doi.org/10.1118/1.595361 0094-2405 Google Scholar**

*A Monte Carlo model for the absorption and flux distribution of light in tissue***,” Lasers Surg. Med., 12 510 –519 (1992). https://doi.org/10.1002/lsm.1900120510 0196-8092 Google Scholar**

*Optical properties of intralipid: a phantom medium for light propagation studies***,” Appl. Opt., 32 (19), 3585 –3595 (1993). https://doi.org/10.1364/AO.32.003585 0003-6935 Google Scholar**

*Analytical model for extracting intrinsic fluorescence in turbid media***,” Phys. Med. Biol., 47 (23), 4271 –4285 (2002). https://doi.org/10.1088/0031-9155/47/23/312 0031-9155 Google Scholar**

*Influence of refractive index matching on the photon diffuse reflectance***,” Appl. Opt., 45 (5), 1062 –1071 (2006). https://doi.org/10.1364/AO.45.001062 0003-6935 Google Scholar**

*Monte Carlo-based inverse model for calculating tissue optical properties. Part I: theory and validation on synthetic phantoms***,” Appl. Opt., 45 (5), 1072 –1078 (2006). https://doi.org/10.1364/AO.45.001072 0003-6935 Google Scholar**

*Monte Carlo-based inverse model for calculating tissue optical properties. Part II: application to breast cancer diagnosis***,” IEEE J. Sel. Top. Quantum Electron., 2 (4), 959 –964 (1996). https://doi.org/10.1109/2944.577322 1077-260X Google Scholar**

*Fluorophore quantitation in tissue-simulating media with confocal detection***,” Appl. Opt., 28 (20), 4289 –4292 (1989). https://doi.org/10.1364/AO.28.004286 0003-6935 Google Scholar**

*Fluorescence spectroscopy of turbid media: autofluorescence of the human aorta***,” Lasers Surg. Med., 21 166 –178 (1997). https://doi.org/10.1002/(SICI)1096-9101(1997)21:2<166::AID-LSM8>3.0.CO;2-O 0196-8092 Google Scholar**

*Propagation of fluorescent light***,” IEEE Trans. Biomed. Eng., 43 (2), 113 –122 (1996). https://doi.org/10.1109/10.481980 0018-9294 Google Scholar**

*Morphological model of human colon tissue fluorescence***,” Opt. Express, 13 (19), 7466 –7482 (2005). https://doi.org/10.1364/OPEX.13.007466 1094-4087 Google Scholar**

*Time-resolved photon migration in bi-layered tissue models***,” Opt. Express, 14 (13), 6157 –6171 (2006). https://doi.org/10.1364/OE.14.006157 1094-4087 Google Scholar**

*Quantitative molecular sensing in biological tissues: an approach to non-invasive optical characterization***,” J. Phys. D: Appl. Phys., 36 (14), 1722 –1728 (2003). https://doi.org/10.1088/0022-3727/36/14/311 0022-3727 Google Scholar**

*Analysis of skin fluorescence distribution by the Monte Carlo simulation***,” J. Biomed. Opt., 9 (2), 339 –346 (2004). https://doi.org/10.1117/1.1645796 1083-3668 Google Scholar**

*Amending of fluorescence sensor signal localization in human skin by matching of the refractive index***,” J. Biomed. Opt., 13 (2), 024017 (2008). https://doi.org/10.1117/1.2907161 1083-3668 Google Scholar**

*Monte-carlo-based model for the extraction of intrinsic fluorescence from turbid media***,” J. Opt. Soc. Am. A, 20 (4), 714 –727 (2003). https://doi.org/10.1364/JOSAA.20.000714 0740-3232 Google Scholar**

*Accelerated monte carlo models to simulate fluorescence spectra from layered tissues***,” Appl. Opt., 32 (4), 426 –434 (1993). https://doi.org/10.1364/AO.32.000426 0003-6935 Google Scholar**

*Condensed monte carlo simulations for the description of light transport***,” Opt. Express, 16 (17), 13188 –13202 (2008). https://doi.org/10.1364/OE.16.013188 1094-4087 Google Scholar**

*Monte carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media***,” 102 –111 (1989). Google Scholar**

*A Monte Carlo model of light propagation in tissue***,” Photochem. Photobiol., 67 (1), 23 –62 (1998). https://doi.org/10.1111/j.1751-1097.1998.tb05161.x 0031-8655 Google Scholar**

*Light distributions from point, line and plane sources for photochemical reactions and fluorescence in turbid biological tissues***,” Astrophys. J., 93 70 –83 (1941). https://doi.org/10.1086/144246 0004-637X Google Scholar**

*Diffuse radiation in the galaxy***,” J. Biomed. Opt., 9 (3), 511 –522 (2004). https://doi.org/10.1117/1.1695559 1083-3668 Google Scholar**

*Analytical model to describe fluorescence spectra of normal and preneoplastic epithelial tissue: comparison with Monte Carlo simulations and clinical measurements***,” Polish J. Med. Phys. Eng., 10 (4), 209 –221 (2004). Google Scholar**

*Fluorescence spectra of some photosensitizers in solutions and in tissue-like media. Experiment and Monte Carlo simulation***,” Appl. Opt., 30 4507 –4514 (1991). https://doi.org/10.1364/AO.30.004507 0003-6935 Google Scholar**

*Light scattering in Intralipid-10% in the wavelength range of 400–1100 nanometers***,” Appl. Opt., 42 (21), 4375 –4380 (2003). https://doi.org/10.1364/AO.42.004375 0003-6935 Google Scholar**

*Optical properties of high-density dispersions of particles: application to intralipid solutions***,” Phys. Med. Biol., 47 3387 –3405 (2002). https://doi.org/10.1088/0031-9155/47/18/308 0031-9155 Google Scholar**

*Quantitative fluorescence lifetime spectroscopy in turbid media: comparison of theoretical, experimental and computational methods***,” Sens. Actuators B, 105 365 –377 (2004). https://doi.org/10.1016/j.snb.2004.06.025 0925-4005 Google Scholar**

*A Monte Carlo simulation of photon propagation in fluorescent poly(ethylene glycol)hydrogel microsensors***,” J. Mod. Opt., 2 (4), 153 –162 (1955). https://doi.org/10.1080/713821040 0950-0340 Google Scholar**

*Reflection by semi-infinite diffusers*